text:
!
!     file hwsplr.f90
!
!     * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
!     *                                                               *
!     *                  copyright (c) 2005 by UCAR                   *
!     *                                                               *
!     *       University Corporation for Atmospheric Research         *
!     *                                                               *
!     *                      all rights reserved                      *
!     *                                                               *
!     *                         Fishpack                              *
!     *                                                               *
!     *                      A Package of Fortran                     *
!     *                                                               *
!     *                Subroutines and Example Programs               *
!     *                                                               *
!     *               for Modeling Geophysical Processes              *
!     *                                                               *
!     *                             by                                *
!     *                                                               *
!     *        John Adams, Paul Swarztrauber and Roland Sweet         *
!     *                                                               *
!     *                             of                                *
!     *                                                               *
!     *         the National Center for Atmospheric Research          *
!     *                                                               *
!     *                Boulder, Colorado  (80307)  U.S.A.             *
!     *                                                               *
!     *                   which is sponsored by                       *
!     *                                                               *
!     *              the National Science Foundation                  *
!     *                                                               *
!     * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
!
!     SUBROUTINE hwsplr(a, b, m, mbdcnd, bda, bdb, c, d, n, nbdcnd, bdc, bdd,
!                       elmbda, f, idimf, pertrb, ierror)
!
!
! DIMENSION OF           bda(n), bdb(n), bdc(m), bdd(m), f(idimf, n+1)
! ARGUMENTS
!
! LATEST REVISION        May 2016
!
! PURPOSE                Solves a finite difference approximation to
!                        the helmholtz equation in polar coordinates.
!                        the equation is
!
!                            (1/r)(d/dr)(r(du/dr)) +
!                            (1/r**2)(d/dtheta)(du/dtheta) +
!                            lambda*u = f(r, theta).
!
! USAGE                  call hwsplr(a, b, m, mbdcnd, bda, bdb, c, d, n,
!                                     nbdcnd, bdc, bdd, elmbda, f, idimf,
!                                     pertrb, ierror, w)
!
! ARGUMENTS
! ON INPUT               a, b
!                          the range of r, i.e., a <= r <= b.
!                          a must be less than b and a must be
!                          non-negative.
!
!                        m
!                          the number of panels into which the
!                          interval (a, b) is subdivided.  hence,
!                          there will be m+1 grid points in the
!                          r-direction given by r(i) = a+(i-1)dr,
!                          for i = 1, 2, ..., m+1,
!                          where dr = (b-a)/m is the panel width.
!                          m must be greater than 3.
!
!                        mbdcnd
!                          indicates the type of boundary condition
!                          at r = a and r = b.
!
!                          = 1  if the solution is specified at
!                               r = a and r = b.
!                          = 2  if the solution is specified at
!                               r = a and the derivative of
!                               the solution with respect to r is
!                               specified at r = b.
!                          = 3  if the derivative of the solution
!                               with respect to r is specified at
!                               r = a (see note below) and r = b.
!                          = 4  if the derivative of the solution
!                               with respect to r is specified at
!                               r = a (see note below) and the
!                               solution is specified at r = b.
!                          = 5  if the solution is unspecified at
!                               r = a = 0 and the solution is
!                               specified at r = b.
!                          = 6  if the solution is unspecified at
!                               r = a = 0 and the derivative of the
!                               solution with respect to r is specified
!                               at r = b.
!
!                          note:
!                          if a = 0, do not use mbdcnd = 3 or 4, but
!                          instead use mbdcnd = 1, 2, 5, or 6  .
!
!                        bda
!                          a one-dimensional array of length n+1 that
!                          specifies the values of the derivative of
!                          the solution with respect to r at r = a.
!
!                          when mbdcnd = 3 or 4,
!                            bda(j) = (d/dr)u(a, theta(j)),
!                            j = 1, 2, ..., n+1  .
!
!                          when mbdcnd has any other value, bda is
!                          a dummy variable.
!
!                        bdb
!                          a one-dimensional array of length n+1 that
!                          specifies the values of the derivative of
!                          the solution with respect to r at r = b.
!
!                          when mbdcnd = 2, 3, or 6,
!                            bdb(j) = (d/dr)u(b, theta(j)),
!                            j = 1, 2, ..., n+1  .
!
!                          when mbdcnd has any other value, bdb is
!                          a dummy variable.
!
!                        c, d
!                          the range of theta, i.e., c <=
!                          theta <= d.  c must be less than d.
!
!                        n
!                          the number of panels into which the
!                          interval (c, d) is subdivided.  hence,
!                          there will be n+1 grid points in the
!                          theta-direction given by
!                          theta(j) = c+(j-1)dtheta for
!                          j = 1, 2, ..., n+1, where
!                          dtheta = (d-c)/n is the panel width.
!                          n must be greater than 3.
!
!                        nbdcnd
!                          indicates the type of boundary conditions
!                          at theta = c and at theta = d.
!
!                          = 0  if the solution is periodic in theta,
!                               i.e., u(i, j) = u(i, n+j).
!                          = 1  if the solution is specified at
!                               theta = c and theta = d
!                               (see note below).
!                          = 2  if the solution is specified at
!                               theta = c and the derivative of the
!                               solution with respect to theta is
!                               specified at theta = d
!                               (see note below).
!                          = 4  if the derivative of the solution
!                               with respect to theta is specified
!                               at theta = c and the solution is
!                               specified at theta = d
!                               (see note below).
!
!                          note:
!                          when nbdcnd = 1, 2, or 4, do not use
!                          mbdcnd = 5 or 6
!                          (the former indicates that the solution
!                          is specified at r = 0, the latter indicates
!                          the solution is unspecified at r = 0).
!                          use instead mbdcnd = 1 or 2  .
!
!                        bdc
!                          a one-dimensional array of length m+1 that
!                          specifies the values of the derivative
!                          of the solution with respect to theta at
!                          theta = c.  when nbdcnd = 3 or 4,
!
!                            bdc(i) = (d/dtheta)u(r(i), c),
!                            i = 1, 2, ..., m+1  .
!
!                          when nbdcnd has any other value, bdc is
!                          a dummy variable.
!
!                        bdd
!                          a one-dimensional array of length m+1 that
!                          specifies the values of the derivative
!                          of the solution with respect to theta at
!                          theta = d.  when nbdcnd = 2 or 3,
!
!                            bdd(i) = (d/dtheta)u(r(i), d),
!                            i = 1, 2, ..., m+1  .
!
!                          when nbdcnd has any other value, bdd is
!                          a dummy variable.
!
!                        elmbda
!                          the constant lambda in the helmholtz
!                          equation.  if lambda < 0, a solution
!                          may not exist.  however, hwsplr will
!                          attempt to find a solution.
!
!                        f
!                          a two-dimensional array, of dimension at
!                          least (m+1)*(n+1), specifying values
!                          of the right side of the helmholtz
!                          equation and boundary data (if any).
!
!                          on the interior, f is defined as follows:
!                          for i = 2, 3, ..., m and j = 2, 3, ..., n
!                          f(i, j) = f(r(i), theta(j)).
!
!                          on the boundaries f is defined as follows:
!                          for j = 1, 2, ..., n+1 and i = 1, 2, ..., m+1
!
!                          mbdcnd   f(1, j)            f(m+1, j)
!                          ------   -------------     -------------
!
!                            1      u(a, theta(j))     u(b, theta(j))
!                            2      u(a, theta(j))     f(b, theta(j))
!                            3      f(a, theta(j))     f(b, theta(j))
!                            4      f(a, theta(j))     u(b, theta(j))
!                            5      f(0, 0)            u(b, theta(j))
!                            6      f(0, 0)            f(b, theta(j))
!
!                          nbdcnd   f(i, 1)            f(i, n+1)
!                          ------   ---------         ---------
!
!                            0      f(r(i), c)         f(r(i), c)
!                            1      u(r(i), c)         u(r(i), d)
!                            2      u(r(i), c)         f(r(i), d)
!                            3      f(r(i), c)         f(r(i), d)
!                            4      f(r(i), c)         u(r(i), d)
!
!                          note:
!                          if the table calls for both the solution
!                          u and the right side f at a corner then
!                          then the solution must be specified.
!
!                        idimf
!                          the row (or first) dimension of the array
!                          f as it appears in the program calling
!                          hwsplr.  this parameter is used to specify
!                          the variable dimension of f.  idimf must
!                          be at least m+1.
!
! ON OUTPUT              f
!                          contains the solution u(i, j) of the finite
!                          difference approximation for the grid point
!                          (r(i), theta(j)),
!                          i = 1, 2, ..., m+1, j = 1, 2, ..., n+1  .
!
!                        pertrb
!                          if a combination of periodic, derivative,
!                          or unspecified boundary conditions is
!                          specified for a poisson equation
!                          (lambda = 0), a solution may not exist.
!                          pertrb is a constant, calculated and
!                          subtracted from f, which ensures that a
!                          solution exists.  hwsplr then computes
!                          this solution, which is a least squares
!                          solution to the original approximation.
!                          this solution plus any constant is also
!                          a solution.  hence, the solution is not
!                          unique.  pertrb should be small compared
!                          to the right side. otherwise, a solution
!                          is obtained to an essentially different
!                          problem.  this comparison should always
!                          be made to insure that a meaningful
!                          solution has been obtained.
!
!                        ierror
!                          An error flag that indicates invalid input
!                          parameters. Except for numbers 0 and 11,
!                          a solution is not attempted.
!
!                          =  0  no error.
!                          =  1  a < 0  .
!                          =  2  a >= b.
!                          =  3  mbdcnd < 1 or mbdcnd > 6  .
!                          =  4  c >= d.
!                          =  5  n <= 3
!                          =  6  nbdcnd < 0 or > 4  .
!                          =  7  a = 0, mbdcnd = 3 or 4  .
!                          =  8  a > 0, mbdcnd >= 5  .
!                          =  9  mbdcnd >= 5, nbdcnd /= 0
!                                and nbdcnd /= 3  .
!                          = 10  idimf < m+1  .
!                          = 11  lambda > 0  .
!                          = 12  m <= 3
!                          = 20 If the dynamic allocation of real and
!                               complex workspace required for solution
!                               fails (for example if n, m are too large
!                               for your computer)
!
!                          Since this is the only means of indicating
!                          a possibly incorrect call to hwsplr, the
!                          user should test ierror after the call.
!
! SPECIAL CONDITIONS     None
!
! I/O                    None
!
! PRECISION              64-bit double precision
!
! REQUIRED files         type_FishpackWorkspace.f90, genbun.f90, type_CyclicReductionUtility.f9090
!
! STANDARD               Fortran 2008
!
! HISTORY                * Written by Roland Sweet at NCAR in the late
!                          1970's. Released on NCAR's public software
!                          libraries in January 1980.
!                        * Revised in June 2004 by John Adams using
!                          Fortran 90 dynamically allocated workspace.
!
! ALGORITHM              The routine defines the finite difference
!                        equations, incorporates boundary data, and
!                        adjusts the right side of singular systems
!                        and then calls genbun to solve the system.
!
! TIMING                 For large m and n, the operation count
!                        is roughly proportional to
!
!                          m*n*log2(n)
!
!                        but also depends on input parameters nbdcnd
!                        and mbdcnd.
!
! ACCURACY               The solution process employed results in a loss
!                        of no more than three significant digits for n
!                        and m as large as 64. More details about
!                        accuracy can be found in the documentation for
!                        subroutine genbun which is the routine that
!                        solves the finite difference equations.
!
! REFERENCES             Swarztrauber, P. and R. Sweet, "Efficient
!                        FORTRAN subprograms for the solution of
!                        elliptic equations"
!                          NCAR TN/IA-109, July, 1975, 138 pp.
!
submodule (centered_helmholtz_solvers) centered_polar_solver

