text:
!
!     file hwsssp.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 hwssp(ts, tf, m, mbdcnd, bdts, bdtf, ps, pf, n, nbdcnd, bdps,
!                      bdpf, elmbda, f, idimf, pertrb, ierror)
!
! DIMENSION OF           bdts(n+1),    bdtf(n+1), bdps(m+1), bdpf(m+1),
! ARGUMENTS              f(idimf, n+1)
!
! LATEST REVISION        May 2016
!
! PURPOSE                Solves a finite difference approximation to
!                        the helmholtz equation in spherical
!                        coordinates and on the surface of the unit
!                        sphere (radius of 1).  the equation is
!
!                          (1/sin(theta))(d/dtheta)(sin(theta)
!                          (du/dtheta)) + (1/sin(theta)**2)(d/dphi)
!                          (du/dphi)  + lambda*u = f(theta, phi)
!
!                        where theta is colatitude and phi is
!                        longitude.
!
! USAGE                  call hwsssp(ts, tf, m, mbdcnd, bdts, bdtf, ps, pf,
!                                     n, nbdcnd, bdps, bdpf, elmbda, f,
!                                     idimf, pertrb, ierror, w)
!
! ARGUMENTS
! ON INPUT               ts, tf
!
!                          the range of theta (colatitude), i.e.,
!                          ts <= theta <= tf. ts must be less
!                          than tf.  ts and tf are in radians.
!                          a ts of zero corresponds to the north
!                          pole and a tf of pi corresponds to
!                          the south pole.
!
!                          * * * important * * *
!
!                          if tf is equal to pi then it must be
!                          computed using the statement
!                          tf = pi_mach(dum). this insures that tf
!                          in the user's program is equal to pi in
!                          this program which permits several tests
!                          of the input parameters that otherwise
!                          would not be possible.
!
!
!                        m
!                          the number of panels into which the
!                          interval (ts, tf) is subdivided.
!                          hence, there will be m+1 grid points in the
!                          theta-direction given by
!                          theta(i) = (i-1)dtheta+ts for
!                          i = 1, 2, ..., m+1, where
!                          dtheta = (tf-ts)/m is the panel width.
!                          m must be greater than 5
!
!                        mbdcnd
!                          indicates the type of boundary condition
!                          at theta = ts and theta = tf.
!
!                          = 1  if the solution is specified at
!                               theta = ts and theta = tf.
!                          = 2  if the solution is specified at
!                               theta = ts and the derivative of
!                               the solution with respect to theta is
!                               specified at theta = tf
!                               (see note 2 below).
!                          = 3  if the derivative of the solution
!                               with respect to theta is specified
!                               specified at theta = ts and
!                               theta = tf (see notes 1, 2 below).
!                          = 4  if the derivative of the solution
!                               with respect to theta is specified
!                               at theta = ts (see note 1 below)
!                               and the solution is specified at
!                               theta = tf.
!                          = 5  if the solution is unspecified at
!                               theta = ts = 0 and the solution
!                               is specified at theta = tf.
!                          = 6  if the solution is unspecified at
!                               theta = ts = 0 and the derivative
!                               of the solution with respect to theta
!                               is specified at theta = tf
!                               (see note 2 below).
!                          = 7  if the solution is specified at
!                               theta = ts and the solution is
!                               is unspecified at theta = tf = pi.
!                          = 8  if the derivative of the solution
!                               with respect to theta is specified
!                               at theta = ts (see note 1 below) and
!                               the solution is unspecified at
!                               theta = tf = pi.
!                          = 9  if the solution is unspecified at
!                               theta = ts = 0 and theta = tf = pi.
!
!                          notes:
!                          if ts = 0, do not use mbdcnd = 3, 4, or 8,
!                          but instead use mbdcnd = 5, 6, or 9  .
!
!                          if tf = pi, do not use mbdcnd = 2, 3, or 6,
!                          but instead use mbdcnd = 7, 8, or 9  .
!
!                        bdts
!                          a one-dimensional array of length n+1 that
!                          specifies the values of the derivative of
!                          the solution with respect to theta at
!                          theta = ts.  when mbdcnd = 3, 4, or 8,
!
!                          bdts(j) = (d/dtheta)u(ts, phi(j)),
!                          j = 1, 2, ..., n+1  .
!
!                          when mbdcnd has any other value, bdts is
!                          a dummy variable.
!
!                        bdtf
!                          a one-dimensional array of length n+1
!                          that specifies the values of the derivative
!                          of the solution with respect to theta at
!                          theta = tf.  when mbdcnd = 2, 3, or 6,
!
!                          bdtf(j) = (d/dtheta)u(tf, phi(j)),
!                          j = 1, 2, ..., n+1  .
!
!                          when mbdcnd has any other value, bdtf is
!                          a dummy variable.
!
!                        ps, pf
!                          the range of phi (longitude), i.e.,
!                          ps <= phi <= pf.  ps must be less
!                          than pf.  ps and pf are in radians.
!                          if ps = 0 and pf = 2*pi, periodic
!                          boundary conditions are usually prescribed.
!
!                          * * * IMPORTANT * * *
!
!                          if pf is equal to 2*pi then it must be
!                          computed using the statement
!                          pf = TWO *pi_mach(dum). this insures that
!                          pf in the users program is equal to
!                          2*pi in this program which permits tests
!                          of the input parameters that otherwise
!                          would not be possible.
!
!                        n
!                          the number of panels into which the
!                          interval (ps, pf) is subdivided.
!                          hence, there will be n+1 grid points
!                          in the phi-direction given by
!                          phi(j) = (j-1)dphi+ps  for
!                          j = 1, 2, ..., n+1, where
!                          dphi = (pf-ps)/n is the panel width.
!                          n must be greater than 4
!
!                        nbdcnd
!                          indicates the type of boundary condition
!                          at phi = ps and phi = pf.
!
!                          = 0  if the solution is periodic in phi,
!                               i.u., u(i, j) = u(i, n+j).
!                          = 1  if the solution is specified at
!                               phi = ps and phi = pf
!                               (see note below).
!                          = 2  if the solution is specified at
!                               phi = ps (see note below)
!                               and the derivative of the solution
!                               with respect to phi is specified
!                               at phi = pf.
!                          = 3  if the derivative of the solution
!                               with respect to phi is specified
!                               at phi = ps and phi = pf.
!                          = 4  if the derivative of the solution
!                               with respect to phi is specified
!                               at ps and the solution is specified
!                               at phi = pf
!
!                          note:
!                          nbdcnd = 1, 2, or 4 cannot be used with
!                          mbdcnd = 5, 6, 7, 8, or 9.  the former indicates
!                          that the solution is specified at a pole, the
!                          latter indicates that the solution is not
!                          specified.  use instead  mbdcnd = 1 or 2.
!
!                        bdps
!                          a one-dimensional array of length m+1 that
!                          specifies the values of the derivative
!                          of the solution with respect to phi at
!                          phi = ps.  when nbdcnd = 3 or 4,
!
!                            bdps(i) = (d/dphi)u(theta(i), ps),
!                            i = 1, 2, ..., m+1  .
!
!                          when nbdcnd has any other value, bdps is
!                          a dummy variable.
!
!                        bdpf
!                          a one-dimensional array of length m+1 that
!                          specifies the values of the derivative
!                          of the solution with respect to phi at
!                          phi = pf.  when nbdcnd = 2 or 3,
!
!                            bdpf(i) = (d/dphi)u(theta(i), pf),
!                            i = 1, 2, ..., m+1  .
!
!                          when nbdcnd has any other value, bdpf is
!                          a dummy variable.
!
!                        elmbda
!                          the constant lambda in the helmholtz
!                          equation.  if lambda > 0, a solution
!                          may not exist.  however, hwsssp will
!                          attempt to find a solution.
!
!                        f
!                          a two-dimensional array that specifies the
!                          value of the right side of the helmholtz
!                          equation and boundary values (if any).
!                          f must be dimensioned at least (m+1)*(n+1).
!
!                          on the interior, f is defined as follows:
!                          for i = 2, 3, ..., m and j = 2, 3, ..., n
!                          f(i, j) = f(theta(i), phi(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(ts, phi(j))      u(tf, phi(j))
!                            2      u(ts, phi(j))      f(tf, phi(j))
!                            3      f(ts, phi(j))      f(tf, phi(j))
!                            4      f(ts, phi(j))      u(tf, phi(j))
!                            5      f(0, ps)           u(tf, phi(j))
!                            6      f(0, ps)           f(tf, phi(j))
!                            7      u(ts, phi(j))      f(pi, ps)
!                            8      f(ts, phi(j))      f(pi, ps)
!                            9      f(0, ps)           f(pi, ps)
!
!                          nbdcnd   f(i, 1)            f(i, n+1)
!                          ------   --------------    --------------
!
!                            0      f(theta(i), ps)    f(theta(i), ps)
!                            1      u(theta(i), ps)    u(theta(i), pf)
!                            2      u(theta(i), ps)    f(theta(i), pf)
!                            3      f(theta(i), ps)    f(theta(i), pf)
!                            4      f(theta(i), ps)    u(theta(i), pf)
!
!                          note:
!                          if the table calls for both the solution u
!                          and the right side f at a corner then the
!                          solution must be specified.
!
!                        idimf
!                          the row (or first) dimension of the array
!                          f as it appears in the program calling
!                          hwsssp.  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
!                          (theta(i), phi(j)),  i = 1, 2, ..., m+1  and
!                          j = 1, 2, ..., n+1  .
!
!                        pertrb
!                          if one specifies a combination of periodic,
!                          derivative or unspecified boundary
!                          conditions 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.  hwsssp then computes
!                          this solution, which is a least squares
!                          solution to the original approximation.
!                          this solution is not unique and is
!                          unnormalized. the value of pertrb should
!                          be small compared to the right side f.
!                          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 8,
!                          a solution is not attempted.
!
!                          = 0  no error
!                          = 1  ts<0 or tf>pi
!                          = 2  ts>=tf
!                          = 3  mbdcnd<1 or mbdcnd>9
!                          = 4  ps<0 or ps>pi+pi
!                          = 5  ps>=pf
!                          = 6  n<5
!                          = 7  m<5
!                          = 8  nbdcnd<0 or nbdcnd>4
!                          = 9  elmbda>0
!                          = 10 idimf<m+1
!                          = 11 nbdcnd equals 1, 2 or 4 and mbdcnd>=5
!                          = 12 ts==0 and mbdcnd equals 3, 4 or 8
!                          = 13 tf==pi and mbdcnd equals 2, 3 or 6
!                          = 14 mbdcnd equals 5, 6 or 9 and ts/=0
!                          = 15 mbdcnd>=7 and tf/=pi
!                          = 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)
!
!
! SPECIAL CONDITIONS     None
!
! I/O                    None
!
! PRECISION              64-bit double precision
!
! REQUIRED files         type_FishpackWorkspace.f90, genbun.f90, type_CyclicReductionUtility.f9090
!
! 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.
!
! PORTABILITY            Fortran 2008
!
! 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             P. N. Swarztrauber, "The direct solution of
!                        the discrete poisson equation on the surface of
!                        a sphere", S.I.A.M. J. Numer. Anal., 15(1974),
!                        pp 212-215.
!
!                        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_spherical_solver

! Parameters confined to the submodule

integer(ip), parameter :: IIWK = 7 ! Size for workspace_indices

contains

    module subroutine hwsssp(ts, tf, m, mbdcnd, bdts, bdtf, ps, pf, n, nbdcnd, &
        bdps, bdpf, 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)     :: ts
        real(wp),    intent(in)     :: tf
        real(wp),    intent(in)     :: ps
        real(wp),    intent(in)     :: pf
        real(wp),    intent(in)     :: elmbda
        real(wp),    intent(out)    :: pertrb
        real(wp),    intent(in)     :: bdts(:)
        real(wp),    intent(in)     :: bdtf(:)
        real(wp),    intent(in)     :: bdps(:)
        real(wp),    intent(in)     :: bdpf(:)
        real(wp),    intent(inout)  :: f(:,:)

        ! Local variables
        type(FishpackWorkspace) :: workspace

        ! Check input arguments
        call hwsssp_check_input_arguments(ts, tf, m, mbdcnd, ps, pf, n, &
            nbdcnd, elmbda, idimf, ierror)

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

        ! Allocate memory
        call hwsssp_initialize_workspace(n, m, workspace)

        ! Solve system
        associate( &
            indx => workspace%workspace_indices, &
            rew => workspace%real_workspace &
            )
            associate( &
                am => rew(indx(1):), &
                bm => rew(indx(2):), &
                cm => rew(indx(3):), &
                sn => rew(indx(4):), &
                ss => rew(indx(5):), &
                sint => rew(indx(6):), &
                d => rew(indx(7):) &
                )
                call hwsssp_lower_routine(ts, tf, m, mbdcnd, bdts, bdtf, ps, pf, n, nbdcnd, &
                    bdps, bdpf, elmbda, f, idimf, pertrb, am, bm, cm, sn, ss, &
                    sint, d, ierror)
            end associate
        end associate

        ! Release memory
        call workspace%destroy()

    end subroutine hwsssp

    subroutine hwsssp_check_input_arguments(ts, tf, m, mbdcnd, ps, pf, n, &
        nbdcnd, elmbda, 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)     :: elmbda
        real(wp),    intent(in)     :: ts
        real(wp),    intent(in)     :: tf
        real(wp),    intent(in)     :: ps
        real(wp),    intent(in)     :: pf

        if (ts < ZERO .or. tf > PI) then
            ierror = 1
        else if (ts >= tf) then
            ierror = 2
        else if (mbdcnd < 1 .or. mbdcnd > 9) then
            ierror = 3
        else if (ps < ZERO .or. pf > TWO_PI) then
            ierror = 4
        else if (ps >= pf) then
            ierror = 5
        else if (n < 5) then
            ierror = 6
        else if (m < 5) then
            ierror = 7
        else if (nbdcnd < 0 .or. nbdcnd > 4) then
            ierror = 8
        else if (elmbda > ZERO) then
            ierror = 9
        else if (idimf < m + 1) then
            ierror = 10
        else if ((nbdcnd == 1 .or. nbdcnd == 2 .or. nbdcnd == 4) .and. mbdcnd >= 5) then
            ierror = 11
        else if (ts == ZERO .and. (mbdcnd==3.or. mbdcnd==4 .or. mbdcnd == 8)) then
            ierror = 12
        else if (tf == PI .and. (mbdcnd==2.or. mbdcnd == 3 .or. mbdcnd == 6)) then
            ierror = 13
        else if ((mbdcnd == 5.or.mbdcnd == 6.or.mbdcnd == 9) .and. ts /= ZERO) then
            ierror = 14
        else if (mbdcnd >= 7 .and. tf /= PI) then
            ierror = 15
        else
            ierror = 0
        end if

    end subroutine hwsssp_check_input_arguments

    subroutine hwsssp_initialize_workspace(n, m, workspace)

        ! Dummy arguments
        integer(ip), intent(in)  :: n
        integer(ip), intent(in)  :: m
        class(FishpackWorkspace), intent(out) :: workspace

        ! Local variables
        integer(ip) :: irwk, icwk

        ! Compute workspace lengths for hwsssp
        irwk = 4 * (n + 1)  &
            + (16 + int(log(real(n+1,kind=wp))/log(TWO), kind=ip)) &
            * (m + 1)
        icwk = 0

        ! Allocate memory
        call workspace%create(irwk, icwk, IIWK)

        ! Set workspace indices
        workspace%workspace_indices = hwsssp_get_workspace_indices(m)

    end subroutine hwsssp_initialize_workspace

    pure function hwsssp_get_workspace_indices(m) &
        result (return_value)

        ! Dummy arguments
        integer(ip), intent(in) :: m
        integer(ip)             :: return_value(IIWK)

        associate( indx => return_value)
            indx(1) = 1
            indx(2) = m+2
            indx(3) = 2*m+3
            indx(4) = 3*m+4
            indx(5) = 4*m+5
            indx(6) = 5*m+6
            indx(7) = 6*m+7
        end associate

    end function hwsssp_get_workspace_indices

    subroutine hwsssp_lower_routine(ts, tf, m, mbdcnd, bdts, bdtf, ps, pf, n, nbdcnd, &
        bdps, bdpf, elmbda, f, idimf, pertrb, am, bm, cm, sn, ss, &
        sint, d, error_flag)

        ! 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
        real(wp),    intent(in)     :: ts
        real(wp),    intent(in)     :: tf
        real(wp),    intent(in)     :: ps
        real(wp),    intent(in)     :: pf
        real(wp),    intent(in)     :: elmbda
        real(wp),    intent(out)    :: pertrb
        real(wp),    intent(in)     :: bdts(:)
        real(wp),    intent(in)     :: bdtf(:)
        real(wp),    intent(in)     :: bdps(:)
        real(wp),    intent(in)     :: bdpf(:)
        real(wp),    intent(inout)  :: f(:,:)
        real(wp),    intent(out)    :: am(:)
        real(wp),    intent(out)    :: bm(:)
        real(wp),    intent(out)    :: cm(:)
        real(wp),    intent(out)    :: sn(:)
        real(wp),    intent(out)    :: ss(:)
        real(wp),    intent(out)    :: sint(:)
        real(wp),    intent(out)    :: d(:)
        integer(ip), intent(out)    :: error_flag

        ! Local variables
        integer(ip) :: mp1, np1, i, inp, isp, mbr, its, itf, itsp, itfm, munk
        integer(ip) :: iid, ii, nbr, jps, jpf, jpsp, jpfm, nunk, ising
        real(wp)    :: fn, fm, dth, half_dth, two_dth, dphi, two_dphi
        real(wp)    :: dphi2, edp2, dth2, cp, wpp, fim1, theta, t1, at, ct, wts, wtf
        real(wp)    :: wps, wpf, fjj, cf, summation, sum1, hne, yhld, sum2, dfn, dnn, dsn
        real(wp)    :: cnp, hld, dfs, dss, dns, csp, rtn, rts, den
        type(CenteredCyclicReductionUtility) :: util

        mp1 = m + 1
        np1 = n + 1
        fn = n
        fm = m
        dth = (tf - ts)/fm
        half_dth = dth/2
        two_dth = TWO * dth
        dphi = (pf - ps)/fn
        two_dphi = TWO * dphi
        dphi2 = dphi**2
        edp2 = elmbda*dphi2
        dth2 = dth**2
        cp = 4.0_wp/(fn*dth2)
        wpp = fn*sin(half_dth)/4

        ! Initialize constants
        dfn = ZERO
        dnn = ZERO
        dsn = ZERO

        do i = 1, mp1
            fim1 = i - 1
            theta = fim1*dth + ts
            sint(i) = sin(theta)

            if (sint(i) == ZERO) cycle

            t1 = ONE/(dth2*sint(i))
            am(i) = t1*sin(theta - half_dth)
            cm(i) = t1*sin(theta + half_dth)
            bm(i) = (-am(i)) - cm(i) + elmbda
        end do
        inp = 0
        isp = 0
        !
        ! boundary condition at theta=ts
        !
        mbr = mbdcnd + 1
        select case (mbr)
            case (1)
                its = 1
            case (2:3, 8)
                at = am(2)
                its = 2
            case (4:5, 9)
                at = am(1)
                its = 1
                cm(1) = am(1) + cm(1)
            case (6:7, 10)
                at = am(2)
                inp = 1
                its = 2
        end select

        !
        ! boundary condition theta=tf
        !
        select case (mbr)
            case (1)
                itf = m
            case (2, 5:6)
                ct = cm(m)
                itf = m
            case (3:4, 7)
                ct = cm(m+1)
                am(m+1) = am(m+1) + cm(m+1)
                itf = m + 1
            case (8:10)
                itf = m
                isp = 1
                ct = cm(m)
        end select

        !
        ! compute homogeneous solution with solution at pole equal to one
        !
        itsp = its + 1
        itfm = itf - 1
        wts = sint(its+1)*am(its+1)/cm(its)
        wtf = sint(itf-1)*cm(itf-1)/am(itf)
        munk = itf - its + 1

        if (isp > 0) then
            d(its) = cm(its)/bm(its)
            do i = itsp, m
                d(i) = cm(i)/(bm(i)-am(i)*d(i-1))
            end do
            ss(m) = -d(m)
            iid = m - its
            do ii = 1, iid
                i = m - ii
                ss(i) = -d(i)*ss(i+1)
            end do
            ss(m+1) = ONE
        end if

        if (inp > 0) then
            sn(1) = ONE
            d(itf) = am(itf)/bm(itf)
            iid = itf - 2
            do ii = 1, iid
                i = itf - ii
                d(i) = am(i)/(bm(i)-cm(i)*d(i+1))
            end do
            sn(2) = -d(2)
            do i = 3, itf
                sn(i) = -d(i)*sn(i-1)
            end do
        end if
        !
        ! boundary conditions at phi=ps
        !
        nbr = nbdcnd + 1
        wps = ONE
        wpf = ONE
        select case (nbr)
            case default
                jps = 1
            case (2:3)
                jps = 2
            case (4:5)
                jps = 1
                wps = HALF
        end select
        !
        ! boundary condition at phi=pf
        !
        select case (nbr)
            case (1)
                jpf = n
            case (2, 5)
                jpf = n
            case (3:4)
                wpf = HALF
                jpf = n + 1
        end select

        jpsp = jps + 1
        jpfm = jpf - 1
        nunk = jpf - jps + 1
        fjj = jpfm - jpsp + 1
        !
        ! scale coefficients for subroutine genbun
        !
        do i = its, itf
            cf = dphi2*sint(i)*sint(i)
            am(i) = cf*am(i)
            bm(i) = cf*bm(i)
            cm(i) = cf*cm(i)
        end do

        am(its) = ZERO
        cm(itf) = ZERO
        ising = 0

        select case (mbr)
            case (1, 4, 7, 9:10)
                select case (nbr)
                    case (1, 4)
                        if (elmbda >= ZERO) then
                            ising = 1
                            summation = wts*wps + wts*wpf + wtf*wps + wtf*wpf

                            if (inp > 0) summation = summation + wpp

                            if (isp > 0) summation = summation + wpp

                            sum1 = sum(sint(itsp:itfm))
                            summation = summation + fjj*(sum1 + wts + wtf)
                            summation = summation + (wps + wpf)*sum1
                            hne = summation
                        end if
                end select
        end select

        select case (mbr)
            case (2:3, 8)
                f(2,jps:jpf) = f(2,jps:jpf) - at*f(1,jps:jpf)
            case (4:5, 9)
                f(1,jps:jpf) = f(1,jps:jpf) + two_dth*bdts(jps:jpf)*at
            case (6:7, 10)
                if (nbdcnd == 3) then
                    yhld = f(1,jps) - 4.0_wp/(fn*dphi*dth2)*(bdpf(2)-bdps(2))
                    f(1,:np1) = yhld
                end if
        end select

        select case (mbr)
            case (2, 5:6)
                f(m,jps:jpf) = f(m,jps:jpf) - ct*f(m+1,jps:jpf)
            case (3:4, 7)
                f(m+1,jps:jpf) = f(m+1,jps:jpf) - two_dth*bdtf(jps:jpf)*ct
            case (8:10)
                if (nbdcnd == 3) then
                    yhld = f(m+1,jps) - 4.0_wp/(fn*dphi*dth2)*(bdpf(m)-bdps(m))
                    f(m+1,:np1) = yhld
                end if
        end select

        select case (nbr)
            case (2:3)
                f(its:itf,2) = f(its:itf,2) &
                    - f(its:itf,1)/(dphi2*sint(its:itf)*sint(its:itf))
            case (4:5)
                f(its:itf,1) = f(its:itf,1) &
                    + two_dphi*bdps(its:itf)/(dphi2*sint(its:itf)*sint(its:itf))
        end select

        select case (nbr)
            case (2, 5)
                f(its:itf,n) = f(its:itf,n) &
                    - f(its:itf,n+1)/(dphi2*sint(its:itf)*sint(its:itf))
            case (3:4)
                f(its:itf,n+1) = f(its:itf,n+1) &
                    - two_dphi*bdpf(its:itf)/(dphi2*sint(its:itf)*sint(its:itf))
        end select

        pertrb = ZERO

        if (ising /= 0) then
            summation = &
                wts*wps*f(its,jps) + wts*wpf*f(its,jpf) &
                + wtf*wps*f(itf,jps) + wtf*wpf*f(itf,jpf)

            if (inp > 0) summation = summation + wpp*f(1,jps)
            if (isp > 0) summation = summation + wpp*f(m+1,jps)

            do i = itsp, itfm
                summation = summation + sint(i)*sum(f(i,jpsp:jpfm))
            end do

            sum1 = sum(f(its,jpsp:jpfm))
            sum2 = sum(f(itf,jpsp:jpfm))
            summation = summation + wts*sum1 + wtf*sum2
            sum1 = ZERO
            sum2 = ZERO
            sum1 = dot_product(sint(itsp:itfm),f(itsp:itfm,jps))
            sum2 = dot_product(sint(itsp:itfm),f(itsp:itfm,jpf))
            summation = summation + wps*sum1 + wpf*sum2
            pertrb = summation/hne
            f(:mp1,:np1) = f(:mp1,:np1) - pertrb
        end if

        ! scale right side for subroutine genbun_lower_routine
        do i = its, itf
            cf = dphi2*sint(i)*sint(i)
            f(i,jps:jpf) = cf*f(i,jps:jpf)
        end do

        call_arguments: associate( &
            a_arg => am(its:), &
            b_arg => bm(its:), &
            c_arg => cm(its:), &
            y_arg => f(its:,jps:) &
            )

            ! Invoke genbun_lower_routine solver
            call util%genbun_lower_routine(nbdcnd, nunk, 1, munk, a_arg, b_arg, c_arg, &
                idimf, y_arg, error_flag, d)
        end associate call_arguments

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

        if (ising > 0 .and. inp > 0 .and. isp <= 0) then
            f(1,:np1) = ZERO
        else if (isp > 0) then
            f(m+1,:np1) = ZERO
        else if (inp > 0) then

            summation = wps*f(its,jps) + wpf*f(its,jpf) + sum(f(its,jpsp:jpfm))
            dfn = cp*summation
            dnn = cp*((wps + wpf + fjj)*(sn(2)-ONE)) + elmbda
            dsn = cp*(wps + wpf + fjj)*sn(m)

            if (isp > 0) then
                cnp = (f(1,1)-dfn)/dnn
                do i = its, itf
                    hld = cnp*sn(i)
                    f(i,jps:jpf) = f(i,jps:jpf) + hld
                end do
                f(1,:np1) = cnp
            end if
        else
            summation = wps*f(itf,jps) + wpf*f(itf,jpf) + sum(f(itf,jpsp:jpfm))
            dfs = cp*summation
            dss = cp*((wps + wpf + fjj)*(ss(m)-ONE)) + elmbda
            dns = cp*(wps + wpf + fjj)*ss(2)

            if (inp <= 0) then
                csp = (f(m+1,1)-dfs)/dss
                do i = its, itf
                    hld = csp*ss(i)
                    f(i,jps:jpf) = f(i,jps:jpf) + hld
                end do
                f(m+1,:np1) = csp
            else
                rtn = f(1,1) - dfn
                rts = f(m+1,1) - dfs
                if (ising > 0) then
                    csp = ZERO
                    cnp = rtn/dnn
                else
                    if (abs(dnn) - abs(dsn) > ZERO) then
                        den = dss - dns*dsn/dnn
                        rts = rts - rtn*dsn/dnn
                        csp = rts/den
                        cnp = (rtn - csp*dns)/dnn
                    else
                        den = dns - dss*dnn/dsn
                        rtn = rtn - rts*dnn/dsn
                        csp = rtn/den
                        cnp = (rts - dss*csp)/dsn
                    end if
                end if
                do i = its, itf
                    hld = cnp*sn(i) + csp*ss(i)
                    f(i,jps:jpf) = f(i,jps:jpf) + hld
                end do
                f(1,:np1) = cnp
                f(m+1,:np1) = csp
            end if
        end if

        if (nbdcnd == 0) f(:mp1,jpf+1) = f(:mp1,jps)

    end subroutine hwsssp_lower_routine

end submodule centered_spherical_solver

max_stars_repo_name: jlokimlin/fishpack

max_issues_repo_name: jlokimlin/fishpack


text:
module temperature_test
    use double_precision_generator_m, only: DOUBLE_PRECISION_GENERATOR
    use double_precision_pair_generator_m, only: DOUBLE_PRECISION_PAIR_GENERATOR
    use double_precision_pair_input_m, only: double_precision_pair_input_t
    use erloff, only: error_list_t
    use iso_varying_string, only: operator(//)
    use non_zero_double_precision_generator_m, only: &
            NON_ZERO_DOUBLE_PRECISION_GENERATOR
    use non_zero_double_precision_pair_generator_m, only: &
            NON_ZERO_DOUBLE_PRECISION_PAIR_GENERATOR
    use quaff, only: &
            temperature_t, &
            fallible_temperature_t, &
            temperature_unit_t, &
            operator(.unit.), &
            parse_temperature, &
            sum, &
            PROVIDED_TEMPERATURE_UNITS, &
            KELVIN
    use quaff_asserts_m, only: assert_equals, assert_equals_within_relative
    use quaff_utilities_m, only: PARSE_ERROR
    use temperature_utilities_m, only: &
            units_input_t, units_pair_input_t, make_units_examples
    use units_examples_m, only: units_examples_t
    use veggies, only: &
            double_precision_input_t, &
            input_t, &
            result_t, &
            test_item_t, &
            test_result_item_t, &
            assert_equals, &
            assert_equals_within_relative, &
            assert_not, &
            assert_that, &
            describe, &
            fail, &
            it

    implicit none
    private
    public :: test_temperature
contains
    function test_temperature() result(tests)
        type(test_item_t) :: tests

        type(units_examples_t) :: examples

        examples = make_units_examples(PROVIDED_TEMPERATURE_UNITS)
        tests = describe( &
                "temperature_t", &
                [ it( &
                        "returns the same value given the same units", &
                        examples%units(), &
                        check_round_trip) &
                , it( &
                        "preserves its value converting to and from a string", &
                        examples%units(), &
                        check_to_and_from_string) &
                , it( &
                        "returns an error trying to parse a bad string", &
                        check_bad_string) &
                , it( &
                        "returns an error trying to parse an unknown unit", &
                        check_bad_unit) &
                , it( &
                        "returns an error trying to parse a bad number", &
                        check_bad_number) &
                , it( &
                        "multiplying by one returns the original temperature", &
                        DOUBLE_PRECISION_GENERATOR, &
                        check_multiply_by_one) &
                , it( &
                        "multiplying by zero returns zero temperature", &
                        DOUBLE_PRECISION_GENERATOR, &
                        check_multiply_by_zero) &
                , it( &
                        "dividing by one returns the original temperature", &
                        DOUBLE_PRECISION_GENERATOR, &
                        check_divide_by_one) &
                , it( &
                        "dividing by itself returns one", &
                        NON_ZERO_DOUBLE_PRECISION_GENERATOR, &
                        check_divide_by_self) &
                , it( &
                        "multiplying and dividing by the same number returns the original temperature", &
                        NON_ZERO_DOUBLE_PRECISION_PAIR_GENERATOR, &
                        check_multiply_divide) &
                , describe( &
                        "operator(==)", &
                        [ it( &
                                "is true for the same temperature", &
                                DOUBLE_PRECISION_GENERATOR, &
                                check_equal_with_same_number) &
                        , it( &
                                "is false for different temperatures", &
                                DOUBLE_PRECISION_GENERATOR, &
                                check_equal_with_different_numbers) &
                        ]) &
                , describe( &
                        "operator(/=)", &
                        [ it( &
                                "is false for the same temperature", &
                                DOUBLE_PRECISION_GENERATOR, &
                                check_not_equal_with_same_number) &
                        , it( &
                                "is true for different temperatures", &
                                DOUBLE_PRECISION_GENERATOR, &
                                check_not_equal_with_different_numbers) &
                        ]) &
                , describe( &
                        "%equal(temperature, within)", &
                        [ it( &
                                "is true for the same temperature even for tiny tolerance", &
                                DOUBLE_PRECISION_GENERATOR, &
                                check_equal_within_with_same_number) &
                        , it( &
                                "is true for sufficiently close values", &
                                DOUBLE_PRECISION_GENERATOR, &
                                check_equal_within_with_close_numbers) &
                        , it( &
                                "is false for sufficiently different numbers", &
                                DOUBLE_PRECISION_GENERATOR, &
                                check_equal_within_with_different_numbers) &
                        ]) &
                , describe( &
                        "operator(>=)", &
                        [ it( &
                                "is true if the lhs is greater than the rhs", &
                                DOUBLE_PRECISION_GENERATOR, &
                                check_greater_than_or_equal_with_greater_number) &
                        , it( &
                                "is true if the lhs is equal to the rhs", &
                                DOUBLE_PRECISION_GENERATOR, &
                                check_greater_than_or_equal_with_same_numbers) &
                        , it( &
                                "is false if the lhs is less than the rhs", &
                                DOUBLE_PRECISION_GENERATOR, &
                                check_greater_than_or_equal_with_lesser_number) &
                        ]) &
                , describe( &
                        "operator(<=)", &
                        [ it( &
                                "is true if the lhs is less than the rhs", &
                                DOUBLE_PRECISION_GENERATOR, &
                                check_less_than_or_equal_with_less_number) &
                        , it( &
                                "is true if the lhs is equal to the rhs", &
                                DOUBLE_PRECISION_GENERATOR, &
                                check_less_than_or_equal_with_same_numbers) &
                        , it( &
                                "is false if the lhs is greater than the rhs", &
                                DOUBLE_PRECISION_GENERATOR, &
                                check_less_than_or_equal_with_greater_number) &
                        ]) &
                , describe( &
                        "operator(>)", &
                        [ it( &
                                "is true if the lhs is greater than the rhs", &
                                DOUBLE_PRECISION_GENERATOR, &
                                check_greater_than_with_greater_number) &
                        , it( &
                                "is false if the lhs is equal to the rhs", &
                                DOUBLE_PRECISION_GENERATOR, &
                                check_greater_than_with_same_numbers) &
                        , it( &
                                "is false if the lhs is less than the rhs", &
                                DOUBLE_PRECISION_GENERATOR, &
                                check_greater_than_with_lesser_number) &
                        ]) &
                , describe( &
                        "operator(<)", &
                        [ it( &
                                "is true if the lhs is less than the rhs", &
                                DOUBLE_PRECISION_GENERATOR, &
                                check_less_than_with_less_number) &
                        , it( &
                                "is false if the lhs is equal to the rhs", &
                                DOUBLE_PRECISION_GENERATOR, &
                                check_less_than_with_same_numbers) &
                        , it( &
                                "is false if the lhs is greater than the rhs", &
                                DOUBLE_PRECISION_GENERATOR, &
                                check_less_than_with_greater_number) &
                        ]) &
                ])
    end function

    function check_round_trip(input) result(result_)
        class(input_t), intent(in) :: input
        type(result_t) :: result_

        select type (input)
        type is (units_input_t)
            result_ = check_round_trip_in(input%unit())
        class default
            result_ = fail("Expected to get a units_input_t")
        end select
    end function

    function check_round_trip_in(units) result(result_)
        class(temperature_unit_t), intent(in) :: units
        type(result_t) :: result_

        type(test_item_t) :: the_test
        type(test_result_item_t) :: the_result

        the_test = it(units%to_string(), DOUBLE_PRECISION_GENERATOR, check_round_trip_)
        the_result = the_test%run()
        result_ = assert_that(the_result%passed(), the_result%verbose_description(.false.))
    contains
        pure function check_round_trip_(input) result(result__)
            class(input_t), intent(in) :: input
            type(result_t) :: result__

            type(temperature_t) :: intermediate

            select type (input)
            type is (double_precision_input_t)
                intermediate = input%input().unit.units
                result__ = assert_equals_within_relative( &
                        input%input(), &
                        intermediate.in.units, &
                        1.0d-12)
            class default
                result__ = fail("Expected to get a double_precision_input_t")
            end select
        end function
    end function

    function check_to_and_from_string(input) result(result_)
        class(input_t), intent(in) :: input
        type(result_t) :: result_

        select type (input)
        type is (units_input_t)
            result_ = check_string_trip(input%unit())
        class default
            result_ = fail("Expected to get an units_input_t")
        end select
    end function

    function check_string_trip(units) result(result_)
        class(temperature_unit_t), intent(in) :: units
        type(result_t) :: result_

        type(test_item_t) :: the_test
        type(test_result_item_t) :: the_result

        the_test = it(units%to_string(), DOUBLE_PRECISION_GENERATOR, do_check)
        the_result = the_test%run()
        result_ = assert_that(the_result%passed(), the_result%verbose_description(.false.))
    contains
        function do_check(input) result(result__)
            class(input_t), intent(in) :: input
            type(result_t) :: result__

            type(error_list_t) :: errors
            type(temperature_t) :: original_temperature
            type(fallible_temperature_t) :: maybe_temperature
            type(temperature_t) :: new_temperature

            select type (input)
            type is (double_precision_input_t)
                original_temperature = input%input().unit.units
                maybe_temperature = parse_temperature( &
                        original_temperature%to_string_in(units))
                new_temperature = maybe_temperature%temperature()
                errors = maybe_temperature%errors()
                result__ = &
                        assert_equals( &
                                original_temperature, &
                                new_temperature) &
                        .and.assert_not(errors%has_any(), errors%to_string())
            class default
                result__ = fail("Expected to get a double_precision_input_t")
            end select
        end function
    end function

    function check_bad_string() result(result_)
        type(result_t) :: result_

        type(error_list_t) :: errors
        type(fallible_temperature_t) :: maybe_temperature

        maybe_temperature = parse_temperature("bad")
        errors = maybe_temperature%errors()
        result_ = assert_that(errors.hasType.PARSE_ERROR, errors%to_string())
    end function

    function check_bad_unit() result(result_)
        type(result_t) :: result_

        type(error_list_t) :: errors
        type(fallible_temperature_t) :: maybe_temperature

        maybe_temperature = parse_temperature("1.0 Kbad", [KELVIN])
        errors = maybe_temperature%errors()
        result_ = assert_that(errors.hasType.PARSE_ERROR, errors%to_string())
    end function

    function check_bad_number() result(result_)
        type(result_t) :: result_

        type(error_list_t) :: errors
        type(fallible_temperature_t) :: maybe_temperature

        maybe_temperature = parse_temperature("bad K")
        errors = maybe_temperature%errors()
        result_ = assert_that(errors.hasType.PARSE_ERROR, errors%to_string())
    end function

    pure function check_multiply_by_one(input) result(result_)
        class(input_t), intent(in) :: input
        type(result_t) :: result_

        type(temperature_t) :: temperature

        select type(input)
        type is (double_precision_input_t)
            temperature = input%input().unit.KELVIN
            result_ = assert_equals(temperature, temperature * 1.0d0)
        class default
            result_ = fail("Expected a double_precision_input_t")
        end select
    end function

    pure function check_multiply_by_zero(input) result(result_)
        class(input_t), intent(in) :: input
        type(result_t) :: result_

        type(temperature_t) :: temperature
        type(temperature_t) :: zero

        select type(input)
        type is (double_precision_input_t)
            temperature = input%input().unit.KELVIN
            zero = 0.0d0.unit.KELVIN
            result_ = assert_equals(zero, temperature * 0.0d0)
        class default
            result_ = fail("Expected a double_precision_input_t")
        end select
    end function

    pure function check_divide_by_one(input) result(result_)
        class(input_t), intent(in) :: input
        type(result_t) :: result_

        type(temperature_t) :: temperature

        select type(input)
        type is (double_precision_input_t)
            temperature = input%input().unit.KELVIN
            result_ = assert_equals(temperature, temperature / 1.0d0)
        class default
            result_ = fail("Expected a double_precision_input_t")
        end select
    end function

    pure function check_divide_by_self(input) result(result_)
        class(input_t), intent(in) :: input
        type(result_t) :: result_

        type(temperature_t) :: temperature

        select type(input)
        type is (double_precision_input_t)
            temperature = input%input().unit.KELVIN
            result_ = assert_equals(1.0d0, temperature / temperature)
        class default
            result_ = fail("Expected a double_precision_input_t")
        end select
    end function

    pure function check_multiply_divide(input) result(result_)
        class(input_t), intent(in) :: input
        type(result_t) :: result_

        type(temperature_t) :: temperature

        select type (input)
        type is (double_precision_pair_input_t)
            temperature = input%first().unit.KELVIN
            result_ = assert_equals(temperature, (temperature * input%second_()) / input%second_())
        class default
            result_ = fail("Expected a double_precision_pair_input_t")
        end select
    end function

    pure function check_equal_with_same_number(input) result(result_)
        class(input_t), intent(in) :: input
        type(result_t) :: result_

        type(temperature_t) :: the_temperature

        select type (input)
        type is (double_precision_input_t)
            the_temperature = input%input().unit.KELVIN
            result_ = assert_that( &
                    the_temperature == the_temperature, &
                    the_temperature%to_string() // " == " // the_temperature%to_string())
        class default
            result_ = fail("Expected to get a double_precision_input_t")
        end select
    end function

    pure function check_equal_with_different_numbers(input) result(result_)
        class(input_t), intent(in) :: input
        type(result_t) :: result_

        type(temperature_t) :: temperature1
        type(temperature_t) :: temperature2

        select type (input)
        type is (double_precision_input_t)
            temperature1 = input%input().unit.KELVIN
            temperature2 = (input%input() + 1.0d0).unit.KELVIN
            result_ = assert_not( &
                    temperature1 == temperature2, &
                    temperature1%to_string() // " == " // temperature2%to_string())
        class default
            result_ = fail("Expected to get a double_precision_input_t")
        end select
    end function

    pure function check_not_equal_with_same_number(input) result(result_)
        class(input_t), intent(in) :: input
        type(result_t) :: result_

        type(temperature_t) :: the_temperature

        select type (input)
        type is (double_precision_input_t)
            the_temperature = input%input().unit.KELVIN
            result_ = assert_not( &
                    the_temperature /= the_temperature, &
                    the_temperature%to_string() // " /= " // the_temperature%to_string())
        class default
            result_ = fail("Expected to get a double_precision_input_t")
        end select
    end function

    pure function check_not_equal_with_different_numbers(input) result(result_)
        class(input_t), intent(in) :: input
        type(result_t) :: result_

        type(temperature_t) :: temperature1
        type(temperature_t) :: temperature2

        select type (input)
        type is (double_precision_input_t)
            temperature1 = input%input().unit.KELVIN
            temperature2 = (input%input() + 1.0d0).unit.KELVIN
            result_ = assert_that( &
                    temperature1 /= temperature2, &
                    temperature1%to_string() // " /= " // temperature2%to_string())
        class default
            result_ = fail("Expected to get a double_precision_input_t")
        end select
    end function

    pure function check_equal_within_with_same_number(input) result(result_)
        class(input_t), intent(in) :: input
        type(result_t) :: result_

        type(temperature_t) :: the_temperature
        type(temperature_t) :: tolerance

        select type (input)
        type is (double_precision_input_t)
            the_temperature = input%input().unit.KELVIN
            tolerance = tiny(1.0d0).unit.KELVIN
            result_ = assert_that( &
                    the_temperature%equal(the_temperature, within = tolerance), &
                    "(" // the_temperature%to_string() // ")%equal(" &
                        // the_temperature%to_string() // ", within = " &
                        // tolerance%to_string() // ")")
        class default
            result_ = fail("Expected to get a double_precision_input_t")
        end select
    end function

    pure function check_equal_within_with_close_numbers(input) result(result_)
        class(input_t), intent(in) :: input
        type(result_t) :: result_

        type(temperature_t) :: temperature1
        type(temperature_t) :: temperature2
        type(temperature_t) :: tolerance

        select type (input)
        type is (double_precision_input_t)
            temperature1 = input%input().unit.KELVIN
            temperature2 = (input%input() + 0.05d0).unit.KELVIN
            tolerance = 0.1d0.unit.KELVIN
            result_ = assert_that( &
                    temperature1%equal(temperature2, within = tolerance), &
                    "(" // temperature1%to_string() // ")%equal(" &
                        // temperature2%to_string() // ", within = " &
                        // tolerance%to_string() // ")")
        class default
            result_ = fail("Expected to get a double_precision_input_t")
        end select
    end function

    pure function check_equal_within_with_different_numbers(input) result(result_)
        class(input_t), intent(in) :: input
        type(result_t) :: result_

        type(temperature_t) :: temperature1
        type(temperature_t) :: temperature2
        type(temperature_t) :: tolerance

        select type (input)
        type is (double_precision_input_t)
            temperature1 = input%input().unit.KELVIN
            temperature2 = (input%input() + 0.2d0).unit.KELVIN
            tolerance = 0.1d0.unit.KELVIN
            result_ = assert_not( &
                    temperature1%equal(temperature2, within = tolerance), &
                    "(" // temperature1%to_string() // ")%equal(" &
                    // temperature2%to_string() // ", within = " &
                    // tolerance%to_string() // ")")
        class default
            result_ = fail("Expected to get a double_precision_input_t")
        end select
    end function

    pure function check_greater_than_or_equal_with_greater_number(input) result(result_)
        class(input_t), intent(in) :: input
        type(result_t) :: result_

        type(temperature_t) :: temperature1
        type(temperature_t) :: temperature2

        select type (input)
        type is (double_precision_input_t)
            temperature1 = input%input().unit. KELVIN
            temperature2 = (input%input() - 1.0d0).unit.KELVIN
            result_ = assert_that( &
                    temperature1 >= temperature2, &
                    temperature1%to_string() // " >= " // temperature2%to_string())
        class default
            result_ = fail("Expected to get a double_precision_input_t")
        end select
    end function

    pure function check_greater_than_or_equal_with_same_numbers(input) result(result_)
        class(input_t), intent(in) :: input
        type(result_t) :: result_

        type(temperature_t) :: temperature1
        type(temperature_t) :: temperature2

        select type (input)
        type is (double_precision_input_t)
            temperature1 = input%input().unit.KELVIN
            temperature2 = input%input().unit.KELVIN
            result_ = assert_that( &
                    temperature1 >= temperature2, &
                    temperature1%to_string() // " >= " // temperature2%to_string())
        class default
            result_ = fail("Expected to get a double_precision_input_t")
        end select
    end function

    pure function check_greater_than_or_equal_with_lesser_number(input) result(result_)
        class(input_t), intent(in) :: input
        type(result_t) :: result_

        type(temperature_t) :: temperature1
        type(temperature_t) :: temperature2

        select type (input)
        type is (double_precision_input_t)
            temperature1 = input%input().unit.KELVIN
            temperature2 = (input%input() + 1.0d0).unit.KELVIN
            result_ = assert_not( &
                    temperature1 >= temperature2, &
                    temperature1%to_string() // " >= " // temperature2%to_string())
        class default
            result_ = fail("Expected to get a double_precision_input_t")
        end select
    end function

    pure function check_less_than_or_equal_with_less_number(input) result(result_)
        class(input_t), intent(in) :: input
        type(result_t) :: result_

        type(temperature_t) :: temperature1
        type(temperature_t) :: temperature2

        select type (input)
        type is (double_precision_input_t)
            temperature1 = input%input().unit. KELVIN
            temperature2 = (input%input() + 1.0d0).unit.KELVIN
            result_ = assert_that( &
                    temperature1 <= temperature2, &
                    temperature1%to_string() // " <= " // temperature2%to_string())
        class default
            result_ = fail("Expected to get a double_precision_input_t")
        end select
    end function

    pure function check_less_than_or_equal_with_same_numbers(input) result(result_)
        class(input_t), intent(in) :: input
        type(result_t) :: result_

        type(temperature_t) :: temperature1
        type(temperature_t) :: temperature2

        select type (input)
        type is (double_precision_input_t)
            temperature1 = input%input().unit.KELVIN
            temperature2 = input%input().unit.KELVIN
            result_ = assert_that( &
                    temperature1 <= temperature2, &
                    temperature1%to_string() // " <= " // temperature2%to_string())
        class default
            result_ = fail("Expected to get a double_precision_input_t")
        end select
    end function

    pure function check_less_than_or_equal_with_greater_number(input) result(result_)
        class(input_t), intent(in) :: input
        type(result_t) :: result_

        type(temperature_t) :: temperature1
        type(temperature_t) :: temperature2

        select type (input)
        type is (double_precision_input_t)
            temperature1 = input%input().unit.KELVIN
            temperature2 = (input%input() - 1.0d0).unit.KELVIN
            result_ = assert_not( &
                    temperature1 <= temperature2, &
                    temperature1%to_string() // " <= " // temperature2%to_string())
        class default
            result_ = fail("Expected to get a double_precision_input_t")
        end select
    end function

    pure function check_greater_than_with_greater_number(input) result(result_)
        class(input_t), intent(in) :: input
        type(result_t) :: result_

        type(temperature_t) :: temperature1
        type(temperature_t) :: temperature2

        select type (input)
        type is (double_precision_input_t)
            temperature1 = input%input().unit. KELVIN
            temperature2 = (input%input() - 1.0d0).unit.KELVIN
            result_ = assert_that( &
                    temperature1 > temperature2, &
                    temperature1%to_string() // " > " // temperature2%to_string())
        class default
            result_ = fail("Expected to get a double_precision_input_t")
        end select
    end function

    pure function check_greater_than_with_same_numbers(input) result(result_)
        class(input_t), intent(in) :: input
        type(result_t) :: result_

        type(temperature_t) :: temperature1
        type(temperature_t) :: temperature2

        select type (input)
        type is (double_precision_input_t)
            temperature1 = input%input().unit.KELVIN
            temperature2 = input%input().unit.KELVIN
            result_ = assert_not( &
                    temperature1 > temperature2, &
                    temperature1%to_string() // " > " // temperature2%to_string())
        class default
            result_ = fail("Expected to get a double_precision_input_t")
        end select
    end function

    pure function check_greater_than_with_lesser_number(input) result(result_)
        class(input_t), intent(in) :: input
        type(result_t) :: result_

        type(temperature_t) :: temperature1
        type(temperature_t) :: temperature2

        select type (input)
        type is (double_precision_input_t)
            temperature1 = input%input().unit.KELVIN
            temperature2 = (input%input() + 1.0d0).unit.KELVIN
            result_ = assert_not( &
                    temperature1 > temperature2, &
                    temperature1%to_string() // " > " // temperature2%to_string())
        class default
            result_ = fail("Expected to get a double_precision_input_t")
        end select
    end function

    pure function check_less_than_with_less_number(input) result(result_)
        class(input_t), intent(in) :: input
        type(result_t) :: result_

        type(temperature_t) :: temperature1
        type(temperature_t) :: temperature2

        select type (input)
        type is (double_precision_input_t)
            temperature1 = input%input().unit. KELVIN
            temperature2 = (input%input() + 1.0d0).unit.KELVIN
            result_ = assert_that( &
                    temperature1 < temperature2, &
                    temperature1%to_string() // " < " // temperature2%to_string())
        class default
            result_ = fail("Expected to get a double_precision_input_t")
        end select
    end function

    pure function check_less_than_with_same_numbers(input) result(result_)
        class(input_t), intent(in) :: input
        type(result_t) :: result_

        type(temperature_t) :: temperature1
        type(temperature_t) :: temperature2

        select type (input)
        type is (double_precision_input_t)
            temperature1 = input%input().unit.KELVIN
            temperature2 = input%input().unit.KELVIN
            result_ = assert_not( &
                    temperature1 < temperature2, &
                    temperature1%to_string() // " <=" // temperature2%to_string())
        class default
            result_ = fail("Expected to get a double_precision_input_t")
        end select
    end function

    pure function check_less_than_with_greater_number(input) result(result_)
        class(input_t), intent(in) :: input
        type(result_t) :: result_

        type(temperature_t) :: temperature1
        type(temperature_t) :: temperature2

        select type (input)
        type is (double_precision_input_t)
            temperature1 = input%input().unit.KELVIN
            temperature2 = (input%input() - 1.0d0).unit.KELVIN
            result_ = assert_not( &
                    temperature1 < temperature2, &
                    temperature1%to_string() // " < " // temperature2%to_string())
        class default
            result_ = fail("Expected to get a double_precision_input_t")
        end select
    end function
end module

max_stars_repo_name: everythingfunctional/quaff

max_issues_repo_name: everythingfunctional/quaff

max_forks_repo_name: everythingfunctional/quaff


text:
module ek_eigenpairs_types_m
  type ek_eigenpairs_local_t ! type number = 1
    double precision, allocatable :: values(:)
    double precision, allocatable :: vectors(:, :)
  end type ek_eigenpairs_local_t

  type ek_eigenpairs_blacs_t ! type number = 2
    double precision, allocatable :: values(:)
    integer :: desc(9)
    double precision, allocatable :: Vectors(:, :)
  end type ek_eigenpairs_blacs_t

  type ek_eigenpairs_types_union_t
    integer :: type_number
    type(ek_eigenpairs_local_t) :: local
    type(ek_eigenpairs_blacs_t) :: blacs
  end type ek_eigenpairs_types_union_t
end module ek_eigenpairs_types_m

max_stars_repo_name: eigenkernel/EigenKernel

max_issues_repo_name: eigenkernel/EigenKernel

max_forks_repo_name: eigenkernel/EigenKernel


text:
subroutine tem_reabcs(a,itask,kflag)
   use Mod_Memor
   use typre

   use Mod_Temperature
   implicit none
   class(TemperatureProblem) :: a
   integer(ip) :: itask
   integer(ip), optional :: kflag

   integer(ip) :: npoin

   !Initializations
   if (itask == 0) then
      if(a%kfl_sourc==2) then
         call a%Mesh%GetNpoin(npoin)
         call a%Memor%alloc(npoin,a%PointwiseSource,'PointwiseSource','tem_reabcs')
      end if

   !Header
   elseif(itask == 1) then
      if(a%Listener%exists('WALLD')) then
         a%delta=a%Listener%getrea('WALLD',0.0_rp,'#Distance to the wall')
      end if

   !Finalization
   elseif(itask == 100) then

   endif

end subroutine

subroutine tem_ReadOnNodes(a)
   use typre
   use Mod_Mesh
   use Mod_MPIObject
   use Mod_Memor
   use Mod_Listen
   use Mod_PhysicalProblem
   use Mod_Temperature
   implicit none
   class(TemperatureProblem) :: a

   !On Nodes
   if(a%kfl_conbc==0) then
      if(a%kfl_sourc==2) then
         a%PointwiseSource(a%gipoin) = a%Listener%param(5)
      end if

   else
      if(a%kfl_sourc==2) a%PointwiseSource(a%gipoin) = a%Listener%param(4)
   end if
end subroutine

subroutine tem_ReadOnBoundaries(a)
   use typre
   use Mod_Mesh
   use Mod_MPIObject
   use Mod_Memor
   use Mod_Listen
   use Mod_PhysicalProblem
   use Mod_Temperature
   implicit none
   class(TemperatureProblem) :: a

   !Robin
   if(a%kfl_fixbo(a%giboun)==3) then

      if(a%kfl_conbc==0) a%kfl_funbo(a%giboun) = int(a%Listener%param(a%gipsta+4))
      if(a%kfl_funbo(a%giboun)/=0) then
         allocate(a%bvnat(a%giboun)%a(6))

         a%bvnat(a%giboun)%a(1)=a%Listener%param(a%gipsta+1)
         a%bvnat(a%giboun)%a(2)=a%Listener%param(a%gipsta+2)
         a%bvnat(a%giboun)%a(3)=a%Listener%param(a%gipsta+3)
         a%bvnat(a%giboun)%a(4)=a%Listener%param(a%gipsta+1)
         a%bvnat(a%giboun)%a(5)=a%Listener%param(a%gipsta+2)
         a%bvnat(a%giboun)%a(6)=a%Listener%param(a%gipsta+3)
      else
         allocate(a%bvnat(a%giboun)%a(3))
         a%bvnat(a%giboun)%a(1)=a%Listener%param(a%gipsta+1)
         a%bvnat(a%giboun)%a(2)=a%Listener%param(a%gipsta+2)
         a%bvnat(a%giboun)%a(3)=a%Listener%param(a%gipsta+3)
      end if
   !Wall law
   else if(a%kfl_fixbo(a%giboun)==4) then

      if(a%kfl_conbc==0) a%kfl_funbo(a%giboun) = int(a%Listener%param(a%gipsta+2))
      if(a%kfl_funbo(a%giboun)/=0) then
         allocate(a%bvnat(a%giboun)%a(2))
         a%bvnat(a%giboun)%a(1)=a%Listener%param(a%gipsta+1)
         a%bvnat(a%giboun)%a(2)=a%Listener%param(a%gipsta+1)
      else
         allocate(a%bvnat(a%giboun)%a(1))
         a%bvnat(a%giboun)%a(1)=a%Listener%param(a%gipsta+1)
      end if
   !Dust Transport
   elseif (a%kfl_fixbo(a%giboun) == 5) then
      !Nothing to be done here

   !Do nothing boundary condition
   elseif (a%kfl_fixbo(a%giboun) == 6) then
      !Do nothing boundary condition
      !Nothing to be done here

   end if

end subroutine

subroutine tem_ReadOnElements(a)
   use typre
   use Mod_Listen
   use Mod_Temperature
   implicit none
   class(TemperatureProblem) :: a

   integer(ip) :: ielem,imaterial

   ielem = a%gielem
   imaterial = int(a%Listener%param(a%gipsta+1))
   if (imaterial <= 0 .or. imaterial > a%NumberOfMaterials) call runend('Wrong material number')
   a%ElementMaterials(ielem) = imaterial
end subroutine

max_stars_repo_name: ciaid-colombia/InsFEM

max_issues_repo_name: ciaid-colombia/InsFEM

max_forks_repo_name: ciaid-colombia/InsFEM


text:
program test_arange
use iso_fortran_env, only:int8, int16, int32, int64, real32, real64, real128
use testing, only:assert, initialize_tests, report_tests
use functional, only:arange

implicit none

logical, dimension(:), allocatable :: tests
logical :: test_failed
integer :: n, ntests

n = 1
ntests = 19
call initialize_tests(tests, ntests)

tests(n) = assert(all(arange(1_int8, 3_int8) == [1_int8, 2_int8, 3_int8]), &
                  'arange,  int8')
n = n + 1

tests(n) = assert(all(arange(1_int16, 3_int16) == [1_int16, 2_int16, 3_int16]), &
                  'arange,  int16')
n = n + 1

tests(n) = assert(all(arange(1_int32, 3_int32) == [1_int32, 2_int32, 3_int32]), &
                  'arange,  int32')
n = n + 1

tests(n) = assert(all(arange(1_int64, 3_int64) == [1_int64, 2_int64, 3_int64]), &
                  'arange,  int64')
n = n + 1

tests(n) = assert(all(arange(1._real32, 3._real32) == [1._real32, 2._real32, 3._real32]), &
                  'arange,  real32')
n = n + 1

tests(n) = assert(all(arange(1._real64, 3._real64) == [1._real64, 2._real64, 3._real64]), &
                  'arange,  real32')
n = n + 1

tests(n) = assert(all(arange(1._real128, 3._real128) == [1._real128, 2._real128, 3._real128]), &
                  'arange,  real128')
n = n + 1

tests(n) = assert(all(arange(1._real128, 3._real128) == [1._real128, 2._real128, 3._real128]), &
                  'arange,  real128')
n = n + 1

tests(n) = assert(all(arange(cmplx(1._real32, 0._real32), &
                             cmplx(3._real32, 0._real32), &
                             cmplx(1._real32, 0._real32))&
                  == [cmplx(1._real32, 0._real32), &
                      cmplx(2._real32, 0._real32), &
                      cmplx(3._real32, 0._real32)]), &
                  'arange,  complex real32')
n = n + 1

tests(n) = assert(all(arange(cmplx(1._real64, 0._real64), &
                             cmplx(3._real64, 0._real64), &
                             cmplx(1._real64, 0._real64))&
                  == [cmplx(1._real64, 0._real64), &
                      cmplx(2._real64, 0._real64), &
                      cmplx(3._real64, 0._real64)]), &
                  'arange,  complex real64')
n = n + 1

tests(n) = assert(all(arange(cmplx(1._real128, 0._real128), &
                             cmplx(3._real128, 0._real128), &
                             cmplx(1._real128, 0._real128))&
                  == [cmplx(1._real128, 0._real128), &
                      cmplx(2._real128, 0._real128), &
                      cmplx(3._real128, 0._real128)]), &
                  'arange,  complex real128')
n = n + 1

tests(n) = assert(all(arange(cmplx(1, 1), cmplx(3, 5), cmplx(1, 2))&
                  == [cmplx(1, 1), cmplx(2, 3), cmplx(3, 5)]), &
                  'arange,  incrementing both parts of complex numbers')
n = n + 1

tests(n) = assert(all(arange(1, 10) == arange(1, 10, 1)), &
                  'arange increment equals 1 when ommited')
n = n + 1

tests(n) = assert(all(arange(1., 10.) == real(arange(1, 10, 1))), &
                  'integer and real arange variants produce same values')
n = n + 1

tests(n) = assert(all(arange(0., 2.4, 0.8) == [0., 0.8, 1.6, 2.4]), &
                  'custom increment value')
n = n + 1

tests(n) = assert(all(arange(3, -1, -1) == [3, 2, 1, 0, -1]), &
                  'negative increment value')
n = n + 1

tests(n) = assert(size(arange(1.0, 1.4, 0.1)) == 5, &
                  'real32-typed arange returns array of expected size')
n = n + 1

tests(n) = assert(size(arange(1.0_real64, 1.4_real64, 0.1_real64)) == 5, &
                  'real64-typed arange returns array of expected size')
n = n + 1

tests(n) = assert(size(arange(1.0_real128, 1.4_real128, 0.1_real128)) == 5, &
                  'real128-typed arange returns array of expected size')
n = n + 1

test_failed = .false.
call report_tests(tests, test_failed)
if(test_failed)stop 1

end program test_arange

max_stars_repo_name: zmoon/functional-fortran

max_issues_repo_name: zmoon/functional-fortran

max_forks_repo_name: zmoon/functional-fortran


text:
Program Hello
    Print *, "Hello World!"
End Program

max_stars_repo_name: NotteShock/Hello-world-fever

max_issues_repo_name: NotteShock/Hello-world-fever

max_forks_repo_name: NotteShock/Hello-world-fever


text:
mile signs, a Music Scene local musician, is Josiah of the Future on guitar, with vocals.  (Josiah is also in The Musical Drum Key, and in the band earthly).  They have some music on http://milesigns.bandcamp.com/ their website.

Traveling throughout the USA, especially Humboldt, Portland, Seattle, Albuquerque, Santa Fe, Denver, Reno, Las Vegas, Los Angeles, and all over California, Josiah was inspired by the commonality of underground music scenes throughout American cities, as well as by the street musicians he met and busked with along the way.  In July 2009, When Josiah and Time Lane when to do an art show in Santa Ana, CA at Neue Transit Studio, back at the motel, Josiah wrote The Watering Hole, which would become the first in the song series to become the album, and the artist mile signs.

Releases from this time period (Feb. 2010) include :

  mile signs
  warm water and sea salt
  the symmetry of a broken heart demo ep
  a way paved

Josiah wrote the bulk of what would become these records in Davis, CA; specifically around the original Mishkas Cafe,  where he would play almost every day during the bulk of 2009  2010.  He played and wrote songs, in front of, across the street, and out back in the alley.

Mile Signs did a minitour of the bay area to debut the album, doing shows at El Rio, Brainwash Cafe in San Francisco, and The Stork Club in Oakland, which resulted in the live recording and release :

  Live at The Stork Club

Upon moving to Seattle, WA (For a year and a half)

Josiah recorded :

  earthly (demo) (2011)

which would become not only an album, but create a movement whereby Josiah has merged his musical projects (Mile Signs & The Musical Drum Key)  together to form a band which encompasses all of his musical styles.  This band called earthly, which is based in Davis, CA.

ext: f

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search