contains

    module subroutine hwsplr(a, b, m, mbdcnd, bda, bdb, c, d, n, nbdcnd, bdc, &
        bdd, elmbda, f, idimf, pertrb, ierror)

        ! Dummy arguments
        integer(ip), intent(in)     :: m
        integer(ip), intent(in)     :: mbdcnd
        integer(ip), intent(in)     :: n
        integer(ip), intent(in)     :: nbdcnd
        integer(ip), intent(in)     :: idimf
        integer(ip), intent(out)    :: ierror
        real(wp),    intent(in)     :: a
        real(wp),    intent(in)     :: b
        real(wp),    intent(in)     :: c
        real(wp),    intent(in)     :: d
        real(wp),    intent(in)     :: elmbda
        real(wp),    intent(out)    :: pertrb
        real(wp),    intent(in)     :: bda(:)
        real(wp),    intent(in)     :: bdb(:)
        real(wp),    intent(in)     :: bdc(:)
        real(wp),    intent(in)     :: bdd(:)
        real(wp),    intent(inout) :: f(:,:)

        ! Local variables
        type(FishpackWorkspace) workspace

        ! Check input arguments
        call hwsplr_check_input_arguments(a, b, m, mbdcnd, c, d, n, nbdcnd, idimf, ierror)

        ! Check error flag
        if (ierror /= 0) return

        ! Allocate memory
        call workspace%initialize_centered_workspace(n, m)

        ! Solve system
        associate( rew => workspace%real_workspace )
            call hwsplr_lower_routine(a, b, m, mbdcnd, bda, bdb, c, d, n, &
                nbdcnd, bdc, bdd, elmbda, f, idimf, pertrb, ierror, rew)
        end associate

        ! Release memory
        call workspace%destroy()

    end subroutine hwsplr

    pure subroutine hwsplr_check_input_arguments(a, b, m, mbdcnd, c, d, n, nbdcnd, idimf, ierror)

        ! Dummy arguments
        integer(ip), intent(in)     :: m
        integer(ip), intent(in)     :: mbdcnd
        integer(ip), intent(in)     :: n
        integer(ip), intent(in)     :: nbdcnd
        integer(ip), intent(in)     :: idimf
        integer(ip), intent(out)    :: ierror
        real(wp),    intent(in)     :: a
        real(wp),    intent(in)     :: b
        real(wp),    intent(in)     :: c
        real(wp),    intent(in)     :: d

        if (a < ZERO) then
            ierror = 1
        else if (a >= b) then
            ierror = 2
        else if (mbdcnd <= 0 .or. mbdcnd >= 7) then
            ierror = 3
        else if (d <= c) then
            ierror = 4
        else if (n <= 3) then
            ierror = 5
        else if (nbdcnd <= -1 .or. 5 <= nbdcnd) then
            ierror = 6
        else if (a == ZERO .and. (mbdcnd==3 .or. mbdcnd==4)) then
            ierror = 7
        else if (a > ZERO .and. 5 <= mbdcnd) then
            ierror = 8
        else if (5 <= mbdcnd .and. nbdcnd /= 0 .and. nbdcnd /= 3) then
            ierror = 9
        else if (idimf < m + 1) then
            ierror = 10
        else if (m <= 3) then
            ierror = 12
        else
            ierror = 0
        end if

    end subroutine hwsplr_check_input_arguments

    subroutine hwsplr_lower_routine(a, b, m, mbdcnd, bda, bdb, c, d, n, nbdcnd, bdc, &
        bdd, elmbda, f, idimf, pertrb, ierror, w)

        ! Dummy arguments
        integer(ip), intent(in)     :: m
        integer(ip), intent(in)     :: mbdcnd
        integer(ip), intent(in)     :: n
        integer(ip), intent(in)     :: nbdcnd
        integer(ip), intent(in)     :: idimf
        integer(ip), intent(out)    :: ierror
        real(wp),    intent(in)     :: a
        real(wp),    intent(in)     :: b
        real(wp),    intent(in)     :: c
        real(wp),    intent(in)     :: d
        real(wp),    intent(in)     :: elmbda
        real(wp),    intent(out)    :: pertrb
        real(wp),    intent(in)     :: bda(:)
        real(wp),    intent(in)     :: bdb(:)
        real(wp),    intent(in)     :: bdc(:)
        real(wp),    intent(in)     :: bdd(:)
        real(wp),    intent(inout)  :: f(:,:)
        real(wp),    intent(inout)  :: w(:)

        ! Local variables
        integer(ip) :: mp1, np1, np, mstart, mstop, munk, nstart, nstop, nunk
        integer(ip) :: id2, id3, id4, id5, id6, ij, i
        integer(ip) :: j, l, lp, k, i1, local_error_flag, iip
        real(wp)    :: dr, half_dr, dr2, dt, dt2
        real(wp)    :: a1, r, s2, a2, s, s1, ypole
        type(CenteredCyclicReductionUtility) :: util

        mp1 = m + 1
        dr = (b - a)/m
        half_dr = dr/2
        dr2 = dr**2
        np1 = n + 1
        dt = (d - c)/n
        dt2 = dt**2
        np = nbdcnd + 1
        !
        ! Define range of indices i and j for unknowns u(i, j).
        !
        mstart = 2
        mstop = mp1

        select case (mbdcnd)
            case (1)
                mstop = m
            case (3)
                mstart = 1
            case (4)
                mstart = 1
                mstop = m
            case (5)
                mstop = m
        end select

        munk = mstop - mstart + 1
        nstart = 1
        nstop = n

        select case (np)
            case (2)
                nstart = 2
            case (3)
                nstart = 2
                nstop = np1
            case (4)
                nstop = np1
        end select

        nunk = nstop - nstart + 1
        !
        ! Define a, b, c coefficients in w-array.
        !
        id2 = munk
        id3 = id2 + munk
        id4 = id3 + munk
        id5 = id4 + munk
        id6 = id5 + munk
        a1 = TWO/dr2
        ij = 0

        if (mbdcnd == 3 .or. mbdcnd == 4) ij = 1

        do i = 1, munk
            r = a + real(i - ij, kind=wp)*dr
            j = id5 + i
            w(j) = r
            j = id6 + i
            w(j) = ONE/r**2
            w(i) = (r - half_dr)/(r*dr2)
            j = id3 + i
            w(j) = (r + half_dr)/(r*dr2)
            j = id2 + i
            w(j) = (-a1) + elmbda
        end do

        select case (mbdcnd)
            case (2, 6)
                w(id2) = a1
            case (3)
                w(id2) = a1
                w(id3+1) = a1
            case (4)
                w(id3+1) = a1
        end select

        select case (mbdcnd)
            case (1:2)
                a1 = w(1)
                f(2, nstart:nstop) = f(2, nstart:nstop) - a1*f(1, nstart:nstop)
            case (3:4)
                a1 = TWO * dr*w(1)
                f(1, nstart:nstop) = f(1, nstart:nstop) + a1*bda(nstart:nstop)
        end select

        select case (mbdcnd)
            case (1, 4:5)
                a1 = w(id4)
                f(m, nstart:nstop) = f(m, nstart:nstop) - a1*f(mp1, nstart:nstop)
            case (2:3, 6)
                a1 = TWO * dr*w(id4)
                f(mp1, nstart:nstop) = f(mp1, nstart:nstop) - a1*bdb(nstart:nstop)
        end select

        !
        ! Enter boundary data for theta-boundaries.
        !
        a1 = ONE/dt2
        l = id5 - mstart + 1
        lp = id6 - mstart + 1

        if (np /= 1) then
            select case (np)
                case (2:3)
                    f(mstart:mstop, 2) = f(mstart:mstop, 2) - a1*w(mstart+lp:mstop+lp)*f &
                        (mstart:mstop, 1)
                case (4:5)
                    a1 = TWO/dt
                    f(mstart:mstop, 1) = f(mstart:mstop, 1) + a1*w(mstart+lp:mstop+lp)* &
                        bdc(mstart:mstop)
            end select

            a1 = ONE/dt2

            select case (np)
                case (2, 5)
                    f(mstart:mstop, n) = f(mstart:mstop, n) &
                        - a1*w(mstart+lp:mstop+lp) * f(mstart:mstop, np1)
                case (3:4)
                    a1 = TWO/dt
                    f(mstart:mstop, np1) = f(mstart:mstop, np1) &
                        - a1 * w(mstart+lp:mstop+lp) * bdd(mstart:mstop)
            end select
        end if

        if ((mbdcnd >= 5) .and. (nbdcnd == 3)) then
            f(1, 1) = f(1, 1) - (bdd(2)-bdc(2))* 4.0_wp/(real(n, kind=wp)*dt*dr2)
        end if
        !
        !     adjust right side of singular problems to insure existence of a
        !     solution.
        !
        pertrb = ZERO

        if_construct: if (elmbda >= ZERO) then
            if (elmbda /= ZERO) then
                ierror = 11
                return
            else
                if (nbdcnd == 0 .or. nbdcnd == 3) then
                    s2 = ZERO
                    select case (mbdcnd)
                        case (1:2, 4:5)
                            exit if_construct
                        case (3)
                            w(id5+1) = HALF * (w(id5+2)-half_dr)
                            s2 = 0.25_wp * dr
                    end select

                    if (nbdcnd == 0) then
                        a2 = ONE
                    else
                        a2 = TWO
                    end if

                    j = id5 + munk
                    w(j) = HALF * (w(j-1)+half_dr)
                    s = ZERO

                    do i = mstart, mstop
                        s1 = ZERO
                        ij = nstart + 1
                        k = nstop - 1
                        s1 = sum(f(i, ij:k))
                        j = i + l
                        s = s + (a2*s1 + f(i, nstart)+f(i, nstop))*w(j)
                    end do

                    s2=real(m, kind=wp)*a+dr*(real((m-1)*(m+1), kind=wp)*HALF+0.25_wp)+s2
                    s1 = (TWO + a2*real(nunk - 2, kind=wp))*s2

                    if (mbdcnd /= 3) then
                        s2 = (real(n, kind=wp)*a2*dr)/8
                        s = s + f(1, 1)*s2
                        s1 = s1 + s2
                    end if

                    pertrb = s/s1
                    f(mstart:mstop, nstart:nstop) = &
                        f(mstart:mstop, nstart:nstop) - pertrb
                end if
            end if
        end if if_construct

        do i = mstart, mstop
            k = i - mstart + 1
            j = i + lp
            a1 = dt2/w(j)
            w(k) = a1*w(k)
            j = id2 + k
            w(j) = a1*w(j)
            j = id3 + k
            w(j) = a1*w(j)
            f(i, nstart:nstop) = a1*f(i, nstart:nstop)
        end do

        w(1) = ZERO
        w(id4) = ZERO
        !
        ! Solve the system of equations.
        !
        i1 = 1
        local_error_flag = 0
        associate( &
            a_arg => w(1:munk), &
            b_arg => w(id2+1:id2+1+munk), &
            c_arg => w(id3+1:id3+1+munk), &
            y_arg => f(mstart:,nstart:nstart+nunk), &
            w_arg => w(id4+1:) &
            )
            call util%genbun_lower_routine(nbdcnd, nunk, i1, munk, a_arg, b_arg, c_arg, &
                idimf, y_arg, local_error_flag, w_arg)
        end associate

        ! Check error flag
        if (local_error_flag /= 0) then
            error stop 'fishpack library: genbun_lower_routine call failed in hwsplr_lower_routine'
        end if

        select case (mbdcnd)
            case (1:4)
                if (nbdcnd == 0) f(mstart:mstop, np1) = f(mstart:mstop, 1)
            case (5)
                j = id5 + munk
                w(j) = w(id2)/w(id3)

                do iip = 3, munk
                    i = munk - iip + 2
                    j = id5 + i
                    lp = id2 + i
                    k = id3 + i
                    w(j) = w(i)/(w(lp)-w(k)*w(j+1))
                end do

                w(id5+1) = -HALF * dt2/(w(id2+1)-w(id3+1)*w(id5+2))

                do i = 2, munk
                    j = id5 + i
                    w(j) = -w(j)*w(j-1)
                end do

                s = ZERO
                s = sum(f(2, nstart:nstop))
                a2 = nunk

                if (nbdcnd /= 0) then
                    s = s - HALF * (f(2, nstart)+f(2, nstop))
                    a2 = a2 - ONE
                end if

                ypole = (0.25_wp *dr2*f(1, 1)-s/a2)/(w(id5+1)-ONE + elmbda*dr2* 0.25_wp)

                do i = mstart, mstop
                    k = l + i
                    f(i, nstart:nstop) = f(i, nstart:nstop) + ypole*w(k)
                end do

                f(1, :np1) = ypole

                if (nbdcnd == 0) f(mstart:mstop, np1) = f(mstart:mstop, 1)

            case (6)
                !
                ! Adjust the solution as necessary for the problems where a = 0.
                !
                if (elmbda == ZERO) then
                    ypole = ZERO
                    f(1, :np1) = ypole
                    if (nbdcnd == 0) f(mstart:mstop, np1) = f(mstart:mstop, 1)
                end if
        end select

    end subroutine hwsplr_lower_routine

end submodule centered_polar_solver
!
! REVISION HISTORY
!
! September 1973    Version 1
! April     1976    Version 2
! January   1978    Version 3
! December  1979    Version 3.1
! February  1985    Documentation upgrade
! November  1988    Version 3.2, FORTRAN 77 changes
! June      2004    Version 5.0, Fortran 90 changes
! May       2016    Fortran 2008 changes
!

max_stars_repo_name: jlokimlin/fishpack

max_issues_repo_name: jlokimlin/fishpack


text:
      program DRSPLO
c>> 1998-01-15 Krogh Generate the plot used in ch02-02.tex
c--S replaces "?": DR?PLO, DR?PL1, dr?plot, ?PLOT, ?ERF, ?ERFC, ?ERFCE
      integer MAXPT
      real             RMAXPT, XA, XB, XSIZE, YSIZE
      parameter (MAXPT=101, RMAXPT=101.E0)
      parameter (XA = 0.0E0, XB = 5.0E0)
      parameter (XSIZE=3.2E0, YSIZE=2.75E0)
      real             DX, OPT(17), X(MAXPT), Y(MAXPT, 3)
      external SERF, SERFC, SERFCE
      real             SERF, SERFC, SERFCE
      integer I
      character COPT*41
      real             PLONY, PLOTXT
      parameter (PLONY=2.E0, PLOTXT=14.E0)
      data OPT / 0.E0, PLONY, RMAXPT, 3.E0, PLOTXT, 1.1E0, .8E0, 0.E0,
     1   PLOTXT, 1.2E0, .4E0, 0.E0, PLOTXT, 1.3E0, .1E0, 0.E0, 0.E0 /
      data COPT /
     1 'F{drsplot.tex}{erf(x)}{erfce(x)}{erfc(x)}'/
c
      DX = (XB - XA) / (MAXPT - 1)
      do 100 I = 1, MAXPT
         X(I) = real(I - 1) * DX
         Y(I,1) = SERF(X(I))
         Y(I,2) = SERFC(X(I))
         Y(I,3) = SERFCE(X(I))
  100 continue
      call SPLOT(XSIZE, YSIZE, X, MAXPT, Y, OPT, COPT)
c
c Don't use code below this point.  This code is included to make
c possible automated checks of the software in different environments.
c
      call DRSPL1
      stop
      end

      subroutine DRSPL1
c            Print .tex output for checkout purposes.
c++ CODE for ~.C. is active
      character BUF*80
      integer I
      open (unit=10, FILE='drsplot.tex', status='OLD')
   10 BUF = ' '
         read(10, '(A)', END=40) BUF
         do 20 I = 80, 2, -1
            if (BUF(I:I) .ne. ' ') go to 30
   20    continue
   30    print '(1X,A)', BUF(1:I)
      go to 10
   40 continue
      close(10)
c++ CODE for .C. is inactive
c%%   char  buf[81];
c%%   FILE  *iofil;
c%%   iofil = fopen("drsplot.tex", "r");
c%%   while (fscanf( iofil, "%[^\n]", buf ) != EOF){
c%%      printf("%s\n", buf);
c%%      if (fscanf( iofil, "%[\n]", buf) == EOF) break;}
c%%   fclose(iofil);
c++ END
      return
      end

ext: f


!------------------------------------------------------------------------------
! BSD 2-Clause License
!
! Copyright (c) 2017-2019, Science and Technology Facilities Council
! All rights reserved.
!
! Redistribution and use in source and binary forms, with or without
! modification, are permitted provided that the following conditions are met:
!
! * Redistributions of source code must retain the above copyright notice, this
!   list of conditions and the following disclaimer.
!
! * Redistributions in binary form must reproduce the above copyright notice,
!   this list of conditions and the following disclaimer in the documentation
!   and/or other materials provided with the distribution.
!
! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
! AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
! IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
! DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE
! FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
! DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
! SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
! CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
! OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
! OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
!------------------------------------------------------------------------------
! Author: A. R. Porter, STFC Daresbury Laboratory

!> Module for describing all aspects of a field (which exists on some
!! grid).
module field_mod
  use kind_params_mod
  use region_mod
  use halo_mod
  use grid_mod
  use gocean_mod, only: gocean_stop
  use tile_mod
  use parallel_mod, only: decomposition_type
  implicit none

  private

  ! Enumeration of grid-point types on the Arakawa C grid. A
  ! field lives on one of these types.
  integer, public, parameter :: GO_U_POINTS   = 0
  integer, public, parameter :: GO_V_POINTS   = 1
  integer, public, parameter :: GO_T_POINTS   = 2
  integer, public, parameter :: GO_F_POINTS   = 3
  !> A field that lives on all grid-points of the grid
  integer, public, parameter :: GO_ALL_POINTS = 4

  !> Abstract interfaces for C and Fortran subroutines to be implemented by
  ! the infrastructure user to retrieve data from, or push data to, devices
  ! using the desired programming model.
  !
  ! All interfaces include:
  !  - A host and a device pointer specifying where the data is copied from/to.
  !  - startx, starty, nx, and ny integers to specify if copying just a sub-
  !  region of the field is enough (the implementation may copy a larger part).
  !  - A blocking boolean to specify if the data copy operation should have
  !  finished on the routine return or just been dispatched.
  !
  abstract interface
    subroutine read_from_device_c_interface(from, to, startx, starty, nx, ny, blocking)
        use iso_c_binding, only: c_ptr, c_int, c_bool
        type(c_ptr), intent(in), value :: from
        type(c_ptr), intent(in), value :: to
        integer(c_int), intent(in), value :: startx, starty, nx, ny
        logical(c_bool), intent(in), value :: blocking
    end subroutine read_from_device_c_interface
  end interface

  abstract interface
    subroutine read_from_device_f_interface(from, to, startx, starty, nx, ny, blocking)
        use iso_c_binding, only: c_ptr
        use kind_params_mod, only: go_wp
        type(c_ptr), intent(in) :: from
        real(go_wp), dimension(:,:), target, intent(inout) :: to
        integer, intent(in) ::  startx, starty, nx, ny
        logical, intent(in) :: blocking
    end subroutine read_from_device_f_interface
  end interface

  abstract interface
    subroutine write_to_device_c_interface(from, to, startx, starty, nx, ny, blocking)
        use iso_c_binding, only: c_ptr, c_int, c_bool
        type(c_ptr), intent(in), value :: from
        type(c_ptr), intent(in), value :: to
        integer(c_int), intent(in), value :: startx, starty, nx, ny
        logical(c_bool), intent(in), value :: blocking
    end subroutine write_to_device_c_interface
  end interface

  abstract interface
    subroutine write_to_device_f_interface(from, to, startx, starty, nx, ny, blocking)
        use iso_c_binding, only: c_ptr
        use kind_params_mod, only: go_wp
        real(go_wp), dimension(:,:), target, intent(in) :: from
        type(c_ptr), intent(in) :: to
        integer, intent(in) ::  startx, starty, nx, ny
        logical, intent(in) :: blocking
    end subroutine write_to_device_f_interface
  end interface

  !> The base field type. Intended to represent a global field
  !! such as the x-component of velocity.
  type, public :: field_type
     !> Which mesh points the field is defined upon
     integer :: defined_on
     !> The grid on which this field is defined
     type(grid_type), pointer :: grid
     !> The internal region of this field
     type(region_type) :: internal
     !> The whole region covered by this field - includes
     !! boundary points
     type(region_type) :: whole
     !> The number of halo regions that this field has.
     !! Halo region values are not computed but copied
     !! from elsewhere.
     integer :: num_halos
     !> Array of objects describing the halos belonging to this field.
     type(halo_type), dimension(:), allocatable :: halo
     !> Whether the data for this field lives in a remote memory space
     !! (e.g. on a GPU)
     logical :: data_on_device
     !> Function pointers to the functions that read/write data from/to
     ! devices.
     procedure(read_from_device_c_interface), POINTER, nopass :: read_from_device_c
     procedure(read_from_device_f_interface), POINTER, nopass :: read_from_device_f
     procedure(write_to_device_c_interface), POINTER, nopass :: write_to_device_c
     procedure(write_to_device_f_interface), POINTER, nopass :: write_to_device_f

  end type field_type

  !> A real, 2D field.
  type, public, extends(field_type) :: r2d_field
     integer :: ntiles
     !> The dimensions of the tiles into which the field
     !! is sub-divided.
     type(tile_type), dimension(:), allocatable :: tile
     !> Array holding the actual field values. Eventually this should
     !! be private so that users are forced to access it via get_data().
     !! This will allow us to guarantee that the data is valid (e.g. in
     !! the case where values are computed on a separate accelerator
     !! device).
     real(go_wp), dimension(:,:), allocatable :: data
     !> Pointer to corresponding buffer on remote device (if any).
     !! This requires variables that are declared to be of this type
     !! have the 'target' attribute.
     type(c_ptr) :: device_ptr
   contains
     !> Setter for the data associated with this field
     procedure, pass :: set_data
     !> Getter for the data associated with this field. Fetches data
     !! from remote accelerator if necessary.
     procedure, pass :: get_data
     procedure, pass :: read_halo_from_device
     procedure, pass :: write_halo_to_device
     procedure, pass :: read_from_device
     procedure, pass :: write_to_device
     procedure, public :: halo_exchange
  end type r2d_field

  !> Interface for the copy_field operation. Overloaded to take
  !! an array or an r2d_field type.
  !! \todo Remove support for raw arrays from this interface.
  interface copy_field
     module procedure copy_2dfield_array, copy_2dfield_array_patch, &
                      copy_2dfield, copy_2dfield_patch
  end interface copy_field

  ! User-defined constructor for r2d_field type objects
  interface r2d_field
     module procedure r2d_field_constructor
  end interface r2d_field

  !> Interface for the field checksum operation. Overloaded to take either
  !! a field object or a 2D, real(go_wp) array.
  interface field_checksum
     module procedure fld_checksum, array_checksum
  end interface field_checksum

  interface free_field
     module procedure r2d_free_field
  end interface

  public copy_field
  public set_field
  public field_checksum
  public free_field

! Grid points on an Arakawa C grid with NE offset (i.e. the U,V and F pts
! immediately to the North and East of a T point share its grid indices)
! are arranged like so:
!
!v(1,ny)----f(1,ny)---v(i-1,ny)--f(i-1,ny)--v(i,ny)----f(i,ny)--v(nx,ny)---f(nx,ny)
!|          |         |          |          |          |        |          |
!|          |         |          |          |          |        |          |
!T[1,ny]----u(1,ny)---T(i-1,ny)--u(i-1,ny)--T(i,ny)----u(i,ny)--T(nx,ny)---u(nx,ny)
!|          |         |          |          |          |        |          |
!|          |         |          |          |          |        |          |
!v(1,j)-----f(1,j)----v(i-1,j)---f(i-1,j)---v(i,j)-----f(i,j)---v(nx,j)----f(nx,j)
!|          |         |          |          |          |        |          |
!|          |         |          |          |          |        |          |
!T[1,j]-----u(1,j)----T(i-1,j)---u(i-1,j)---T(i,j)-----u(i,j)---T(nx,j)----u(nx,j)
!|          |         |          |          |          |        |          |
!|          |         |          |          |          |        |          |
!v(1,j-1)---f(1,j-1)--v(i-1,j-1)-f(i-1,j-1)-v(i,j-1)---f(i,j-1)-v(nx,j-1)--f(nx,j-1)
!|          |         |          |          |          |        |          |
!|          |         |          |          |          |        |          |
!T[1,j-1]---u(1,j-1)--T(i-1,j-1)-u(i-1,j-1)-T(i,j-1)---u(i,j-1)-T(nx,j-1)--u(nx,j-1)
!|          |         |          |          |          |        |          |
!|          |         |          |          |          |        |          |
!v(1,1)-----f(1,1)----v(i-1,1)---f(i-1,1)---v(i,1)-----f(i,1)---v(nx,1)----f(nx,1)
!|          |         |          |          |          |        |          |
!|          |         |          |          |          |        |          |
!T[1,1]     u(1,1)    T(i-1,1)---u(i-1,1)---T(i,1)-----u(i,1)---T(nx,1)----u(nx,1)

  !> The no. of cols/rows used to define boundary data in the absence
  !! of periodic BCs.
  !! \todo remove this parameter and determine it from
  !! the supplied T mask.
  integer, public, parameter :: NBOUNDARY = 1

  !> Whether or not to 'tile' (sub-divide) field arrays
  logical, public, parameter :: TILED_FIELDS = .FALSE.

  !> The tiling decomposition to use for each field
  type(decomposition_type), save :: decomp

  !> Whether or not the tiling decomposition has been created yet
  logical, save :: tiling_initialised = .false.

contains

  !===================================================

  function r2d_field_constructor(grid,    &
                                 grid_points, &
                                 do_tile) result(self)
    use parallel_mod, only: go_decompose, on_master
!$    use omp_lib, only : omp_get_max_threads
    implicit none
    ! Arguments
    !> Pointer to the grid on which this field lives
    type(grid_type), intent(in), target  :: grid
    !> Which grid-point type the field is defined on
    integer,         intent(in)          :: grid_points
    !> If the field should be tiled among all threads, or if only
    !> a single field should be allocated (which is not currently
    !> supported by PSyclone)
    logical, intent(in), optional :: do_tile
    ! Local declarations
    type(r2d_field), target :: self
    logical :: local_do_tiling = .false.
    integer :: ierr
    character(len=8) :: fld_type
    integer :: ji, jj
    !> The upper bounds actually used to allocate arrays (as opposed
    !! to the limits carried around with the field)
    integer :: upper_x_bound, upper_y_bound
    integer :: itile, nthreads, ntilex, ntiley

    if (present(do_tile)) then
       local_do_tiling = do_tile
    endif

    ! Set this field's grid pointer to point to the grid pointed to
    ! by the supplied grid_ptr argument
    self%grid => grid

    !> The data associated with this device is currently local
    !! to where we're executing
    self%data_on_device = .FALSE.

    !> The function pointers are initialized with NULL
    nullify(self%read_from_device_c)
    nullify(self%read_from_device_f)
    nullify(self%write_to_device_c)
    nullify(self%write_to_device_f)

    ! Set-up the limits of the 'internal' region of this field
    !
    call set_field_bounds(self,fld_type,grid_points)

    ! Dimensions of the grid of tiles.
    if(.not. get_grid_dims(ntilex, ntiley) )then
       ntilex = 1
       ntiley = 1
    end if

    if (local_do_tiling) then
       nthreads = 1
       !> \TODO #38 We don't actually support building dl_esm_inf
       !! with OpenMP enabled!
!$     nthreads = omp_get_max_threads()

       if (.not. tiling_initialised) then
!$        write (*,"(/'Have ',I3,' OpenMP threads available.')") nthreads
          decomp = go_decompose(self%internal%nx, self%internal%ny, &
               nthreads, ntilex, ntiley)
          tiling_initialised = .true.
       end if

       self%ntiles = nthreads
       allocate(self%tile(self%ntiles), Stat=ierr)
       if(ierr /= 0)then
          call gocean_stop('r2d constructor failed to allocate tiling structures')
       end if
       do itile = 1, self%ntiles
          self%tile(itile)%whole = decomp%subdomains(itile)%global
          self%tile(itile)%internal = decomp%subdomains(itile)%internal
       end do
    else
       ! We're not tiling this field so set ntiles to zero
       self%ntiles = 0
    end if

    ! We allocate *all* fields to have the same extent as that
    ! of the grid. This enables the (Cray) compiler
    ! to safely evaluate code within if blocks that are
    ! checking for conditions at the boundary of the domain.
    ! Hence we use self%whole%{x,y}stop + 1...
    upper_x_bound = self%grid%nx
    upper_y_bound = self%grid%ny

    if (on_master()) then
        write(*, "('Allocating ',(A),' field with bounds: (',I1,':',I4, "// &
            "',',I1,':',I4,'), internal region is (',I1,':',I4, ',',I1,':',I4,')' )") &
            TRIM(ADJUSTL(fld_type)), 1, upper_x_bound, 1, upper_y_bound, &
             self%internal%xstart, self%internal%xstop, &
             self%internal%ystart, self%internal%ystop
    endif

    ! Allocating with a lower bound != 1 causes problems whenever
    ! array passed as assumed-shape dummy argument because lower
    ! bounds default to 1 in called unit (unless array declared as
    ! allocatable). However, all loops will be in the generated, middle
    ! layer and the generator knows the array bounds. This may give us the
    ! ability to solve this problem (by passing array bounds to the
    ! kernels).
    allocate(self%data(1:upper_x_bound, 1:upper_y_bound), &
                       Stat=ierr)
    if(ierr /= 0)then
       call gocean_stop('r2d_field_constructor: ERROR: failed to '// &
                        'allocate field')
    end if

    ! Since we're allocating the arrays to be larger than strictly
    ! required we explicitly set all elements to -999 in case the code
    ! does access 'out-of-bounds' elements during speculative
    ! execution. If we're running with OpenMP this also gives
    ! us the opportunity to do a 'first touch' policy to aid with
    ! memory<->thread locality...
    if (self%ntiles > 0) then
!$OMP PARALLEL DO schedule(runtime), default(none), &
!$OMP private(itile,ji,jj), shared(self)
       do itile = 1, self%ntiles
          do jj = self%tile(itile)%whole%ystart, self%tile(itile)%whole%ystop
             do ji = self%tile(itile)%whole%xstart, self%tile(itile)%whole%xstop
                self%data(ji,jj) = 0.0_go_wp
             end do
          end do
       end do
!$OMP END PARALLEL DO
    else
       self%data(:,:) = 0.0_go_wp
    end if
  end function r2d_field_constructor

  !===================================================
  ! Frees the data allocated for this array. Does nothing
  ! if the data was never allocated.
  subroutine r2d_free_field(fld)
    implicit none
    type(r2d_field), intent(inout) :: fld
    ! Arguments
    !> fld: Pointer to the fld which data is to be freed
    if (allocated(fld%data)) then
       deallocate(fld%data)
    end if
  end subroutine r2d_free_field

  !===================================================

  subroutine read_from_device(self, startx, starty, nx, ny, blocking)
    !> Update the host data with the device data.
    ! The optional startx, starty, nx and ny parameters can be provided to
    ! read just from a sub-region of the field.
    ! The blocking optional parameter specifies if the data transfer should
    ! have finished on the return from this routine.

    class(r2d_field), target :: self
    integer, optional, intent(in) :: startx, starty, nx, ny
    logical, optional, intent(in) :: blocking
    integer :: local_startx, local_starty, local_nx, local_ny
    logical :: local_blocking

    if(self%data_on_device)then

       if (present(startx)) then
           local_startx = startx
       else
           local_startx = 1
       endif

       if (present(starty)) then
           local_starty = starty
       else
           local_starty = 1
       endif

       if (present(nx)) then
           local_nx = nx
       else
           local_nx = size(self%data, 1)
       endif

       if (present(ny)) then
           local_ny = ny
       else
           local_ny = size(self%data, 2)
       endif

       if (present(blocking)) then
           local_blocking = blocking
       else
           local_blocking = .true.
       endif

       if(associated(self%read_from_device_c))then
          call self%read_from_device_c(self%device_ptr, C_LOC(self%data), &
                int(local_startx, c_int), int(local_starty, c_int), &
                int(local_nx, c_int), int(local_ny, c_int), &
                logical(local_blocking, c_bool))
        else if(associated(self%read_from_device_f))then
          call self%read_from_device_f(self%device_ptr, self%data, &
                local_startx, local_starty, local_nx, local_ny, local_blocking)
        else
          call gocean_stop("ERROR: Data is on a device but no instructions " // &
              "about how to retrieve the data have been provided.")
       endif
    end if
  end subroutine read_from_device

  subroutine write_to_device(self, startx, starty, nx, ny, blocking)
    !> Update the device data with host data.
    ! The optional startx, starty, nx and ny parameters can be provided to
    ! write just a sub-region of the field to the device.
    ! The blocking optional parameter specifies if the data transfer should
    ! have finished on the return from this routine.

    class(r2d_field), target :: self
    integer, optional, intent(in) :: startx, starty, nx, ny
    logical, optional, intent(in) :: blocking
    integer :: local_startx, local_starty, local_nx, local_ny
    logical :: local_blocking

    if(self%data_on_device)then

       if (present(startx)) then
           local_startx = startx
       else
           local_startx = 1
       endif

       if (present(starty)) then
           local_starty = starty
       else
           local_starty = 1
       endif

       if (present(nx)) then
           local_nx = nx
       else
           local_nx = size(self%data, 1)
       endif

       if (present(ny)) then
           local_ny = ny
       else
           local_ny = size(self%data, 2)
       endif

       if (present(blocking)) then
           local_blocking = blocking
       else
           local_blocking = .true.
       endif

       if(associated(self%write_to_device_c))then
          call self%write_to_device_c(C_LOC(self%data), self%device_ptr, &
                int(local_startx, c_int), int(local_starty, c_int), &
                int(local_nx, c_int), int(local_ny, c_int), &
                logical(local_blocking, c_bool))
       else if(associated(self%write_to_device_f))then
          call self%write_to_device_f(self%data, self%device_ptr, &
                local_startx, local_starty, local_nx, local_ny, local_blocking)
        else
          call gocean_stop("ERROR: Data is on a device but no instructions " // &
              "about how to write new data have been provided.")
       endif
    endif
  end subroutine write_to_device

  function get_data(self) result(dptr)
    !> Getter for the data associated with a field. If the data is on a
    ! device it ensures that the host copy is up-to-date with that on
    ! any accelerator device.
    class(r2d_field), target :: self
    real(go_wp), dimension(:,:), pointer :: dptr

    ! Synchronise the whole array
    call self%read_from_device()

    ! Return a reference to the data
    dptr => self%data
  end function get_data

  !===================================================

  function set_data(self, array) result(flag)
    !> Setter for the data associated with a field. If data is on a
    !! remote accelerator device then the device copy is updated too.
    implicit none
    class(r2d_field) :: self
    integer :: flag
    real(go_wp), dimension(:,:) :: array
    self%data = array

    ! Synchronise the whole array
    call self%write_to_device()

    flag = 0
  end function set_data

  !===================================================

  subroutine set_field_bounds(fld, fld_type, grid_points)
    implicit none
    !> The field we're working on. We use class as we want this
    ! to be polymorphic as this routine doesn't care whether the
    ! field is tiled or not.
    class(field_type), intent(inout) :: fld
    !> Which grid-point type the field is defined on
    integer,          intent(in)    :: grid_points
    !> Character string describing the grid-point type
    character(len=8), intent(out)   :: fld_type

    select case(grid_points)

    case(GO_U_POINTS)
       write(fld_type, "('C-U')")
       call cu_field_init(fld)
    case(GO_V_POINTS)
       write(fld_type, "('C-V')")
       call cv_field_init(fld)
    case(GO_T_POINTS)
       write(fld_type, "('C-T')")
       call ct_field_init(fld)
    case(GO_F_POINTS)
       write(fld_type, "('C-F')")
       call cf_field_init(fld)
    case(GO_ALL_POINTS)
       write(fld_type, "('C-All')")
       call field_init(fld)
    case default
       call gocean_stop('r2d_field_constructor: ERROR: invalid '//&
                        'specifier for type of mesh points')
    end select

    ! Compute and store dimensions of internal region of field
    fld%internal%nx = fld%internal%xstop - fld%internal%xstart + 1
    fld%internal%ny = fld%internal%ystop - fld%internal%ystart + 1

    ! In addition to the 'internal region' of the field, we may have
    ! external points that define B.C.'s or that act as halos. Here
    ! we store the full extent of the field, inclusive of such
    ! points.
    !> \todo Replace the use of NBOUNDARY here with info. computed
    !! from the T-point mask.
    if(fld%grid%boundary_conditions(1) /= GO_BC_PERIODIC)then
       fld%whole%xstart = fld%internal%xstart - NBOUNDARY
       fld%whole%xstop  = fld%internal%xstop  + NBOUNDARY
    else
       fld%whole%xstart = fld%internal%xstart - NBOUNDARY
       fld%whole%xstop  = fld%internal%xstop  + NBOUNDARY
    end if
    if(fld%grid%boundary_conditions(2) /= go_BC_PERIODIC)then
       fld%whole%ystart = fld%internal%ystart - NBOUNDARY
       fld%whole%ystop  = fld%internal%ystop  + NBOUNDARY
    else
       fld%whole%ystart = fld%internal%ystart - NBOUNDARY
       fld%whole%ystop  = fld%internal%ystop  + NBOUNDARY
    end if

    fld%whole%nx = fld%whole%xstop - fld%whole%xstart + 1
    fld%whole%ny = fld%whole%ystop - fld%whole%ystart + 1

  end subroutine set_field_bounds

  !===================================================

  subroutine field_init(fld)
    implicit none
    class(field_type), intent(inout) :: fld
    ! Locals
    integer :: M, N

    fld%defined_on = GO_ALL_POINTS

    M = fld%grid%nx
    N = fld%grid%ny

    ! An 'all points' field is defined upon every point in the grid
    fld%internal%xstart = 1
    fld%internal%xstop  = M
    fld%internal%ystart = 1
    fld%internal%ystop  = N

    ! We have no halo regions
    fld%num_halos = 0

  end subroutine field_init

  !===================================================

  subroutine cu_field_init(fld)
    implicit none
    class(field_type), intent(inout) :: fld

    fld%defined_on = GO_U_POINTS

    select case(fld%grid%offset)

    case(GO_OFFSET_SW)
       call cu_sw_init(fld)

    case(GO_OFFSET_NE)
       call cu_ne_init(fld)

    case default
       call gocean_stop('cu_field_init: ERROR - unsupported grid offset!')

    end select

  end subroutine cu_field_init

  !===================================================

  subroutine cu_sw_init(fld)
    implicit none
    class(field_type), intent(inout) :: fld

    ! Set up a field defined on U points when the grid has
    ! a South-West offset:

    !   vi-1j+1--fij+1---vij+1---fi+1j+1
    !   |        |       |       |
    !   |        |       |       |
    !   Ti-1j----uij-----Tij-----ui+1j
    !   |        |       |       |
    !   |        |       |       |
    !   vi-1j----fij-----vij-----fi+1j
    !   |        |       |       |
    !   |        |       |       |
    !   Ti-1j-1--uij-1---Tij-1---ui+1j-1

    !
    if(fld%grid%boundary_conditions(1) == GO_BC_PERIODIC)then
       ! When implementing periodic boundary conditions, all mesh
       ! point types have the same extents as the grid of T points. We
       ! then have a halo of width grid_mod::HALO_WIDTH_X on either
       ! side of the domain.
       ! When updating a quantity on U points we write to:
       ! (using 'x' to indicate a location that is written and 'a' an
       ! additional point that shallow dispenses with):
       !
       ! i=istart i=istop
       !  o  o  o  o  a
       !  o  x  x  x  a  j=ystop
       !  o  x  x  x  a
       !  o  x  x  x  a  j=ystart
       !  a  a  a  a  a

       fld%internal%xstart = fld%grid%subdomain%internal%xstart
       fld%internal%xstop  = fld%grid%subdomain%internal%xstop
    else
       ! When updating a quantity on U points with this offset convention
       ! we write to (using 'x' to indicate a location that is written,
       !                    'b' a boundary point and
       !                    'o' a point that is external to the domain):
       !
       ! i=1         i=M
       !  o  b  b  b  b   j=N
       !  o  b  x  x  b
       !  o  b  x  x  b
       !  o  b  x  x  b
       !  o  b  b  b  b   j=1
       fld%internal%xstart = fld%grid%subdomain%internal%xstart + 1
       fld%internal%xstop  = fld%grid%subdomain%internal%xstop
    end if

    fld%internal%ystart = fld%grid%subdomain%internal%ystart
    fld%internal%ystop  = fld%grid%subdomain%internal%ystop

    ! When applying periodic (wrap-around) boundary conditions (PBCs)
    ! we must fill the regions marked with 'b' above.
    ! This looks like (using x to indicate a location that is written
    ! first and y a location that is written second):
    !
    !  i=2      i=M
    ! _ y  y  y  y   j=N
    ! /|x  o  o  o
    ! | x  o  o  o
    ! \ x  o  o  o
    !  \x_ o  o  o   j=1
    !   |\______/

    ! In array notation this looks like:
    !
    ! (2    , 1:N-1) = (M  , 1:N-1)
    ! (2:M, N) = (2:M, 1)

    call init_periodic_bc_halos(fld)

  end subroutine cu_sw_init

  !===================================================

  subroutine cu_ne_init(fld)
    implicit none
    class(field_type), intent(inout) :: fld

    ! Set up a field defined on U points when the grid types have
    ! a North-East offset relative to the T point.

    if(fld%grid%boundary_conditions(1) /= GO_BC_PERIODIC)then
       ! If we do not have periodic boundary conditions then we do
       ! not need to allow for boundary points here - they are
       ! already contained within the region defined by T mask.
       ! The T mask has been used to determine the grid%subdomain
       ! which describes the area on the grid that is actually being
       ! modelled (as opposed to having values supplied from B.C.'s etc.)
       fld%internal%xstart = fld%grid%subdomain%internal%xstart
       ! Now that all fields are allocated to be the same size, we always have
       ! as many U-points as T-points
       fld%internal%xstop  = fld%grid%subdomain%internal%xstop
    else
      call gocean_stop('ERROR: cu_ne_init: implement periodic boundary conditions!')
    end if
    if(fld%grid%boundary_conditions(2) /= GO_BC_PERIODIC)then
      fld%internal%ystart = fld%grid%subdomain%internal%ystart
      fld%internal%ystop  = fld%grid%subdomain%internal%ystop
    else
      call gocean_stop('ERROR: cu_ne_init: implement periodic BCs!')
    end if

!> \todo Is this concept of halo definitions useful?
    fld%num_halos = 0

  end subroutine cu_ne_init

  !===================================================

  subroutine cv_field_init(fld)
    implicit none
    class(field_type), intent(inout) :: fld

    fld%defined_on = GO_V_POINTS

    select case(fld%grid%offset)

    case(GO_OFFSET_SW)
       call cv_sw_init(fld)

    case(GO_OFFSET_NE)
       call cv_ne_init(fld)

    case default
       call gocean_stop('cv_field_init: ERROR - unsupported grid offset!')

    end select

  end subroutine cv_field_init

  !===================================================

  subroutine cv_sw_init(fld)
    implicit none
    class(field_type), intent(inout) :: fld

    if(fld%grid%boundary_conditions(2) == GO_BC_PERIODIC)then
       ! When implementing periodic boundary conditions, all
       ! mesh point types have the same extents as the grid of
       ! T points. We then have a halo of width 1 on either side
       ! of the domain.
       fld%internal%xstart = fld%grid%subdomain%internal%xstart
       fld%internal%xstop  = fld%grid%subdomain%internal%xstop

       fld%internal%ystart = fld%grid%subdomain%internal%ystart
       fld%internal%ystop  = fld%grid%subdomain%internal%ystop
    else
       ! When updating a quantity on V points we write to:
       ! (using x to indicate a location that is written):
       !
       ! i=1      i=M
       !  b  b  b  b   j=N
       !  b  x  x  b
       !  b  x  x  b
       !  b  b  b  b
       !  o  o  o  o   j=1
       ! We're not offset from the T points in the x dimension so
       ! we have the same x bounds.
       fld%internal%xstart = fld%grid%subdomain%internal%xstart
       fld%internal%xstop  = fld%grid%subdomain%internal%xstop

       fld%internal%ystart = fld%grid%subdomain%internal%ystart + 1
       fld%internal%ystop  = fld%grid%subdomain%internal%ystop
       call gocean_stop('cv_sw_init: IMPLEMENT non-periodic BCs!')
    endif

    ! When applying periodic (wrap-around) boundary conditions (PBCs)
    ! we must fill the regions marked with 'b' above.
    ! This looks like (using x to indicate a location that is written
    ! first and y a location that is written second):
    !
    !  i=1     i=M
    !  -o  o  o  y   j=N
    ! / o  o  o  y
    ! | o  o  o  y
    ! \ o  o  o  y
    !  \x  x  x _y   j=2
    !    \______/|

    ! In array notation this looks like:
    ! First row = last internal row
    ! field(1:M    ,1:1  ) = field(1:M,N-1:N-1)
    ! Last col = first internal col
    ! field(M:M,1:N) = field(2:2,  1:N)

    call init_periodic_bc_halos(fld)

  end subroutine cv_sw_init

  !===================================================

  subroutine cv_ne_init(fld)
    implicit none
    class(field_type), intent(inout) :: fld

    if(fld%grid%boundary_conditions(1) /= GO_BC_PERIODIC)then
      ! If we do not have periodic boundary conditions then we do
      ! not need to allow for boundary points here - they are
      ! already contained within the region.
      fld%internal%xstart = fld%grid%subdomain%internal%xstart
      fld%internal%xstop  = fld%grid%subdomain%internal%xstop
    else
      call gocean_stop('ERROR: cv_ne_init: implement periodic BCs!')
    end if

    if(fld%grid%boundary_conditions(2) /= GO_BC_PERIODIC)then
       fld%internal%ystart = fld%grid%subdomain%internal%ystart
       ! Since we allocate all fields with the same extents, we always have
       ! as many V points as T points
       fld%internal%ystop  = fld%grid%subdomain%internal%ystop
    else
      call gocean_stop('ERROR: cv_ne_init: implement periodic BCs!')
    end if

  end subroutine cv_ne_init

  !===================================================

  subroutine ct_field_init(fld)
    implicit none
    class(field_type), intent(inout) :: fld

    fld%defined_on = GO_T_POINTS

    select case(fld%grid%offset)

    case(GO_OFFSET_SW)
       call ct_sw_init(fld)

    case(GO_OFFSET_NE)
       call ct_ne_init(fld)

    case default
       call gocean_stop('ct_field_init: ERROR - unsupported grid offset!')

    end select

  end subroutine ct_field_init

  !===================================================

  subroutine ct_sw_init(fld)
    implicit none
    class(field_type), intent(inout) :: fld

    ! When updating a quantity on T points we write to:
    ! (using x to indicate a location that is written):
    !
    ! i=1      i=M
    !  b  b  b  b   j=N
    !  b  x  x  b
    !  b  x  x  b
    !  b  b  b  b   j=1

    fld%internal%xstart = fld%grid%subdomain%internal%xstart
    fld%internal%xstop  = fld%grid%subdomain%internal%xstop
    fld%internal%ystart = fld%grid%subdomain%internal%ystart
    fld%internal%ystop  = fld%grid%subdomain%internal%ystop

    ! When applying periodic (wrap-around) boundary conditions
    ! (PBCs) we must fill the regions marked with 'b' above.
    ! This looks like (using x to indicate a location that is
    ! written first and y a location that is written second):
    !
    !  i=1      i=M
    ! _ y  y  y  y   j=N
    ! /|o  o  o  x
    ! | o  o  o  x
    ! \ o  o  o  x
    !  \o  o  o _x   j=1
    !    \______/|

    ! In array notation this looks like:
    ! Last col = first col
    ! field(M:M,  1:N-1  ) = field(1:1  ,1:N-1)
    ! Last row = first row
    ! field(1:M,N:N) = field(1:M,1:1)

    call init_periodic_bc_halos(fld)

  end subroutine ct_sw_init

  !===================================================

  subroutine ct_ne_init(fld)
    implicit none
    class(field_type), intent(inout) :: fld

    ! When updating a quantity on T points with a NE offset
    ! we write to (using x to indicate a location that is written):
    !
    ! i=1          Nx
    !  b  b  b  b  b Ny
    !  b  x  x  x  b
    !  b  x  x  x  b
    !  b  x  x  x  b
    !  b  b  b  b  b j=1

    if(fld%grid%boundary_conditions(1) /= GO_BC_PERIODIC)then
      ! If we do not have periodic boundary conditions then we do
      ! not need to allow for boundary points here - they are
      ! already contained within the region.
      ! Start and stop are just the same as those calculated from the T mask
      ! earlier because this is a field on T points.
      fld%internal%xstart = fld%grid%subdomain%internal%xstart
      fld%internal%xstop  = fld%grid%subdomain%internal%xstop
    else
      call gocean_stop('ERROR: ct_ne_init: implement periodic BCs!')
    end if

    if(fld%grid%boundary_conditions(2) /= GO_BC_PERIODIC)then
      ! Start and stop are just the same as those calculated from the T mask
      ! earlier because this is a field on T points.
      fld%internal%ystart = fld%grid%subdomain%internal%ystart
      fld%internal%ystop  = fld%grid%subdomain%internal%ystop
    else
      call gocean_stop('ERROR: ct_ne_init: implement periodic BCs!')
    end if

  end subroutine ct_ne_init

  !===================================================

  subroutine cf_field_init(fld)
    implicit none
    class(field_type), intent(inout) :: fld

    fld%defined_on = GO_F_POINTS

    select case(fld%grid%offset)

    case(GO_OFFSET_SW)
       call cf_sw_init(fld)

    case(GO_OFFSET_NE)
       call cf_ne_init(fld)

    case default
       call gocean_stop('cf_field_init: ERROR - unsupported grid offset!')

    end select

  end subroutine cf_field_init

  !===================================================

  subroutine cf_sw_init(fld)
    implicit none
    class(field_type), intent(inout) :: fld

    ! When updating a quantity on F points we write to:
    ! (using x to indicate a location that is written):
    !
    ! i=1         i=M
    !  o  b  b  b  b   j=N
    !  o  b  x  x  b
    !  o  b  x  x  b
    !  o  b  b  b  b
    !  o  o  o  o  o   j=1
    if(fld%grid%boundary_conditions(1) == GO_BC_PERIODIC)then
       fld%internal%xstart = fld%grid%subdomain%internal%xstart
       fld%internal%xstop  = fld%grid%subdomain%internal%xstop !internal%xstart + fld%grid%simulation_domain%nx - 1
    else
       fld%internal%xstart = fld%grid%subdomain%internal%xstart + 1
       fld%internal%xstop  = fld%grid%subdomain%internal%xstop
       ! I think these are correct but we stop because I've not properly
       ! gone through the coding.
       call gocean_stop('cf_sw_init: CHECK non-periodic BCs!')
    end if

    if(fld%grid%boundary_conditions(2) == GO_BC_PERIODIC)then
       fld%internal%ystart = fld%grid%subdomain%internal%ystart
       fld%internal%ystop  = fld%grid%subdomain%internal%ystop !fld%internal%ystart + fld%grid%simulation_domain%ny - 1
    else
       fld%internal%ystart = fld%grid%subdomain%internal%ystart + 1
       fld%internal%ystop  = fld%grid%subdomain%internal%ystop
       ! I think these are correct but we stop because I've not properly
       ! gone through the coding.
       call gocean_stop('cf_sw_init: CHECK non-periodic BCs!')
    end if

    ! When applying periodic (wrap-around) boundary conditions
    ! (PBCs) we must fill the regions marked with 'b' above.
    ! This looks like (using x to indicate a location that is
    ! written first and y a location that is written second):
    !
    !  i=2      i=M
    ! .-x  o  o  o   j=N
    ! | x  o  o  o
    ! | x  o  o  o
    ! | x  o  o  o
    ! ->y_ y  y  y   j=2
    !   |\______/

    ! In array notation this looks like:
    ! First col = last col
    ! field(2:2, 2:N) = field(M:M, 2:N)
    ! First row = last row
    ! field(2:M, 2:2) = field(2:M, N:N)

    call init_periodic_bc_halos(fld)

  end subroutine cf_sw_init

  !===================================================

  subroutine cf_ne_init(fld)
    implicit none
    class(field_type), intent(inout) :: fld

    ! When updating a quantity on F points we write to:
    ! (using x to indicate a location that is written
    !        b a boundary point - defined by ext. b.c.
    !        o a point that is external to the domain):
    !
    ! i=1       Nx
    !  o  o  o  o   Ny
    !  b  b  b  o
    !  b  x  b  o
    !  b  x  b  o
    !  b  b  b  o   j=1

    if(fld%grid%boundary_conditions(1) /= GO_BC_PERIODIC)then
      ! If we do not have periodic boundary conditions then we do
      ! not need to allow for boundary points here - they are
      ! already contained within the region.
      fld%internal%xstart = fld%grid%subdomain%internal%xstart
      fld%internal%xstop  = fld%grid%subdomain%internal%xstop
    else
      call gocean_stop('ERROR: cf_ne_init: implement periodic BCs!')
      stop
    end if

    if(fld%grid%boundary_conditions(2) /= GO_BC_PERIODIC)then
      fld%internal%ystart = fld%grid%subdomain%internal%ystart
      fld%internal%ystop  = fld%grid%subdomain%internal%ystop
    else
      call gocean_stop('ERROR: cf_ne_init: implement periodic BCs!')
    end if

  end subroutine cf_ne_init

  !===================================================

  SUBROUTINE copy_2dfield_array(field_in, field_out)
    IMPLICIT none
    REAL(go_wp), INTENT(in),  DIMENSION(:,:) :: field_in
    REAL(go_wp), INTENT(out), DIMENSION(:,:) :: field_out

    field_out(:,:) = field_in(:,:)

  end subroutine copy_2dfield_array

  !===================================================

  !> Copy from one patch in an array to another patch
  !! Will be superseded by interface using field object
  !! instead of array data.
  subroutine copy_2dfield_array_patch(field, src, dest)
    implicit none
    real(go_wp),       intent(inout), dimension(:,:) :: field
    type(region_type), intent(in)                    :: src, dest

    field(dest%xstart:dest%xstop,dest%ystart:dest%ystop) = &
     field(src%xstart:src%xstop ,src%ystart:src%ystop)

  end subroutine copy_2dfield_array_patch

  !===================================================

  SUBROUTINE copy_2dfield(field_in, field_out)
    IMPLICIT none
    type(r2d_field), intent(in)    :: field_in
    type(r2d_field), intent(inout) :: field_out
    integer :: it, ji, jj

    if (field_out%ntiles == 0) then
       field_out%data(:,:) = field_in%data(:,:)
    else
       !> \TODO #38 We don't actually support building dl_esm_inf
       !! with OpenMP enabled!
       !$OMP PARALLEL DO SCHEDULE(RUNTIME), DEFAULT(none), &
       !$OMP SHARED(field_out, field_in), PRIVATE(ji,jj)
       do it = 1, field_out%ntiles, 1
          do jj= field_out%tile(it)%whole%ystart, field_out%tile(it)%whole%ystop
             do ji = field_out%tile(it)%whole%xstart, field_out%tile(it)%whole%xstop
                field_out%data(ji,jj) = field_in%data(ji,jj)
             end do
          end do
       end do
       !$OMP END PARALLEL DO
    end if
  end subroutine copy_2dfield

  !===================================================

  !> Copy from one patch to another in a field
  subroutine copy_2dfield_patch(field, src, dest)
    implicit none
    type(r2d_field), intent(inout) :: field
    type(region_type),    intent(in)    :: src, dest

    field%data(dest%xstart:dest%xstop,dest%ystart:dest%ystop) = &
     field%data(src%xstart:src%xstop ,src%ystart:src%ystop)

  end subroutine copy_2dfield_patch

  !===================================================

  SUBROUTINE set_field(fld, val)
    implicit none
    class(field_type), INTENT(out) :: fld
    real(go_wp), INTENT(in) :: val

    select type(fld)
    type is (r2d_field)
       fld%data = val
    class default
    end select

  END SUBROUTINE set_field

  !===================================================

  !> Compute the checksum of the internal values of the supplied field
  !> object. If the data is in a remote device, the field%get_data()
  !> call will first bring it to the host.
  function fld_checksum(field) result(val)
    implicit none
    type(r2d_field), intent(in) :: field
    real(go_wp) :: val

    val = array_checksum(field%get_data(),           &
                         field%internal%xstart, field%internal%xstop, &
                         field%internal%ystart, field%internal%ystop)
    return

  end function fld_checksum

  !===================================================

  !> Provide access to exchange_generic for halo swaps for this field.
  !! If the data is in a remote device, it will be synchronised with the host
  !! before and after the exchange.
  !! The depth value is currently ignored and hardwired to 1: the device
  !! synchonization and the halo_exchange communicators need to be updated
  !! with the provided depth.
  !!
  !! TODO #58: This method could benefit from an asynchronous execution model.
  subroutine halo_exchange(self, depth)
    use parallel_comms_mod, only: Iplus, Iminus, Jplus, Jminus, NONE, &
         exchange_generic
    implicit none
    class(r2d_field), target, intent(inout) :: self
    integer, intent(in) :: depth
    ! Locals
    integer :: exch  !> Handle for exchange

    ! Make sure the data from each halo is in the host before the communication
    call self%read_halo_from_device(comm=JPlus, blocking=.false.)
    call self%read_halo_from_device(comm=Jminus, blocking=.false.)
    call self%read_halo_from_device(comm=IPlus, blocking=.false.)
    call self%read_halo_from_device(comm=IMinus, blocking=.true.)

    ! Execute the halo exchange
    call exchange_generic(b2=self%data, handle=exch, &
                          comm1=JPlus, comm2=Jminus, comm3=IPlus, comm4=IMinus)

    ! Synchronize the data back to the accelerator device if necessary
    call self%write_halo_to_device(comm=JPlus, blocking=.false.)
    call self%write_halo_to_device(comm=Jminus, blocking=.false.)
    call self%write_halo_to_device(comm=IPlus, blocking=.false.)
    call self%write_halo_to_device(comm=IMinus, blocking=.true.)

  end subroutine halo_exchange

  subroutine read_halo_from_device(self, comm, blocking)
    use parallel_comms_mod, only: isrcsend, jsrcsend, nxsend, nysend
    implicit none
    class(r2d_field), target, intent(inout) :: self
    integer, intent(in) :: comm
    logical, intent(in) :: blocking

    ! Get location and dimension to read data for the given communicator
    if (isrcsend(comm) > 0) then
        call self%read_from_device(isrcsend(comm), jsrcsend(comm), &
            nxsend(comm), nysend(comm), blocking)
    endif
  end subroutine read_halo_from_device

  subroutine write_halo_to_device(self, comm, blocking)
    use parallel_comms_mod, only: idesrecv, jdesrecv, nxrecv, nyrecv
    implicit none
    class(r2d_field), target, intent(inout) :: self
    integer, intent(in) :: comm
    logical, intent(in) :: blocking

    ! Get location and dimension to write data for the given communicator
    if (idesrecv(comm) > 0) then
        call self%write_to_device(idesrecv(comm), jdesrecv(comm), &
            nxrecv(comm), nyrecv(comm), blocking)
    endif
  end subroutine write_halo_to_device
  !===================================================

  !> Compute the checksum of ALL of the elements of supplied array. Performs a
  !! global sum if built with distributed-memory support.
  function array_checksum(field, &
                          xstart, xstop, &
                          ystart, ystop) result(val)
    use parallel_comms_mod, only: global_sum
    implicit none
    real(go_wp), dimension(:,:), intent(in) :: field
    integer, optional, intent(in) :: xstart, xstop, ystart, ystop
    real(go_wp) :: val

    if( present(xstart) )then
       val = SUM( ABS(field(xstart:xstop,ystart:ystop) ) )
    else
       val = SUM( ABS(field(:,:)) )
    end if

    ! Sum over all processors if running with distributed memory
    call global_sum(val)

  end function array_checksum

  !===================================================

  subroutine init_periodic_bc_halos(fld)
    implicit none
    class(field_type), intent(inout) :: fld
    ! Locals
    integer :: ihalo

    ! Check whether we have PBCs in x AND y dimensions
    fld%num_halos = 0
    if( fld%grid%boundary_conditions(1) == GO_BC_PERIODIC )then
       fld%num_halos = fld%num_halos + 2
    end if
    if( fld%grid%boundary_conditions(2) == GO_BC_PERIODIC )then
       fld%num_halos = fld%num_halos + 2
    end if

    allocate( fld%halo(fld%num_halos) )

    ihalo = 0
    if( fld%grid%boundary_conditions(1) == GO_BC_PERIODIC )then
       ! E-most column set to W-most internal column
       ihalo = ihalo + 1
       fld%halo(ihalo)%dest%xstart = fld%internal%xstop + 1
       fld%halo(ihalo)%dest%xstop  = fld%internal%xstop + 1
       fld%halo(ihalo)%dest%ystart = fld%internal%ystart
       fld%halo(ihalo)%dest%ystop  = fld%internal%ystop

       fld%halo(ihalo)%source%xstart = fld%internal%xstart
       fld%halo(ihalo)%source%xstop  = fld%internal%xstart
       fld%halo(ihalo)%source%ystart = fld%internal%ystart
       fld%halo(ihalo)%source%ystop  = fld%internal%ystop

       ! W-most column set to E-most internal column
       ihalo = ihalo + 1
       fld%halo(ihalo)%dest%xstart = fld%internal%xstart-1
       fld%halo(ihalo)%dest%xstop  = fld%internal%xstart-1
       fld%halo(ihalo)%dest%ystart = fld%internal%ystart
       fld%halo(ihalo)%dest%ystop  = fld%internal%ystop

       fld%halo(ihalo)%source%xstart = fld%internal%xstop
       fld%halo(ihalo)%source%xstop  = fld%internal%xstop
       fld%halo(ihalo)%source%ystart = fld%internal%ystart
       fld%halo(ihalo)%source%ystop  = fld%internal%ystop
    end if

    if( fld%grid%boundary_conditions(2) == GO_BC_PERIODIC )then
       ! N-most row set to S-most internal row
       ihalo = ihalo + 1
       fld%halo(ihalo)%dest%xstart = fld%internal%xstart - 1
       fld%halo(ihalo)%dest%xstop  = fld%internal%xstop  + 1
       fld%halo(ihalo)%dest%ystart = fld%internal%ystop+1
       fld%halo(ihalo)%dest%ystop  = fld%internal%ystop+1

       fld%halo(ihalo)%source%xstart = fld%internal%xstart - 1
       fld%halo(ihalo)%source%xstop  = fld%internal%xstop  + 1
       fld%halo(ihalo)%source%ystart = fld%internal%ystart
       fld%halo(ihalo)%source%ystop  = fld%internal%ystart

       ! S-most row set to N-most internal row
       ihalo = ihalo + 1
       fld%halo(ihalo)%dest%xstart = fld%internal%xstart - 1
       fld%halo(ihalo)%dest%xstop  = fld%internal%xstop  + 1
       fld%halo(ihalo)%dest%ystart = fld%internal%ystart - 1
       fld%halo(ihalo)%dest%ystop  = fld%internal%ystart - 1

       fld%halo(ihalo)%source%xstart = fld%internal%xstart - 1
       fld%halo(ihalo)%source%xstop  = fld%internal%xstop  + 1
       fld%halo(ihalo)%source%ystart = fld%internal%ystop
       fld%halo(ihalo)%source%ystop  = fld%internal%ystop
    end if

  end subroutine init_periodic_bc_halos

  !==============================================

  !> Routine to get the dimensions of the OpenMP tiling grid.
  !! Reads the GOCEAN_OMP_GRID environment variable which
  !! should have the format "NxM" where N is nx and M is ny.
  !! Returns false if the environment variable is not set
  !! or does not conform to this format.
  function get_grid_dims(nx, ny) result(success)
    implicit none
    integer, intent(inout) :: nx, ny
    logical :: success
    character(len=20) :: lstr
    integer :: idx, ierr

    success = .FALSE.

    call get_environment_variable(NAME='GOCEAN_OMP_GRID', VALUE=lstr, &
                                  STATUS=ierr)

    if(ierr /= 0)return

    ! We expect the string to have the format 'AxB' where A and B are
    ! integers.
    idx = index(lstr, 'x')
    if(idx == 0)then
       write (*,"(/'shallow_omp_mod::get_grid_dims: failed to parse ' &
                 &  'GOCEAN_OMP_GRID string: ',(A))") TRIM(lstr)
       write (*,"('   -  will use defaults for dimensions of tiling grid')")
       return
    endif

    read(lstr(1:idx-1),*,iostat=ierr) nx
    if(ierr /= 0)return

    read(lstr(idx+1:),*,iostat=ierr) ny
    if(ierr == 0)success = .TRUE.

  end function get_grid_dims

end module field_mod


! { dg-do compile }
! { dg-options "-Wunused-parameter" }
! PR66605
MODULE test
  IMPLICIT NONE
  INTEGER, PARAMETER :: wp = KIND(1.0D0)
CONTAINS
SUBROUTINE sub (neq, time, y, dydt)
  IMPLICIT NONE
    INTEGER :: neq
    REAL(WP) :: time, y(neq), dydt(neq)

    dydt(1) = 1.0 / y(1)
END SUBROUTINE sub
END MODULE


text:
subroutine foo(a, n, m, b)
  implicit none

  real, intent(in) :: a(n, m)
  integer, intent(in) :: n, m
  real, intent(out) :: b(size(a, 1))

  integer :: i

  do i = 1, size(b)
    b(i) = sum(a(i,:))
  enddo
end subroutine

subroutine trans(x,y)
  implicit none
  real, intent(in), dimension(:,:) :: x
  real, intent(out), dimension( size(x,2), size(x,1) ) :: y
  integer :: N, M, i, j
  N = size(x,1)
  M = size(x,2)
  DO i=1,N
     do j=1,M
        y(j,i) = x(i,j)
     END DO
  END DO
end subroutine trans

subroutine flatten(x,y)
  implicit none
  real, intent(in), dimension(:,:) :: x
  real, intent(out), dimension( size(x) ) :: y
  integer :: N, M, i, j, k
  N = size(x,1)
  M = size(x,2)
  k = 1
  DO i=1,N
     do j=1,M
        y(k) = x(i,j)
        k = k + 1
     END DO
  END DO
end subroutine flatten

max_stars_repo_name: rocketbot-cl/recognition

max_issues_repo_name: rocketbot-cl/recognition

max_forks_repo_name: rocketbot-cl/recognition


text:
module stdlib_stats_distribution_PRNG
    use stdlib_kinds, only: int8, int16, int32, int64
    use stdlib_error, only: error_stop
    implicit none
    private
    integer, parameter :: MAX_INT_BIT_SIZE = bit_size(1_int64)
    integer(int64) :: st(4) ! internal states for xoshiro256ss function
    integer(int64) :: si = 614872703977525537_int64 ! default seed value
    logical :: seed_initialized = .false.

    public :: random_seed
    public :: dist_rand

    interface dist_rand
    !! Version experimental
    !!
    !! Generation of random integers with different kinds
    !! ([Specification](../page/specs/stdlib_stats_distribution_PRNG.html#
    !! description))
        module procedure dist_rand_iint8
        module procedure dist_rand_iint16
        module procedure dist_rand_iint32
        module procedure dist_rand_iint64
    end interface dist_rand

    interface random_seed
    !! Version experimental
    !!
    !! Set seed value for random number generator
    !! ([Specification](../page/specs/stdlib_stats_distribution_PRNG.html#
    !! description))
    !!
        module procedure random_distribution_seed_iint8
        module procedure random_distribution_seed_iint16
        module procedure random_distribution_seed_iint32
        module procedure random_distribution_seed_iint64
    end interface random_seed

    contains

    function dist_rand_iint8(n) result(res)
    !! Random integer generation for various kinds
    !! result = [-2^k, 2^k - 1], k = 7, 15, 31, 63, depending on input kind
    !! Result will be operated by bitwise operators to generate desired integer
    !! and real pseudorandom numbers
    !!
        integer(int8), intent(in) :: n
        integer(int8) :: res
        integer :: k

        k = MAX_INT_BIT_SIZE - bit_size(n)
        if(k < 0) call error_stop("Error(dist_rand): Integer bit size is"       &
                  //" greater than 64bit")
        res = shiftr(xoshiro256ss( ), k)
    end function dist_rand_iint8

    function dist_rand_iint16(n) result(res)
    !! Random integer generation for various kinds
    !! result = [-2^k, 2^k - 1], k = 7, 15, 31, 63, depending on input kind
    !! Result will be operated by bitwise operators to generate desired integer
    !! and real pseudorandom numbers
    !!
        integer(int16), intent(in) :: n
        integer(int16) :: res
        integer :: k

        k = MAX_INT_BIT_SIZE - bit_size(n)
        if(k < 0) call error_stop("Error(dist_rand): Integer bit size is"       &
                  //" greater than 64bit")
        res = shiftr(xoshiro256ss( ), k)
    end function dist_rand_iint16

    function dist_rand_iint32(n) result(res)
    !! Random integer generation for various kinds
    !! result = [-2^k, 2^k - 1], k = 7, 15, 31, 63, depending on input kind
    !! Result will be operated by bitwise operators to generate desired integer
    !! and real pseudorandom numbers
    !!
        integer(int32), intent(in) :: n
        integer(int32) :: res
        integer :: k

        k = MAX_INT_BIT_SIZE - bit_size(n)
        if(k < 0) call error_stop("Error(dist_rand): Integer bit size is"       &
                  //" greater than 64bit")
        res = shiftr(xoshiro256ss( ), k)
    end function dist_rand_iint32

    function dist_rand_iint64(n) result(res)
    !! Random integer generation for various kinds
    !! result = [-2^k, 2^k - 1], k = 7, 15, 31, 63, depending on input kind
    !! Result will be operated by bitwise operators to generate desired integer
    !! and real pseudorandom numbers
    !!
        integer(int64), intent(in) :: n
        integer(int64) :: res
        integer :: k

        k = MAX_INT_BIT_SIZE - bit_size(n)
        if(k < 0) call error_stop("Error(dist_rand): Integer bit size is"       &
                  //" greater than 64bit")
        res = shiftr(xoshiro256ss( ), k)
    end function dist_rand_iint64

    function xoshiro256ss( ) result (res)
    ! Generate random 64-bit integers
    ! Written in 2018 by David Blackman and Sebastiano Vigna (vigna@acm.org)
    ! http://prng.di.unimi.it/xoshiro256starstar.c
    !
    ! This is xoshiro256** 1.0, one of our all-purpose, rock-solid
    ! generators. It has excellent (sub-ns) speed, a state (256 bits) that is
    ! large enough for any parallel application, and it passes all tests we
    ! are aware of.
    !
    ! The state must be seeded so that it is not everywhere zero. If you have
    ! a 64-bit seed, we suggest to seed a splitmix64 generator and use its
    ! output to fill st.
    !
    ! Fortran 90 version translated from C by Jim-215-Fisher
    !
        integer(int64) :: res, t

        if(.not. seed_initialized) call random_distribution_seed_iint64(si,t)
        res = rol64(st(2) * 5, 7) * 9
        t = shiftl(st(2), 17)
        st(3) = ieor(st(3), st(1))
        st(4) = ieor(st(4), st(2))
        st(2) = ieor(st(2), st(3))
        st(1) = ieor(st(1), st(4))
        st(3) = ieor(st(3), t)
        st(4) = rol64(st(4), 45)
    end function xoshiro256ss

    pure function rol64(x, k) result(res)
        integer(int64), intent(in) :: x
        integer, intent(in) :: k
        integer(int64) :: t1, t2, res

        t1 = shiftr(x, (64 - k))
        t2 = shiftl(x, k)
        res = ior(t1, t2)
    end function rol64

    function splitmix64(s) result(res)
    ! Written in 2015 by Sebastiano Vigna (vigna@acm.org)
    ! This is a fixed-increment version of Java 8's SplittableRandom
    ! generator.
    ! See http://dx.doi.org/10.1145/2714064.2660195 and
    ! http://docs.oracle.com/javase/8/docs/api/java/util/SplittableRandom.html
    !
    ! It is a very fast generator passing BigCrush, and it can be useful if
    ! for some reason you absolutely want 64 bits of state.
    !
    ! Fortran 90 translated from C by Jim-215-Fisher
    !
        integer(int64) :: res
        integer(int64), intent(in), optional :: s
        integer(int64) :: int01 = -7046029254386353131_int64,                  &
                          int02 = -4658895280553007687_int64,                  &
                          int03 = -7723592293110705685_int64
    ! Values are converted from C unsigned integer of 0x9e3779b97f4a7c15,
    ! 0xbf58476d1ce4e5b9, 0x94d049bb133111eb

        if(present(s)) si = s
        res = si
        si = res + int01
        res = ieor(res, shiftr(res, 30)) * int02
        res = ieor(res, shiftr(res, 27)) * int03
        res = ieor(res, shiftr(res, 31))
    end function splitmix64

    subroutine random_distribution_seed_iint8(put, get)
    !! Set seed value for random number generator
    !!
        integer(int8), intent(in) :: put
        integer(int8), intent(out) :: get
        integer(int64) :: tmp
        integer :: i

        tmp = splitmix64(int(put, kind = int64))
        do i = 1, 10
            tmp = splitmix64( )
        end do
        do i = 1, 4
            tmp = splitmix64( )
            st(i) = tmp
        end do
        get = int(tmp, kind = int8)
        seed_initialized = .true.
    end subroutine random_distribution_seed_iint8

    subroutine random_distribution_seed_iint16(put, get)
    !! Set seed value for random number generator
    !!
        integer(int16), intent(in) :: put
        integer(int16), intent(out) :: get
        integer(int64) :: tmp
        integer :: i

        tmp = splitmix64(int(put, kind = int64))
        do i = 1, 10
            tmp = splitmix64( )
        end do
        do i = 1, 4
            tmp = splitmix64( )
            st(i) = tmp
        end do
        get = int(tmp, kind = int16)
        seed_initialized = .true.
    end subroutine random_distribution_seed_iint16

    subroutine random_distribution_seed_iint32(put, get)
    !! Set seed value for random number generator
    !!
        integer(int32), intent(in) :: put
        integer(int32), intent(out) :: get
        integer(int64) :: tmp
        integer :: i

        tmp = splitmix64(int(put, kind = int64))
        do i = 1, 10
            tmp = splitmix64( )
        end do
        do i = 1, 4
            tmp = splitmix64( )
            st(i) = tmp
        end do
        get = int(tmp, kind = int32)
        seed_initialized = .true.
    end subroutine random_distribution_seed_iint32

    subroutine random_distribution_seed_iint64(put, get)
    !! Set seed value for random number generator
    !!
        integer(int64), intent(in) :: put
        integer(int64), intent(out) :: get
        integer(int64) :: tmp
        integer :: i

        tmp = splitmix64(int(put, kind = int64))
        do i = 1, 10
            tmp = splitmix64( )
        end do
        do i = 1, 4
            tmp = splitmix64( )
            st(i) = tmp
        end do
        get = int(tmp, kind = int64)
        seed_initialized = .true.
    end subroutine random_distribution_seed_iint64

end module stdlib_stats_distribution_PRNG

max_stars_repo_name: LKedward/stdlib-fpm

max_issues_repo_name: LKedward/stdlib-fpm

max_forks_repo_name: LKedward/stdlib-fpm


text:
program grid_3d

  integer :: i, j, k

  integer, parameter :: nx = 2
  integer, parameter :: ny = 2
  integer, parameter :: nz = 24
  integer, parameter :: nzonesx = 64
  integer, parameter :: nzonesy = 64
  integer, parameter :: nzonesz = 768

  integer :: ix, iy, iz
  integer :: nlevs, ngrids

 99 format(i1)
100 format('         ((',i4,',',i4,',',i4,') (',i4,',',i4,',',i4,') ('i4,',',i4,',',i4,'))')
101 format('   ((',i4,',',i4,',',i4,') (',i4,',',i4,',',i4,') ('i4,',',i4,',',i4,')) ', i8)

  ix = 0
  nlevs = 1
  ngrids = nx*ny*nz

  write (*, 99) nlevs
  write (*,101) 0,0,0,nzonesx-1,nzonesy-1,nzonesz-1,0,0,0,ngrids
  do i = 1, nx

     iy = 0
     do j = 1, ny

        iz = 0
        do k = 1, nz

           write (*,100) ix,iy,iz, ix+nzonesx/nx-1,iy+nzonesy/ny-1,iz+nzonesz/nz-1,0,0,0

           iz = iz + nzonesz/nz
        enddo

        iy = iy + nzonesy/ny
     enddo

     ix = ix + nzonesx/nx
  enddo

end program grid_3d

max_stars_repo_name: sailoridy/MAESTRO

max_issues_repo_name: sailoridy/MAESTRO

max_forks_repo_name: sailoridy/MAESTRO


text:
MODULE  M1
  IMPLICIT NONE
  TYPE mmm
    COMPLEX(KIND=8), DIMENSION(:,:), POINTER :: data
  END TYPE mmm

CONTAINS

  SUBROUTINE S(ma,mb,mc)
      TYPE(mmm), POINTER :: ma,mb,mc
      COMPLEX(KIND=8), DIMENSION(:, :), &
        POINTER                                :: a, b, c
      INTEGER :: i,j
      a=>ma%data
      b=>mb%data
      c=>mc%data
      DO i=1,size(a,1)
      DO j=1,size(a,2)
         c(i,j)=a(i,j)*b(i,j)
      ENDDO
      ENDDO
  END SUBROUTINE

END MODULE M1

max_stars_repo_name: maurizioabba/rose

max_issues_repo_name: sujankh/rose-matlab

max_forks_repo_name: sujankh/rose-matlab


text:
C
C
*+ stack_init

       subroutine stack_init(nmaps,smap,nbuffers,sbuffer,status)
C      ----------------------------------------------------------
C
C Initialise the stack
C
C Given:
C     number of and size of map stack entries
        integer      nmaps
        integer      smap
C     number of and size of buffer entries
        integer      nbuffers
        integer      sbuffer
C
C Returned:
C     Status
        integer      status
C
C Initialise the map stack. Two functions are performed by the routine:
C    (a) define the total number of map slots and buffers and sizes
C    (b) initialise all the pointers and map_status block
C
C If the number of maps specified in the call to this routine is less
C than zero then only function (b) is carried out -- in this way the
C stack can be re-initialised to the same size as in the initial
C definition.
*-
       include 'mapcat_stack.inc'

       integer   i, j, is, ip, irange

       if (status.ne.0) return

       if (nmaps.gt.0) then
         number_map_entries   = nmaps
         number_stack_entries = number_redtape_entries +
     *                          number_map_entries +
     *                          number_work_entries
         size_map_entry       = smap
         size_work_entry      = smap
         number_map_buffers   = nbuffers
         size_map_buffer      = sbuffer
       end if
       do i=1,max_active_maps
         active_list(1,i) = 0
       end do
       do i=1,number_map_buffers
         buffer_status(i) = 0
       end do
       is = 0
       ip = 1
       do i = 1,number_redtape_entries
         is = is + 1
         do j = 1,length_stack_status - 3
           stack_status(j,is) = 0
         end do
         stack_status(ip_type,is)     = type_redtape
         stack_status(ip_pointer,is)  = ip
         stack_status(ip_size,is)     = size_redtape_entry
         ip = ip + stack_status(ip_size,is)
       end do
       ip = 1
       do i = 1,number_map_entries
         is = is + 1
         do j = 1,length_stack_status - 3
           stack_status(j,is) = 0
         end do
         stack_status(ip_type,is)     = type_map
         stack_status(ip_pointer,is)  = ip
         stack_status(ip_size,is)     = size_map_entry
         ip = ip + stack_status(ip_size,is)
       end do
       ip = 1
       if (number_work_entries.ge.1) then
         irange = number_work_entries
         do i = 1,irange
           is = is + 1
           do j = 1,length_stack_status - 3
             stack_status(j,is) = 0
           end do
           stack_status(ip_type,is)     = type_work
           stack_status(ip_pointer,is)  = ip
           stack_status(ip_size,is)  = size_work_entry
           ip = ip + stack_status(ip_size,is)
         end do
       end if
       call mapcat_err(status,'stack_init',
     *                 'Stack Initialisation Failed')
       end

max_stars_repo_name: CavendishAstrophysics/anmap

max_issues_repo_name: CavendishAstrophysics/anmap

max_forks_repo_name: CavendishAstrophysics/anmap