text:
!WRF:MODEL_LAYER:PHYSICS
!
module module_sf_sfclayrev

    real, parameter :: vconvc = 1.
    real, parameter :: czo = 0.0185
    real, parameter :: ozo = 1.59e-5

    real, dimension(0:1000), save :: psim_stab, psim_unstab, psih_stab, psih_unstab

contains

    !-------------------------------------------------------------------
    subroutine sfclayrev(u3d, v3d, t3d, qv3d, p3d, dz8w, &
            cp, g, rovcp, r, xlv, psfc, chs, chs2, cqs2, cpm, &
            znt, ust, pblh, mavail, zol, mol, regime, psim, psih, &
            fm, fh, &
            xland, hfx, qfx, lh, tsk, flhc, flqc, qgh, qsfc, rmol, &
            u10, v10, th2, t2, q2, &
            gz1oz0, wspd, br, isfflx, dx, &
            svp1, svp2, svp3, svpt0, ep1, ep2, &
            karman, eomeg, stbolt, &
            p1000mb, &
            ids, ide, jds, jde, kds, kde, &
            ims, ime, jms, jme, kms, kme, &
            its, ite, jts, jte, kts, kte, &
            ustm, ck, cka, cd, cda, isftcflx, iz0tlnd, scm_force_flux)
        !-------------------------------------------------------------------
        implicit none
        !-------------------------------------------------------------------
        !   changes in v3.7 over water surfaces:
        !          1. for znt/cd, replacing constant ozo with 0.11*1.5e-5/ust(i)
        !             the coare 3.5 (edson et al. 2013) formulation is also available
        !          2. for vconv, reducing magnitude by half
        !          3. for ck, replacing carlson-boland with coare 3
        !-------------------------------------------------------------------
        !-- u3d         3d u-velocity interpolated to theta points (m/s)
        !-- v3d         3d v-velocity interpolated to theta points (m/s)
        !-- t3d         temperature (k)
        !-- qv3d        3d water vapor mixing ratio (kg/kg)
        !-- p3d         3d pressure (pa)
        !-- dz8w        dz between full levels (m)
        !-- cp          heat capacity at constant pressure for dry air (j/kg/k)
        !-- g           acceleration due to gravity (m/s^2)
        !-- rovcp       r/cp
        !-- r           gas constant for dry air (j/kg/k)
        !-- xlv         latent heat of vaporization for water (j/kg)
        !-- psfc        surface pressure (pa)
        !-- znt         roughness length (m)
        !-- ust         u* in similarity theory (m/s)
        !-- ustm        u* in similarity theory (m/s) without vconv correction
        !               used to couple with tke scheme
        !-- pblh        pbl height from previous time (m)
        !-- mavail      surface moisture availability (between 0 and 1)
        !-- zol         z/l height over monin-obukhov length
        !-- mol         t* (similarity theory) (k)
        !-- regime      flag indicating pbl regime (stable, unstable, etc.)
        !-- psim        similarity stability function for momentum
        !-- psih        similarity stability function for heat
        !-- fm          integrated stability function for momentum
        !-- fh          integrated stability function for heat
        !-- xland       land mask (1 for land, 2 for water)
        !-- hfx         upward heat flux at the surface (w/m^2)
        !-- qfx         upward moisture flux at the surface (kg/m^2/s)
        !-- lh          net upward latent heat flux at surface (w/m^2)
        !-- tsk         surface temperature (k)
        !-- flhc        exchange coefficient for heat (w/m^2/k)
        !-- flqc        exchange coefficient for moisture (kg/m^2/s)
        !-- chs         heat/moisture exchange coefficient for lsm (m/s)
        !-- qgh         lowest-level saturated mixing ratio
        !-- qsfc        ground saturated mixing ratio
        !-- u10         diagnostic 10m u wind
        !-- v10         diagnostic 10m v wind
        !-- th2         diagnostic 2m theta (k)
        !-- t2          diagnostic 2m temperature (k)
        !-- q2          diagnostic 2m mixing ratio (kg/kg)
        !-- gz1oz0      log(z/z0) where z0 is roughness length
        !-- wspd        wind speed at lowest model level (m/s)
        !-- br          bulk richardson number in surface layer
        !-- isfflx      isfflx=1 for surface heat and moisture fluxes
        !-- dx          horizontal grid size (m)
        !-- svp1        constant for saturation vapor pressure (kpa)
        !-- svp2        constant for saturation vapor pressure (dimensionless)
        !-- svp3        constant for saturation vapor pressure (k)
        !-- svpt0       constant for saturation vapor pressure (k)
        !-- ep1         constant for virtual temperature (r_v/r_d - 1) (dimensionless)
        !-- ep2         constant for specific humidity calculation
        !               (r_d/r_v) (dimensionless)
        !-- karman      von karman constant
        !-- eomeg       angular velocity of earth's rotation (rad/s)
        !-- stbolt      stefan-boltzmann constant (w/m^2/k^4)
        !-- ck          enthalpy exchange coeff at 10 meters
        !-- cd          momentum exchange coeff at 10 meters
        !-- cka         enthalpy exchange coeff at the lowest model level
        !-- cda         momentum exchange coeff at the lowest model level
        !-- isftcflx    =0, (charnock and carlson-boland); =1, ahw ck, cd, =2 garratt
        !-- iz0tlnd     =0 carlson-boland, =1 czil_new
        !-- ids         start index for i in domain
        !-- ide         end index for i in domain
        !-- jds         start index for j in domain
        !-- jde         end index for j in domain
        !-- kds         start index for k in domain
        !-- kde         end index for k in domain
        !-- ims         start index for i in memory
        !-- ime         end index for i in memory
        !-- jms         start index for j in memory
        !-- jme         end index for j in memory
        !-- kms         start index for k in memory
        !-- kme         end index for k in memory
        !-- its         start index for i in tile
        !-- ite         end index for i in tile
        !-- jts         start index for j in tile
        !-- jte         end index for j in tile
        !-- kts         start index for k in tile
        !-- kte         end index for k in tile
        !-------------------------------------------------------------------
        integer, intent(in) :: ids, ide, jds, jde, kds, kde, &
                ims, ime, jms, jme, kms, kme, &
                its, ite, jts, jte, kts, kte
        !
        integer, intent(in) :: isfflx
        real, intent(in) :: svp1, svp2, svp3, svpt0
        real, intent(in) :: ep1, ep2, karman, eomeg, stbolt
        real, intent(in) :: p1000mb
        !
        real, dimension(ims:ime, kms:kme, jms:jme), &
                intent(in) :: dz8w

        real, dimension(ims:ime, kms:kme, jms:jme), &
                intent(in) :: qv3d, &
                p3d, &
                t3d

        real, dimension(ims:ime, jms:jme), &
                intent(in) :: mavail, &
                pblh, &
                xland, &
                tsk
        real, dimension(ims:ime, jms:jme), &
                intent(out) :: u10, &
                v10, &
                th2, &
                t2, &
                q2, &
                qsfc

        !
        real, dimension(ims:ime, jms:jme), &
                intent(inout) :: regime, &
                hfx, &
                qfx, &
                lh, &
                mol, rmol
        !m the following 5 are change to memory size
        !
        real, dimension(ims:ime, jms:jme), &
                intent(inout) :: gz1oz0, wspd, br, &
                psim, psih, fm, fh

        real, dimension(ims:ime, kms:kme, jms:jme), &
                intent(in) :: u3d, &
                v3d

        real, dimension(ims:ime, jms:jme), &
                intent(in) :: psfc

        real, dimension(ims:ime, jms:jme), &
                intent(inout) :: znt, &
                zol, &
                ust, &
                cpm, &
                chs2, &
                cqs2, &
                chs

        real, dimension(ims:ime, jms:jme), &
                intent(inout) :: flhc, flqc

        real, dimension(ims:ime, jms:jme), &
                intent(inout) :: &
                qgh

        real, intent(in) :: cp, g, rovcp, r, xlv, dx

        real, optional, dimension(ims:ime, jms:jme), &
                intent(out) :: ck, cka, cd, cda

        real, optional, dimension(ims:ime, jms:jme), &
                intent(inout) :: ustm

        integer, optional, intent(in) :: isftcflx, iz0tlnd
        integer, optional, intent(in) :: scm_force_flux
        ! local vars

        real, dimension(its:ite) :: u1d, &
                v1d, &
                qv1d, &
                p1d, &
                t1d

        real, dimension(its:ite) :: dz8w1d

        integer :: i, j

        do j = jts, jte
            do i = its, ite
                dz8w1d(i) = dz8w(i, 1, j)
            enddo

            do i = its, ite
                u1d(i) = u3d(i, 1, j)
                v1d(i) = v3d(i, 1, j)
                qv1d(i) = qv3d(i, 1, j)
                p1d(i) = p3d(i, 1, j)
                t1d(i) = t3d(i, 1, j)
            enddo

            !  sending array starting locations of optional variables may cause
            !  troubles, so we explicitly change the call.

            call sfclayrev1d(j, u1d, v1d, t1d, qv1d, p1d, dz8w1d, &
                    cp, g, rovcp, r, xlv, psfc(ims, j), chs(ims, j), chs2(ims, j), &
                    cqs2(ims, j), cpm(ims, j), pblh(ims, j), rmol(ims, j), &
                    znt(ims, j), ust(ims, j), mavail(ims, j), zol(ims, j), &
                    mol(ims, j), regime(ims, j), psim(ims, j), psih(ims, j), &
                    fm(ims, j), fh(ims, j), &
                    xland(ims, j), hfx(ims, j), qfx(ims, j), tsk(ims, j), &
                    u10(ims, j), v10(ims, j), th2(ims, j), t2(ims, j), &
                    q2(ims, j), flhc(ims, j), flqc(ims, j), qgh(ims, j), &
                    qsfc(ims, j), lh(ims, j), &
                    gz1oz0(ims, j), wspd(ims, j), br(ims, j), isfflx, dx, &
                    svp1, svp2, svp3, svpt0, ep1, ep2, karman, eomeg, stbolt, &
                    p1000mb, &
                    ids, ide, jds, jde, kds, kde, &
                    ims, ime, jms, jme, kms, kme, &
                    its, ite, jts, jte, kts, kte                          &
                    #if ( em_core == 1 )
, isftcflx, iz0tlnd, scm_force_flux, &
                    ustm(ims, j), ck(ims, j), cka(ims, j), &
                    cd(ims, j), cda(ims, j)                               &
                    #endif
)
        enddo

    end subroutine sfclayrev

    !-------------------------------------------------------------------
    subroutine sfclayrev1d(j, ux, vx, t1d, qv1d, p1d, dz8w1d, &
            cp, g, rovcp, r, xlv, psfcpa, chs, chs2, cqs2, cpm, pblh, rmol, &
            znt, ust, mavail, zol, mol, regime, psim, psih, fm, fh, &
            xland, hfx, qfx, tsk, &
            u10, v10, th2, t2, q2, flhc, flqc, qgh, &
            qsfc, lh, gz1oz0, wspd, br, isfflx, dx, &
            svp1, svp2, svp3, svpt0, ep1, ep2, &
            karman, eomeg, stbolt, &
            p1000mb, &
            ids, ide, jds, jde, kds, kde, &
            ims, ime, jms, jme, kms, kme, &
            its, ite, jts, jte, kts, kte, &
            isftcflx, iz0tlnd, scm_force_flux, &
            ustm, ck, cka, cd, cda)
        !-------------------------------------------------------------------
        implicit none
        !-------------------------------------------------------------------
        !-------------------------------------------------------------------
        !   changes in v3.7 over water surfaces:
        !          1. for znt/cd, replacing constant ozo with 0.11*1.5e-5/ust(i)
        !             the coare 3.5 (edson et al. 2013) formulation is also available
        !          2. for vconv, reducing magnitude by half
        !          3. for ck, replacing carlson-boland with coare 3
        !-------------------------------------------------------------------
        !-- ux          1d u-velocity interpolated to theta points (m/s)
        !-- vx         1d v-velocity interpolated to theta points (m/s)
        !-- t1d         temperature (k)
        !-- qv1d        3d water vapor mixing ratio (kg/kg)
        !-- p1d         3d pressure (pa)
        !-- dz8w1d        dz between full levels (m)
        !-- cp          heat capacity at constant pressure for dry air (j/kg/k)
        !-- g           acceleration due to gravity (m/s^2)
        !-- rovcp       r/cp
        !-- r           gas constant for dry air (j/kg/k)
        !-- xlv         latent heat of vaporization for water (j/kg)
        !-- psfcpa        surface pressure (pa)
        !-- chs         heat/moisture exchange coefficient for lsm (m/s)
        !-- chs2        ???
        !-- cqs2          ???
        !-- cpm         heat capacity at a constant pressure for moist air cpm(i) = cp * (1. + 0.8 * qx(i))
        !-- pblh        pbl height from previous time (m)
        !-- rmol        1 / monin obukhov length
        !-- znt         z naught, zo, roughness length (m)
        !-- ust         u* in similarity theory (m/s)
        !-- ustm        u* in similarity theory (m/s) without vconv correction
        !               used to couple with tke scheme
        !-- mavail      surface moisture availability (between 0 and 1)
        !-- zol         z/l height over monin-obukhov length
        !-- mol         t* (similarity theory) (k) (Actually this might be monin-obukhov length but am not sure)
        !-- regime      flag indicating pbl regime (stable, unstable, etc.)
        !-- psim        similarity stability function for momentum
        !-- psih        similarity stability function for heat
        !-- fm          integrated stability function for momentum
        !-- fh          integrated stability function for heat
        !-- xland       land mask (1 for land, 2 for water)
        !-- hfx         upward heat flux at the surface (w/m^2)
        !-- qfx         upward moisture flux at the surface (kg/m^2/s)
        !-- lh          net upward latent heat flux at surface (w/m^2)
        !-- tsk         surface temperature (k)
        !-- flhc        exchange coefficient for heat (w/m^2/k)
        !-- flqc        exchange coefficient for moisture (kg/m^2/s)
        !-- qgh         lowest-level saturated mixing ratio
        !-- qsfc        ground saturated mixing ratio
        !-- u10         diagnostic 10m u wind
        !-- v10         diagnostic 10m v wind
        !-- th2         diagnostic 2m theta (k)
        !-- t2          diagnostic 2m temperature (k)
        !-- q2          diagnostic 2m mixing ratio (kg/kg)
        !-- gz1oz0      log(z/z0) where z0 is roughness length
        !-- wspd        wind speed at lowest model level (m/s)
        !-- br          bulk richardson number in surface layer
        !-- isfflx      isfflx=1 for surface heat and moisture fluxes
        !-- dx          horizontal grid size (m)
        !-- svp1        constant for saturation vapor pressure (kpa)
        !-- svp2        constant for saturation vapor pressure (dimensionless)
        !-- svp3        constant for saturation vapor pressure (k)
        !-- svpt0       constant for saturation vapor pressure (k)
        !-- ep1         constant for virtual temperature (r_v/r_d - 1) (dimensionless)
        !-- ep2         constant for specific humidity calculation
        !               (r_d/r_v) (dimensionless)
        !-- karman      von karman constant
        !-- eomeg       angular velocity of earth's rotation (rad/s)
        !-- stbolt      stefan-boltzmann constant (w/m^2/k^4)
        !-- ck          enthalpy exchange coeff at 10 meters
        !-- cd          momentum exchange coeff at 10 meters
        !-- cka         enthalpy exchange coeff at the lowest model level
        !-- cda         momentum exchange coeff at the lowest model level
        !-- isftcflx    =0, (charnock and carlson-boland); =1, ahw ck, cd, =2 garratt
        !-- iz0tlnd     =0 carlson-boland, =1 czil_new
        !-- ids         start index for i in domain
        !-- ide         end index for i in domain
        !-- jds         start index for j in domain
        !-- jde         end index for j in domain
        !-- kds         start index for k in domain
        !-- kde         end index for k in domain
        !-- ims         start index for i in memory
        !-- ime         end index for i in memory
        !-- jms         start index for j in memory
        !-- jme         end index for j in memory
        !-- kms         start index for k in memory
        !-- kme         end index for k in memory
        !-- its         start index for i in tile
        !-- ite         end index for i in tile
        !-- jts         start index for j in tile
        !-- jte         end index for j in tile
        !-- kts         start index for k in tile
        !-- kte         end index for k in tile
        !-------------------------------------------------------------------
        real, parameter :: xka = 2.4e-5
        real, parameter :: prt = 1.

        integer, intent(in) :: ids, ide, jds, jde, kds, kde, &
                ims, ime, jms, jme, kms, kme, &
                its, ite, jts, jte, kts, kte, &
                j
        !
        integer, intent(in) :: isfflx
        real, intent(in) :: svp1, svp2, svp3, svpt0
        real, intent(in) :: ep1, ep2, karman, eomeg, stbolt
        real, intent(in) :: p1000mb

        !
        real, dimension(ims:ime), &
                intent(in) :: mavail, &
                pblh, &
                xland, &
                tsk
        !
        real, dimension(ims:ime), &
                intent(in) :: psfcpa

        real, dimension(ims:ime), &
                intent(inout) :: regime, &
                hfx, &
                qfx, &
                mol, rmol
        !m the following 5 are changed to memory size---
        !
        real, dimension(ims:ime), &
                intent(inout) :: gz1oz0, wspd, br, &
                psim, psih, fm, fh

        real, dimension(ims:ime), &
                intent(inout) :: znt, &
                zol, &
                ust, &
                cpm, &
                chs2, &
                cqs2, &
                chs

        real, dimension(ims:ime), &
                intent(inout) :: flhc, flqc

        real, dimension(ims:ime), &
                intent(inout) :: &
                qgh

        real, dimension(ims:ime), &
                intent(out) :: u10, v10, &
                th2, t2, q2, qsfc, lh

        real, intent(in) :: cp, g, rovcp, r, xlv, dx

        ! module-local variables, defined in subroutine sfclay
        real, dimension(its:ite), intent(in) :: dz8w1d

        real, dimension(its:ite), intent(in) :: ux, &
                vx, &
                qv1d, &
                p1d, &
                t1d

        real, optional, dimension(ims:ime), &
                intent(out) :: ck, cka, cd, cda
        real, optional, dimension(ims:ime), &
                intent(inout) :: ustm

        integer, optional, intent(in) :: isftcflx, iz0tlnd
        integer, optional, intent(in) :: scm_force_flux

        ! local vars

        real, dimension(its:ite) :: za, &
                thvx, zqkl, &
                zqklp1, &
                thx, qx, &
                psih2, &
                psim2, &
                psih10, &
                psim10, &
                denomq, &
                denomq2, &
                denomt2, &
                wspdi, &
                gz2oz0, &
                gz10oz0
        !
        real, dimension(its:ite) :: &
                rhox, govrth, &
                tgdsa
        !
        real, dimension(its:ite) :: scr3, scr4
        real, dimension(its:ite) :: thgb, psfc
        !
        integer :: kl

        integer :: n, i, k, kk, l, nzol, nk, nzol2, nzol10

        real :: pl, thcon, tvcon, e1
        real :: zl, tskv, dthvdz, dthvm, vconv, rzol, rzol2, rzol10, zol2, zol10
        real :: dtg, psix, dtthx, psix10, psit, psit2, psiq, psiq2, psiq10
        real :: fluxc, vsgd, z0q, visc, restar, czil, gz0ozq, gz0ozt
        real :: zw, zn1, zn2
        !
        ! .... paj ...
        !
        real :: zolzz, zol0
        !     real    :: zolri,zolri2
        !     real    :: psih_stable,psim_stable,psih_unstable,psim_unstable
        !     real    :: psih_stable_full,psim_stable_full,psih_unstable_full,psim_unstable_full
        real :: zl2, zl10, z0t
        real, dimension(its:ite) :: pq, pq2, pq10
        ! thgb      potential skin temperature (K)
        ! tskv
        !-------------------------------------------------------------------
        kl = kte

        do i = its, ite
            ! psfc cb
            psfc(i) = psfcpa(i) / 1000.
        enddo
        !
        !----convert ground temperature to potential temperature:
        !
        do 5 i = its, ite
            tgdsa(i) = tsk(i)
            ! psfc cb
            !        thgb(i)=tsk(i)*(100./psfc(i))**rovcp
            thgb(i) = tsk(i) * (p1000mb / psfcpa(i))**rovcp
        5 continue
        !
        !-----decouple flux-form variables to give u,v,t,theta,theta-vir.,
        !     t-vir., qv, and qc at cross points and at ktau-1.
        !
        !     *** note ***
        !         the boundary winds may not be adequately affected by friction,
        !         so use only interior values of ux and vx to calculate
        !         tendencies.
        !
        10 continue

        !     do 24 i=its,ite
        !        ux(i)=u1d(i)
        !        vx(i)=v1d(i)
        !  24 continue

        26 continue

        !.....scr3(i,k) store temperature,
        !     scr4(i,k) store virtual temperature.

        do 30 i = its, ite
            ! pl cb
            pl = p1d(i) / 1000.
            scr3(i) = t1d(i)
            !         thcon=(100./pl)**rovcp
            thcon = (p1000mb * 0.001 / pl)**rovcp
            thx(i) = scr3(i) * thcon
            scr4(i) = scr3(i)
            thvx(i) = thx(i)
            qx(i) = 0.
        30 continue
        !
        do i = its, ite
            qgh(i) = 0.
            flhc(i) = 0.
            flqc(i) = 0.
            cpm(i) = cp
        enddo
        !
        !     if(idry.eq.1)goto 80
        do 50 i = its, ite
            qx(i) = qv1d(i)
            tvcon = (1. + ep1 * qx(i))
            thvx(i) = thx(i) * tvcon
            scr4(i) = scr3(i) * tvcon
        50 continue
        !
        do 60 i = its, ite
            e1 = svp1 * exp(svp2 * (tgdsa(i) - svpt0) / (tgdsa(i) - svp3))
            !  for land points qsfc can come from previous time step
            if(xland(i).gt.1.5.or.qsfc(i).le.0.0)qsfc(i) = ep2 * e1 / (psfc(i) - e1)
            ! qgh changed to use lowest-level air temp consistent with myjsfc change
            ! q2sat = qgh in lsm
            e1 = svp1 * exp(svp2 * (t1d(i) - svpt0) / (t1d(i) - svp3))
            pl = p1d(i) / 1000.
            qgh(i) = ep2 * e1 / (pl - e1)
            cpm(i) = cp * (1. + 0.8 * qx(i))
        60 continue
        80 continue

        !-----compute the height of full- and half-sigma levels above ground
        !     level, and the layer thicknesses.

        do 90 i = its, ite
            zqklp1(i) = 0.
            rhox(i) = psfc(i) * 1000. / (r * scr4(i))
        90 continue
        !
        do 110 i = its, ite
            zqkl(i) = dz8w1d(i) + zqklp1(i)
        110 continue
        !
        do 120 i = its, ite
            za(i) = 0.5 * (zqkl(i) + zqklp1(i))
        120 continue
        !
        do 160 i = its, ite
            govrth(i) = g / thx(i)
        160 continue

        !-----calculate bulk richardson no. of surface layer, according to
        !     akb(1976), eq(12).

        do 260 i = its, ite
            gz1oz0(i) = alog((za(i) + znt(i)) / znt(i))   ! log((z+z0)/z0)
            gz2oz0(i) = alog((2. + znt(i)) / znt(i))      ! log((2+z0)/z0)
            gz10oz0(i) = alog((10. + znt(i)) / znt(i))    ! log((10+z0)z0)
            if((xland(i) - 1.5).ge.0)then
                zl = znt(i)
            else
                zl = 0.01
            endif
            wspd(i) = sqrt(ux(i) * ux(i) + vx(i) * vx(i))

            tskv = thgb(i) * (1. + ep1 * qsfc(i))
            dthvdz = (thvx(i) - tskv)
            !  convective velocity scale vc and subgrid-scale velocity vsg
            !  following beljaars (1994, qjrms) and mahrt and sun (1995, mwr)
            !                                ... hong aug. 2001
            !
            !       vconv = 0.25*sqrt(g/tskv*pblh(i)*dthvm)
            !      use beljaars over land, old mm5 (wyngaard) formula over water
            if (xland(i).lt.1.5) then
                fluxc = max(hfx(i) / rhox(i) / cp                    &
                        + ep1 * tskv * qfx(i) / rhox(i), 0.)
                vconv = vconvc * (g / tgdsa(i) * pblh(i) * fluxc)**.33
            else
                if(-dthvdz.ge.0)then
                    dthvm = -dthvdz
                else
                    dthvm = 0.
                endif
                !       vconv = 2.*sqrt(dthvm)
                ! v3.7: reducing contribution in calm conditions
                vconv = sqrt(dthvm)
            endif
            ! mahrt and sun low-res correction
            vsgd = 0.32 * (max(dx / 5000. - 1., 0.))**.33
            wspd(i) = sqrt(wspd(i) * wspd(i) + vconv * vconv + vsgd * vsgd)
            wspd(i) = amax1(wspd(i), 0.1)
            br(i) = govrth(i) * za(i) * dthvdz / (wspd(i) * wspd(i))
            !  if previously unstable, do not let into regimes 1 and 2
            if(mol(i).lt.0.)br(i) = amin1(br(i), 0.0)
            !jdf
            rmol(i) = -govrth(i) * dthvdz * za(i) * karman
            !jdf

        260 continue

        !
        !-----diagnose basic parameters for the appropriated stability class:
        !
        !
        !     the stability classes are determined by br (bulk richardson no.)
        !     and hol (height of pbl/monin-obukhov length).
        !
        !     criteria for the classes are as follows:
        !
        !        1. br .ge. 0.0;
        !               represents nighttime stable conditions (regime=1),
        !
        !        3. br .eq. 0.0
        !               represents forced convection conditions (regime=3),
        !
        !        4. br .lt. 0.0
        !               represents free convection conditions (regime=4).
        !
        !ccccc

        do 320 i = its, ite
            !
            if (br(i).gt.0) then
                if (br(i).gt.250.0) then
                    zol(i) = zolri(250.0, za(i), znt(i))
                else
                    zol(i) = zolri(br(i), za(i), znt(i))
                endif
            endif
            !
            if (br(i).lt.0) then
                if(ust(i).lt.0.001)then
                    zol(i) = br(i) * gz1oz0(i)
                else
                    if (br(i).lt.-250.0) then
                        zol(i) = zolri(-250.0, za(i), znt(i))
                    else
                        zol(i) = zolri(br(i), za(i), znt(i))
                    endif
                endif
            endif
            !
            ! ... paj: compute integrated similarity functions.
            !
            zolzz = zol(i) * (za(i) + znt(i)) / za(i) ! (z+z0/l
            zol10 = zol(i) * (10. + znt(i)) / za(i)   ! (10+z0)/l
            zol2 = zol(i) * (2. + znt(i)) / za(i)     ! (2+z0)/l
            zol0 = zol(i) * znt(i) / za(i)          ! z0/l
            zl2 = (2.) / za(i) * zol(i)             ! 2/l
            zl10 = (10.) / za(i) * zol(i)           ! 10/l

            if((xland(i) - 1.5).lt.0.)then
                zl = (0.01) / za(i) * zol(i)   ! (0.01)/l
            else
                zl = zol0                     ! z0/l
            endif

            if(br(i).lt.0.)goto 310  ! go to unstable regime (class 4)
            if(br(i).eq.0.)goto 280  ! go to neutral regime (class 3)
            !
            !-----class 1; stable (nighttime) conditions:
            !
            regime(i) = 1.
            !
            ! ... paj: psim and psih. follows cheng and brutsaert 2005 (cb05).
            !
            psim(i) = psim_stable(zolzz) - psim_stable(zol0)
            psih(i) = psih_stable(zolzz) - psih_stable(zol0)
            !
            psim10(i) = psim_stable(zol10) - psim_stable(zol0)
            psih10(i) = psih_stable(zol10) - psih_stable(zol0)
            !
            psim2(i) = psim_stable(zol2) - psim_stable(zol0)
            psih2(i) = psih_stable(zol2) - psih_stable(zol0)
            !
            ! ... paj: preparations to compute psiq. follows cb05+carlson boland jam 1978.
            !
            pq(i) = psih_stable(zol(i)) - psih_stable(zl)
            pq2(i) = psih_stable(zl2) - psih_stable(zl)
            pq10(i) = psih_stable(zl10) - psih_stable(zl)
            !
            !       1.0 over monin-obukhov length
            rmol(i) = zol(i) / za(i)
            !

            goto 320
            !
            !-----class 3; forced convection:
            !
            280   regime(i) = 3.
            psim(i) = 0.0
            psih(i) = psim(i)
            psim10(i) = 0.
            psih10(i) = psim10(i)
            psim2(i) = 0.
            psih2(i) = psim2(i)
            !
            ! paj: preparations to compute psiq.
            !
            pq(i) = psih(i)
            pq2(i) = psih2(i)
            pq10(i) = 0.
            !
            zol(i) = 0.
            rmol(i) = zol(i) / za(i)

            goto 320
            !
            !-----class 4; free convection:
            !
            310   continue
            regime(i) = 4.
            !
            ! ... paj: psim and psih ...
            !
            psim(i) = psim_unstable(zolzz) - psim_unstable(zol0)
            psih(i) = psih_unstable(zolzz) - psih_unstable(zol0)
            !
            psim10(i) = psim_unstable(zol10) - psim_unstable(zol0)
            psih10(i) = psih_unstable(zol10) - psih_unstable(zol0)
            !
            psim2(i) = psim_unstable(zol2) - psim_unstable(zol0)
            psih2(i) = psih_unstable(zol2) - psih_unstable(zol0)
            !
            ! ... paj: preparations to compute psiq
            !
            pq(i) = psih_unstable(zol(i)) - psih_unstable(zl)
            pq2(i) = psih_unstable(zl2) - psih_unstable(zl)
            pq10(i) = psih_unstable(zl10) - psih_unstable(zl)
            !
            !---limiot psih and psim in the case of thin layers and high roughness
            !---  this prevents denominator in fluxes from getting too small
            psih(i) = amin1(psih(i), 0.9 * gz1oz0(i))
            psim(i) = amin1(psim(i), 0.9 * gz1oz0(i))
            psih2(i) = amin1(psih2(i), 0.9 * gz2oz0(i))
            psim10(i) = amin1(psim10(i), 0.9 * gz10oz0(i))
            !
            ! ahw: mods to compute ck, cd
            psih10(i) = amin1(psih10(i), 0.9 * gz10oz0(i))

            rmol(i) = zol(i) / za(i)

        320 continue
        !
        !-----compute the frictional velocity:
        !     za(1982) eqs(2.60),(2.61).
        !
        do 330 i = its, ite
            dtg = thx(i) - thgb(i)
            psix = gz1oz0(i) - psim(i)
            psix10 = gz10oz0(i) - psim10(i)

            !     lower limit added to prevent large flhc in soil model
            !     activates in unstable conditions with thin layers or high z0
            !       psit=amax1(gz1oz0(i)-psih(i),2.)
            psit = gz1oz0(i) - psih(i)
            psit2 = gz2oz0(i) - psih2(i)
            !
            if((xland(i) - 1.5).ge.0)then
                zl = znt(i)
            else
                zl = 0.01
            endif
            !
            psiq = alog(karman * ust(i) * za(i) / xka + za(i) / zl) - pq(i)
            psiq2 = alog(karman * ust(i) * 2. / xka + 2. / zl) - pq2(i)

            ! ahw: mods to compute ck, cd
            psiq10 = alog(karman * ust(i) * 10. / xka + 10. / zl) - pq10(i)

            ! v3.7: using fairall 2003 to compute z0q and z0t over water:
            !       adapted from module_sf_mynn.f
            if ((xland(i) - 1.5).ge.0.) then
                visc = (1.32 + 0.009 * (scr3(i) - 273.15)) * 1.e-5
                restar = ust(i) * znt(i) / visc
                z0t = (5.5e-5) * (restar**(-0.60))
                z0t = min(z0t, 1.0e-4)
                z0t = max(z0t, 2.0e-9)
                z0q = z0t

                psiq = max(alog((za(i) + z0q) / z0q) - psih(i), 2.)
                psit = max(alog((za(i) + z0t) / z0t) - psih(i), 2.)
                psiq2 = max(alog((2. + z0q) / z0q) - psih2(i), 2.)
                psit2 = max(alog((2. + z0t) / z0t) - psih2(i), 2.)
                psiq10 = max(alog((10. + z0q) / z0q) - psih10(i), 2.)
            endif

            if (present(isftcflx)) then
                if (isftcflx.eq.1 .and. (xland(i) - 1.5).ge.0.) then
                    ! v3.1
                    !             z0q = 1.e-4 + 1.e-3*(max(0.,ust(i)-1.))**2
                    ! hfip1
                    !             z0q = 0.62*2.0e-5/ust(i) + 1.e-3*(max(0.,ust(i)-1.5))**2
                    ! v3.2
                    z0q = 1.e-4
                    !
                    ! ... paj: recompute psih for z0q
                    !
                    zolzz = zol(i) * (za(i) + z0q) / za(i)    ! (z+z0q)/l
                    zol10 = zol(i) * (10. + z0q) / za(i)   ! (10+z0q)/l
                    zol2 = zol(i) * (2. + z0q) / za(i)     ! (2+z0q)/l
                    zol0 = zol(i) * z0q / za(i)          ! z0q/l
                    !
                    if (zol(i).gt.0.) then
                        psih(i) = psih_stable(zolzz) - psih_stable(zol0)
                        psih10(i) = psih_stable(zol10) - psih_stable(zol0)
                        psih2(i) = psih_stable(zol2) - psih_stable(zol0)
                    else
                        if (zol(i).eq.0) then
                            psih(i) = 0.
                            psih10(i) = 0.
                            psih2(i) = 0.
                        else
                            psih(i) = psih_unstable(zolzz) - psih_unstable(zol0)
                            psih10(i) = psih_unstable(zol10) - psih_unstable(zol0)
                            psih2(i) = psih_unstable(zol2) - psih_unstable(zol0)
                        endif
                    endif
                    !
                    psiq = alog((za(i) + z0q) / z0q) - psih(i)
                    psit = psiq
                    psiq2 = alog((2. + z0q) / z0q) - psih2(i)
                    psiq10 = alog((10. + z0q) / z0q) - psih10(i)
                    psit2 = psiq2
                endif
                if (isftcflx.eq.2 .and. (xland(i) - 1.5).ge.0.) then
                    ! ahw: garratt formula: calculate roughness reynolds number
                    !        kinematic viscosity of air (linear approc to
                    !                 temp dependence at sea level)
                    ! gz0ozt and gz0ozq are based off formulas from brutsaert (1975), which
                    ! garratt (1992) used with values of k = 0.40, pr = 0.71, and sc = 0.60
                    visc = (1.32 + 0.009 * (scr3(i) - 273.15)) * 1.e-5
                    !!            visc=1.5e-5
                    restar = ust(i) * znt(i) / visc
                    gz0ozt = 0.40 * (7.3 * sqrt(sqrt(restar)) * sqrt(0.71) - 5.)
                    !
                    ! ... paj: compute psih for z0t for temperature ...
                    !
                    z0t = znt(i) / exp(gz0ozt)
                    !
                    zolzz = zol(i) * (za(i) + z0t) / za(i)    ! (z+z0t)/l
                    zol10 = zol(i) * (10. + z0t) / za(i)   ! (10+z0t)/l
                    zol2 = zol(i) * (2. + z0t) / za(i)     ! (2+z0t)/l
                    zol0 = zol(i) * z0t / za(i)          ! z0t/l
                    !
                    if (zol(i).gt.0.) then
                        psih(i) = psih_stable(zolzz) - psih_stable(zol0)
                        psih10(i) = psih_stable(zol10) - psih_stable(zol0)
                        psih2(i) = psih_stable(zol2) - psih_stable(zol0)
                    else
                        if (zol(i).eq.0) then
                            psih(i) = 0.
                            psih10(i) = 0.
                            psih2(i) = 0.
                        else
                            psih(i) = psih_unstable(zolzz) - psih_unstable(zol0)
                            psih10(i) = psih_unstable(zol10) - psih_unstable(zol0)
                            psih2(i) = psih_unstable(zol2) - psih_unstable(zol0)
                        endif
                    endif
                    !
                    !              psit=gz1oz0(i)-psih(i)+restar2
                    !              psit2=gz2oz0(i)-psih2(i)+restar2
                    psit = alog((za(i) + z0t) / z0t) - psih(i)
                    psit2 = alog((2. + z0t) / z0t) - psih2(i)
                    !
                    gz0ozq = 0.40 * (7.3 * sqrt(sqrt(restar)) * sqrt(0.60) - 5.)
                    z0q = znt(i) / exp(gz0ozq)
                    !
                    zolzz = zol(i) * (za(i) + z0q) / za(i)    ! (z+z0q)/l
                    zol10 = zol(i) * (10. + z0q) / za(i)   ! (10+z0q)/l
                    zol2 = zol(i) * (2. + z0q) / za(i)     ! (2+z0q)/l
                    zol0 = zol(i) * z0q / za(i)          ! z0q/l
                    !
                    if (zol(i).gt.0.) then
                        psih(i) = psih_stable(zolzz) - psih_stable(zol0)
                        psih10(i) = psih_stable(zol10) - psih_stable(zol0)
                        psih2(i) = psih_stable(zol2) - psih_stable(zol0)
                    else
                        if (zol(i).eq.0) then
                            psih(i) = 0.
                            psih10(i) = 0.
                            psih2(i) = 0.
                        else
                            psih(i) = psih_unstable(zolzz) - psih_unstable(zol0)
                            psih10(i) = psih_unstable(zol10) - psih_unstable(zol0)
                            psih2(i) = psih_unstable(zol2) - psih_unstable(zol0)
                        endif
                    endif
                    !
                    psiq = alog((za(i) + z0q) / z0q) - psih(i)
                    psiq2 = alog((2. + z0q) / z0q) - psih2(i)
                    psiq10 = alog((10. + z0q) / z0q) - psih10(i)
                    !              psiq=gz1oz0(i)-psih(i)+2.28*sqrt(sqrt(restar))-2.
                    !              psiq2=gz2oz0(i)-psih2(i)+2.28*sqrt(sqrt(restar))-2.
                    !              psiq10=gz10oz0(i)-psih(i)+2.28*sqrt(sqrt(restar))-2.
                endif
            endif
            if(present(ck) .and. present(cd) .and. present(cka) .and. present(cda)) then
                ck(i) = (karman / psix10) * (karman / psiq10)
                cd(i) = (karman / psix10) * (karman / psix10)
                cka(i) = (karman / psix) * (karman / psiq)
                cda(i) = (karman / psix) * (karman / psix)
            endif
            if (present(iz0tlnd)) then
                if (iz0tlnd.eq.1 .and. (xland(i) - 1.5).le.0.) then
                    zl = znt(i)
                    !             czil related changes for land
                    visc = (1.32 + 0.009 * (scr3(i) - 273.15)) * 1.e-5
                    restar = ust(i) * zl / visc
                    !             modify czil according to chen & zhang, 2009

                    czil = 10.0 ** (-0.40 * (zl / 0.07))
                    !
                    ! ... paj: compute phish for z0t over land
                    !
                    z0t = znt(i) / exp(czil * karman * sqrt(restar))
                    !
                    zolzz = zol(i) * (za(i) + z0t) / za(i)    ! (z+z0t)/l
                    zol10 = zol(i) * (10. + z0t) / za(i)   ! (10+z0t)/l
                    zol2 = zol(i) * (2. + z0t) / za(i)     ! (2+z0t)/l
                    zol0 = zol(i) * z0t / za(i)          ! z0t/l
                    !
                    if (zol(i).gt.0.) then
                        psih(i) = psih_stable(zolzz) - psih_stable(zol0)
                        psih10(i) = psih_stable(zol10) - psih_stable(zol0)
                        psih2(i) = psih_stable(zol2) - psih_stable(zol0)
                    else
                        if (zol(i).eq.0) then
                            psih(i) = 0.
                            psih10(i) = 0.
                            psih2(i) = 0.
                        else
                            psih(i) = psih_unstable(zolzz) - psih_unstable(zol0)
                            psih10(i) = psih_unstable(zol10) - psih_unstable(zol0)
                            psih2(i) = psih_unstable(zol2) - psih_unstable(zol0)
                        endif
                    endif
                    !
                    psiq = alog((za(i) + z0t) / z0t) - psih(i)
                    psiq2 = alog((2. + z0t) / z0t) - psih2(i)
                    psit = psiq
                    psit2 = psiq2
                    !
                    !              psit=gz1oz0(i)-psih(i)+czil*karman*sqrt(restar)
                    !              psiq=gz1oz0(i)-psih(i)+czil*karman*sqrt(restar)
                    !              psit2=gz2oz0(i)-psih2(i)+czil*karman*sqrt(restar)
                    !              psiq2=gz2oz0(i)-psih2(i)+czil*karman*sqrt(restar)

                endif
            endif
            ! to prevent oscillations average with old value
            ust(i) = 0.5 * ust(i) + 0.5 * karman * wspd(i) / psix
            ! tke coupling: compute ust without vconv for use in tke scheme
            wspdi(i) = sqrt(ux(i) * ux(i) + vx(i) * vx(i))
            if (present(ustm)) then
                ustm(i) = 0.5 * ustm(i) + 0.5 * karman * wspdi(i) / psix
            endif

            u10(i) = ux(i) * psix10 / psix
            v10(i) = vx(i) * psix10 / psix
            th2(i) = thgb(i) + dtg * psit2 / psit
            q2(i) = qsfc(i) + (qx(i) - qsfc(i)) * psiq2 / psiq
            t2(i) = th2(i) * (psfcpa(i) / p1000mb)**rovcp
            !
            if((xland(i) - 1.5).lt.0.)then
                ust(i) = amax1(ust(i), 0.001)
            endif
            mol(i) = karman * dtg / psit / prt
            denomq(i) = psiq
            denomq2(i) = psiq2
            denomt2(i) = psit2
            fm(i) = psix
            fh(i) = psit
        330 continue
        !
        335 continue

        !-----compute the surface sensible and latent heat fluxes:
        if (present(scm_force_flux)) then
            if (scm_force_flux.eq.1) goto 350
        endif
        do i = its, ite
            qfx(i) = 0.
            hfx(i) = 0.
        enddo
        350 continue

        if (isfflx.eq.0) goto 410

        !-----over water, alter roughness length (znt) according to wind (ust).

        do 360 i = its, ite
            if((xland(i) - 1.5).ge.0)then
                !         znt(i)=czo*ust(i)*ust(i)/g+ozo
                ! since v3.7 (ref: ec physics document for cy36r1)
                znt(i) = czo * ust(i) * ust(i) / g + 0.11 * 1.5e-5 / ust(i)
                ! v3.9: add limit as in isftcflx = 1,2
                znt(i) = min(znt(i), 2.85e-3)
                ! coare 3.5 (edson et al. 2013)
                !         czc = 0.0017*wspd(i)-0.005
                !         czc = min(czc,0.028)
                !         znt(i)=czc*ust(i)*ust(i)/g+0.11*1.5e-5/ust(i)
                ! ahw: change roughness length, and hence the drag coefficients ck and cd
                if (present(isftcflx)) then
                    if (isftcflx.ne.0) then
                        !               znt(i)=10.*exp(-9.*ust(i)**(-.3333))
                        !               znt(i)=10.*exp(-9.5*ust(i)**(-.3333))
                        !               znt(i)=znt(i) + 0.11*1.5e-5/amax1(ust(i),0.01)
                        !               znt(i)=0.011*ust(i)*ust(i)/g+ozo
                        !               znt(i)=max(znt(i),3.50e-5)
                        ! ahw 2012:
                        zw = min((ust(i) / 1.06)**(0.3), 1.0)
                        zn1 = 0.011 * ust(i) * ust(i) / g + ozo
                        zn2 = 10. * exp(-9.5 * ust(i)**(-.3333)) + &
                                0.11 * 1.5e-5 / amax1(ust(i), 0.01)
                        znt(i) = (1.0 - zw) * zn1 + zw * zn2
                        znt(i) = min(znt(i), 2.85e-3)
                        znt(i) = max(znt(i), 1.27e-7)
                    endif
                endif
                zl = znt(i)
            else
                zl = 0.01
            endif
            flqc(i) = rhox(i) * mavail(i) * ust(i) * karman / denomq(i)
            !       flqc(i)=rhox(i)*mavail(i)*ust(i)*karman/(   &
            !               alog(karman*ust(i)*za(i)/xka+za(i)/zl)-psih(i))
            dtthx = abs(thx(i) - thgb(i))
            if(dtthx.gt.1.e-5)then
                flhc(i) = cpm(i) * rhox(i) * ust(i) * mol(i) / (thx(i) - thgb(i))
                !         write(*,1001)flhc(i),cpm(i),rhox(i),ust(i),mol(i),thx(i),thgb(i),i
                1001   format(f8.5, 2x, f12.7, 2x, f12.10, 2x, f12.10, 2x, f13.10, 2x, f12.8, f12.8, 2x, i3)
            else
                flhc(i) = 0.
            endif
        360 continue

        !
        !-----compute surface moist flux:
        !
        !     if(idry.eq.1)goto 390
        !
        if (present(scm_force_flux)) then
            if (scm_force_flux.eq.1) goto 405
        endif

        do 370 i = its, ite
            qfx(i) = flqc(i) * (qsfc(i) - qx(i))
            qfx(i) = amax1(qfx(i), 0.)
            lh(i) = xlv * qfx(i)
        370 continue

        !-----compute surface heat flux:
        !
        390 continue
        do 400 i = its, ite
            if(xland(i) - 1.5.gt.0.)then
                hfx(i) = flhc(i) * (thgb(i) - thx(i))
                !         if ( present(isftcflx) ) then
                !            if ( isftcflx.ne.0 ) then
                ! ahw: add dissipative heating term (commented out in 3.6.1)
                !               hfx(i)=hfx(i)+rhox(i)*ustm(i)*ustm(i)*wspdi(i)
                !            endif
                !         endif
            elseif(xland(i) - 1.5.lt.0.)then
                hfx(i) = flhc(i) * (thgb(i) - thx(i))
                hfx(i) = amax1(hfx(i), -250.)
            endif
        400 continue

        405 continue

        do i = its, ite
            if((xland(i) - 1.5).ge.0)then
                zl = znt(i)
            else
                zl = 0.01
            endif
            !v3.1.1
            !         chs(i)=ust(i)*karman/(alog(karman*ust(i)*za(i) &
            !                /xka+za(i)/zl)-psih(i))
            chs(i) = ust(i) * karman / denomq(i)
            !        gz2oz0(i)=alog(2./znt(i))
            !        psim2(i)=-10.*gz2oz0(i)
            !        psim2(i)=amax1(psim2(i),-10.)
            !        psih2(i)=psim2(i)
            ! v3.1.1
            !         cqs2(i)=ust(i)*karman/(alog(karman*ust(i)*2.0  &
            !               /xka+2.0/zl)-psih2(i))
            !         chs2(i)=ust(i)*karman/(gz2oz0(i)-psih2(i))
            cqs2(i) = ust(i) * karman / denomq2(i)
            chs2(i) = ust(i) * karman / denomt2(i)
        enddo

        410 continue
        !jdf
        !     do i=its,ite
        !       if(ust(i).ge.0.1) then
        !         rmol(i)=rmol(i)*(-flhc(i))/(ust(i)*ust(i)*ust(i))
        !       else
        !         rmol(i)=rmol(i)*(-flhc(i))/(0.1*0.1*0.1)
        !       endif
        !     enddo
        !jdf

        !
    end subroutine sfclayrev1d

    !====================================================================
    subroutine sfclayrevinit

        integer :: n
        real :: zolf

        do n = 0, 1000
            ! stable function tables
            zolf = float(n) * 0.01
            psim_stab(n) = psim_stable_full(zolf)
            psih_stab(n) = psih_stable_full(zolf)

            ! unstable function tables
            zolf = -float(n) * 0.01
            psim_unstab(n) = psim_unstable_full(zolf)
            psih_unstab(n) = psih_unstable_full(zolf)

        enddo

    end subroutine sfclayrevinit

    function zolri(ri, z, z0)
        !
        if (ri.lt.0.)then
            x1 = -5.
            x2 = 0.
        else
            x1 = 0.
            x2 = 5.
        endif
        !
        fx1 = zolri2(x1, ri, z, z0)
        fx2 = zolri2(x2, ri, z, z0)
        do while (abs(x1 - x2) > 0.01)
            if(abs(fx2).lt.abs(fx1))then
                x1 = x1 - fx1 / (fx2 - fx1) * (x2 - x1)
                fx1 = zolri2(x1, ri, z, z0)
                zolri = x1
            else
                x2 = x2 - fx2 / (fx2 - fx1) * (x2 - x1)
                fx2 = zolri2(x2, ri, z, z0)
                zolri = x2
            endif
            !
        enddo
        !

        return
    end function

    !
    ! -----------------------------------------------------------------------
    !
    function zolri2(zol2, ri2, z, z0)
        !
        if(zol2 * ri2 .lt. 0.)zol2 = 0.  ! limit zol2 - must be same sign as ri2
        !
        zol20 = zol2 * z0 / z ! z0/l
        zol3 = zol2 + zol20 ! (z+z0)/l
        !
        if (ri2.lt.0) then
            psix2 = log((z + z0) / z0) - (psim_unstable(zol3) - psim_unstable(zol20))
            psih2 = log((z + z0) / z0) - (psih_unstable(zol3) - psih_unstable(zol20))
        else
            psix2 = log((z + z0) / z0) - (psim_stable(zol3) - psim_stable(zol20))
            psih2 = log((z + z0) / z0) - (psih_stable(zol3) - psih_stable(zol20))
        endif
        !
        zolri2 = zol2 * psih2 / psix2**2 - ri2
        !
        return
    end function
    !
    ! ... integrated similarity functions ...
    !
    function psim_stable_full(zolf)
        psim_stable_full = -6.1 * log(zolf + (1 + zolf**2.5)**(1. / 2.5))
        return
    end function

    function psih_stable_full(zolf)
        psih_stable_full = -5.3 * log(zolf + (1 + zolf**1.1)**(1. / 1.1))
        return
    end function

    function psim_unstable_full(zolf)
        x = (1. - 16. * zolf)**.25
        psimk = 2 * alog(0.5 * (1 + x)) + alog(0.5 * (1 + x * x)) - 2. * atan(x) + 2. * atan(1.)
        !
        ym = (1. - 10. * zolf)**0.33
        psimc = (3. / 2.) * log((ym**2. + ym + 1.) / 3.) - sqrt(3.) * atan((2. * ym + 1) / sqrt(3.)) + 4. * atan(1.) / sqrt(3.)
        !
        psim_unstable_full = (psimk + zolf**2 * (psimc)) / (1 + zolf**2.)

        return
    end function

    function psih_unstable_full(zolf)
        y = (1. - 16. * zolf)**.5
        psihk = 2. * log((1 + y) / 2.)
        !
        yh = (1. - 34. * zolf)**0.33
        psihc = (3. / 2.) * log((yh**2. + yh + 1.) / 3.) - sqrt(3.) * atan((2. * yh + 1) / sqrt(3.)) + 4. * atan(1.) / sqrt(3.)
        !
        psih_unstable_full = (psihk + zolf**2 * (psihc)) / (1 + zolf**2.)

        return
    end function

    ! look-up table functions
    function psim_stable(zolf)
        integer :: nzol
        real :: rzol
        nzol = int(zolf * 100.)
        rzol = zolf * 100. - nzol
        if(nzol + 1 .le. 1000)then
            psim_stable = psim_stab(nzol) + rzol * (psim_stab(nzol + 1) - psim_stab(nzol))
        else
            psim_stable = psim_stable_full(zolf)
        endif
        return
    end function

    function psih_stable(zolf)
        integer :: nzol
        real :: rzol
        nzol = int(zolf * 100.)
        rzol = zolf * 100. - nzol
        if(nzol + 1 .le. 1000)then
            psih_stable = psih_stab(nzol) + rzol * (psih_stab(nzol + 1) - psih_stab(nzol))
        else
            psih_stable = psih_stable_full(zolf)
        endif
        return
    end function

    function psim_unstable(zolf)
        integer :: nzol
        real :: rzol
        nzol = int(-zolf * 100.)
        rzol = -zolf * 100. - nzol
        if(nzol + 1 .le. 1000)then
            psim_unstable = psim_unstab(nzol) + rzol * (psim_unstab(nzol + 1) - psim_unstab(nzol))
        else
            psim_unstable = psim_unstable_full(zolf)
        endif
        return
    end function

    function psih_unstable(zolf)
        integer :: nzol
        real :: rzol
        nzol = int(-zolf * 100.)
        rzol = -zolf * 100. - nzol
        if(nzol + 1 .le. 1000)then
            psih_unstable = psih_unstab(nzol) + rzol * (psih_unstab(nzol + 1) - psih_unstab(nzol))
        else
            psih_unstable = psih_unstable_full(zolf)
        endif
        return
    end function

    !-------------------------------------------------------------------

end module module_sf_sfclayrev

!
! ----------------------------------------------------------
!

max_stars_repo_name: NCAR/mlsurfacelayer

max_issues_repo_name: NCAR/mlsurfacelayer

max_forks_repo_name: NCAR/mlsurfacelayer


text:
      subroutine fcnkb(t,dsize,drdt)
      character*100 id,revision,revdate
      data id/'$RCSfile: fcnkb.f,v $'/
      data revision/'$Revision: 1.7 $'/
      data revdate/'$Date: 1994/03/14 00:34:07 $'/
      include "const.h"
      include "salcon.h"
      include "arrays.h"
      include "interp.h"

      real*8 dsize(*),drdt(*),t
      real*8 sw,es,ew,qs,temp,falpha,fbeta,s,cr
      real*8 rhol,c,e,rho,ck,c2,c7,denom,dqvdt

      integer i

c
c     initialize diffusivities D,Ktemp and lv through
c     common block const.h
c
      nrhs=nrhs+1
      call constants(dsize(dmax+2)-273.15,dsize(dmax+4))
      sw=0.
      temp=dsize(dmax+2)
      es=ew(dsize(dmax+2))*100.
      qs=eps*es/(dsize(dmax+4) - es)
      dqvdt=0.
      do 10 i=1,dmax
         falpha=dsize(i)/(dsize(i) + lalpha)
         fbeta=dsize(i)/(dsize(i) + lbeta)
         if(i .gt. dmaxa) then
            c2=c2nacl
            c7=c7nacl
         else
            c2=c2am
            c7=c7am
         endif
c
c     check to see whether droplet is < 0.2 um, has a radius
c     less than rcrit and the saturation is less than scrit.
c     then get the equilibrium radius using the interpolator
c
         S=dsize(dmax+3)
c        if((dsize(i) .lt. 2.e-7) .and.
c    $        (dsize(i) .lt. splinearray(i,2,npoints)) .and.
c    $        (1+S) .lt. splinearray(i,1,npoints)) then
c           call psplint(i,1+S,dsize(i))
c write(*,*)i,dsize(i),inter
c           drdt(i)=0.
c        else
            rhol=(dsize(i)**3*pi43*1.e3+csize(i)*c7)/(dsize(i)**3*pi43)
            Ck= (2*sigma/(Rv*dsize(dmax+2)*rhol*dsize(i)))
            Cr= c2*csize(i)/(pi43*dsize(i)**3*rhol-csize(i))
            denom=rhol*
     $           (Rv*temp/(fbeta*D*es)+Lv**2/(falpha*Ktemp*Rv*temp**2))
            drdt(i)=1./dsize(i)*(S - Ck + Cr)/denom
c        endif
         dqvdt=dqvdt-pi4*rnr(i)*dsize(i)**2*drdt(i)*1.e3
 10   continue
c
c
c
      e=dsize(dmax+1)*dsize(dmax+4)/(eps + dsize(dmax+1))
      rho=(dsize(dmax+4) - e)/(rd*dsize(dmax+2))*(1.+
     c     dsize(dmax+7)) + e/(rv*dsize(dmax+2))
c     c=cpd+dsize(dmax+1)*cv+dsize(dmax+7)*cw
	c=1005.
      drdt(dmax+1)=dqvdt
      drdt(dmax+2)=-lv/c*drdt(dmax+1)
c
c try to take out the velocity effect
c
c     -g/c*dsize(dmax+5)
c

c
c old version
c
c$$$      drdt(dmax+3)=dsize(dmax+4)*rv/(es*rd)*drdt(dmax+1)
c$$$     c     -(1.+dsize(dmax+3))
c$$$     c     *(lv/(rv*dsize(dmax+2)**2)*drdt(dmax+2)
c$$$     c     +g/(rd*dsize(dmax+2))*dsize(dmax+5))
c
c new version  (94/10/26)
c
      drdt(dmax+3)=(1.+dsize(dmax+3))*(Rd/(dsize(dmax+1)*
     c (Rv*dsize(dmax+1)+Rd))*drdt(dmax+1)
     c     -g*rho*dsize(dmax+5)/dsize(dmax+4)
     c     -lv/(Rv*dsize(dmax+2)**2)*drdt(dmax+2))
      drdt(dmax+4)= -1.*rho*g*dsize(dmax+5)
      drdt(dmax+5)=0.
      drdt(dmax+6)=dsize(dmax+5)
      drdt(dmax+7)= -1.*dqvdt
      return
      end

ext: f

max_stars_repo_name: phaustin/empm

max_issues_repo_name: phaustin/empm

max_forks_repo_name: phaustin/empm


text:
program test
  integer, parameter :: b = 80
  character(len=b) :: msg1 = "hello world"
  character(len=80) :: msg2
  character :: ch1 = "c"
  character :: ch2
  integer :: a
  msg2 = "Hello world"
  ch2 = "C"
  a = 6 + 4
end program test

max_stars_repo_name: clementval/fc

max_issues_repo_name: clementval/fc

max_forks_repo_name: clementval/fc


text:
      subroutine qqb_w_twdk_gs(p,msq)
************************************************************************
*     Author: J.M. Campbell                                            *
*     January, 2005.                                                   *
************************************************************************
c---Matrix element SUBTRACTION squared averaged over initial colors and spins
c     q(-p1)+qbar(-p2) -->  W + t [or t~](p5) + g(p6)
c                           |
c                            -->l(p3)+a(p4)

      implicit none
      include 'constants.f'
      include 'ptilde.f'
      include 'qqgg.f'
      include 'nores.f'
      include 'nwz.f'
      integer j,k,nd,iq
      double precision p(mxpart,4),msq(maxd,-nf:nf,-nf:nf)
      double precision
     & msq16_2(-nf:nf,-nf:nf),msq26_1(-nf:nf,-nf:nf),
     & msq26_5(-nf:nf,-nf:nf),msq16_5(-nf:nf,-nf:nf),
     & msq56_1(-nf:nf,-nf:nf),msq56_2(-nf:nf,-nf:nf),
     & msq26_5v(-nf:nf,-nf:nf),msq26_1v(-nf:nf,-nf:nf),
     & msq16_2v(-nf:nf,-nf:nf),msq16_5v(-nf:nf,-nf:nf),
     & dummyv(-nf:nf,-nf:nf),
     & sub16_2(4),sub26_1(4),sub16_5(4),sub26_5(4),
     & sub56_1(4),sub56_2(4),sub56_1v,sub56_2v,
     & sub26_5v,sub26_1v,sub16_5v,sub16_2v
      external qqb_w_twdk,qqb_w_twdk_gvec,donothing_gvec

c--- Note that the subtractions here are very similar to the ones
c--- for W+1 jet, except that the top mass means that we must
c--- separate the initial-final and final-initial dipoles

c--- set up lepton variables depending on whether it's t or tbar
      if     (nwz .eq. -1) then
        iq=+1 ! quark initial state
      elseif (nwz .eq. +1) then
        iq=-1 ! antiquark initial state
      else
        write(6,*) 'Error in qqb_w_twdk_gs, nwz is not +1 or -1 :  ',nwz
      stop
      endif

      ndmax=6

      qqproc=.true.
      qgproc=.false.
      gqproc=.false.
      ggproc=.false.

c--- calculate the initial-initial dipoles
      call dips_mass(3,p,1,6,2,sub16_2,sub16_2v,msq16_2,msq16_2v,
     . qqb_w_twdk,qqb_w_twdk_gvec)
      call dips_mass(4,p,2,6,1,sub26_1,sub26_1v,msq26_1,msq26_1v,
     . qqb_w_twdk,qqb_w_twdk_gvec)

c--- now the initial-final ones
      call dips_mass(5,p,1,6,5,sub16_5,sub16_5v,msq16_5,msq16_5v,
     . qqb_w_twdk,qqb_w_twdk_gvec)
      call dips_mass(6,p,2,6,5,sub26_5,sub26_5v,msq26_5,msq26_5v,
     . qqb_w_twdk,qqb_w_twdk_gvec)

c--- now the final-initial ones
      call dips_mass(1,p,5,6,1,sub56_1,sub56_1v,msq56_1,dummyv,
     . qqb_w_twdk,donothing_gvec)
      call dips_mass(2,p,5,6,2,sub56_2,sub56_2v,msq56_2,dummyv,
     . qqb_w_twdk,donothing_gvec)

      do j=-nf,nf
      do k=-nf,nf
      do nd=1,ndmax
        msq(nd,j,k)=0d0
      enddo
      enddo
      enddo

      do j=-nf,nf,nf
      do k=-nf,nf,nf

      if     ((k .eq. 0).and. (j .ne. 0)) then
c--- q-g and qb-g cases
      msq(3,j,k)=xn*msq16_2(j,k)*sub16_2(qq)
      msq(4,j,k)=xn*(msq26_1(j,k)*sub26_1(gg)+msq26_1v(j,k)*sub26_1v)
      msq(5,j,k)=-msq16_5(j,k)*sub16_5(qq)/xn
      msq(1,j,k)=-msq56_1(j,k)*sub56_1(qq)/xn
      msq(6,j,k)=xn*(msq26_5(j,k)*sub26_5(gg)+msq26_5v(j,k)*sub26_5v)
      msq(2,j,k)=xn*msq56_2(j,k)*sub56_2(qq)

      elseif ((j .eq. 0).and.(k.ne.0)) then
c--- g-q and g-qb cases
      msq(3,j,k)=xn*(msq16_2(j,k)*sub16_2(gg)+msq16_2v(j,k)*sub16_2v)
      msq(4,j,k)=xn*msq26_1(j,k)*sub26_1(qq)
      msq(5,j,k)=xn*(msq16_5(j,k)*sub16_5(gg)+msq16_5v(j,k)*sub16_5v)
      msq(1,j,k)=xn*msq56_1(j,k)*sub56_1(qq)
      msq(6,j,k)=-msq26_5(j,k)*sub26_5(qq)/xn
      msq(2,j,k)=-msq56_2(j,k)*sub56_2(qq)/xn

      elseif ((j .eq. 0).and.(k .eq. 0)) then
      if (nores .eqv. .false.) then
      msq(3,j,k)=(msq16_2(-5,k)+msq16_2(+5,k))*sub16_2(qg)*2d0*tr
      msq(4,j,k)=(msq26_1(k,-5)+msq26_1(k,+5))*sub26_1(qg)*2d0*tr
      endif

      endif

      enddo
      enddo

      do j=1,4
c--- subtraction terms for 4-quark matrix elements
      msq(3,j,+5*iq)=2d0*cf*(
     .          msq16_2(0,+5*iq)*sub16_2(gq)+msq16_2v(0,+5*iq)*sub16_2v)
      msq(4,+5*iq,j)=2d0*cf*(
     .          msq26_1(+5*iq,0)*sub26_1(gq)+msq26_1v(+5*iq,0)*sub26_1v)

      msq(3,-j,+5*iq)=msq(3,+j,+5*iq)
      msq(4,+5*iq,-j)=msq(4,+5*iq,+j)
      enddo

      return
      end

max_stars_repo_name: tmartini/JHUGen

max_issues_repo_name: tmartini/JHUGen

max_forks_repo_name: tmartini/JHUGen


text:
C Various CERNLIB mathematical routines needed for the optimum interval,
C as modified from 2006_src.tar.gz to not call other cernlib routines.
C See http://cernlib.web.cern.ch/cernlib/download/2006_source/tar/
C
      FUNCTION DFREQ(X)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
C
C
      DIMENSION P1(0:3),Q1(0:3),P2(0:7),Q2(0:7),P3(0:4),Q3(0:4)

      PARAMETER(Z1 = 1, HF = Z1/2)
      PARAMETER(C1 = 0.56418 95835 47756 29D0)
      PARAMETER(W2 = 1.41421 35623 73095 05D0, RW2 = 1/W2)

      DATA (P1(I),Q1(I),I=0,3)
     +/+2.42667 95523 05317 5D+2, +2.15058 87586 98612 0D+2,
     1 +2.19792 61618 29415 2D+1, +9.11649 05404 51490 1D+1,
     2 +6.99638 34886 19135 5D+0, +1.50827 97630 40778 7D+1,
     3 -3.56098 43701 81538 5D-2, +1/

      DATA (P2(I),Q2(I),I=0,7)
     +/+3.00459 26102 01616 01D+2, +3.00459 26095 69832 93D+2,
     1 +4.51918 95371 18729 42D+2, +7.90950 92532 78980 27D+2,
     2 +3.39320 81673 43436 87D+2, +9.31354 09485 06096 21D+2,
     3 +1.52989 28504 69404 04D+2, +6.38980 26446 56311 67D+2,
     4 +4.31622 27222 05673 53D+1, +2.77585 44474 39876 43D+2,
     5 +7.21175 82508 83093 66D+0, +7.70001 52935 22947 30D+1,
     6 +5.64195 51747 89739 71D-1, +1.27827 27319 62942 35D+1,
     7 -1.36864 85738 27167 07D-7, +1/

      DATA (P3(I),Q3(I),I=0,4)
     +/-2.99610 70770 35421 74D-3, +1.06209 23052 84679 18D-2,
     1 -4.94730 91062 32507 34D-2, +1.91308 92610 78298 41D-1,
     2 -2.26956 59353 96869 30D-1, +1.05167 51070 67932 07D+0,
     3 -2.78661 30860 96477 88D-1, +1.98733 20181 71352 56D+0,
     4 -2.23192 45973 41846 86D-2, +1/

      V=RW2*ABS(X)
      IF(V .LT. HF) THEN
       Y=V**2
       AP=P1(3)
       AQ=Q1(3)
       DO 1 I = 2,0,-1
       AP=P1(I)+Y*AP
    1  AQ=Q1(I)+Y*AQ
       H=V*AP/AQ
       HC=1-H
      ELSEIF(V .LT. 4) THEN
       AP=P2(7)
       AQ=Q2(7)
       DO 2 I = 6,0,-1
       AP=P2(I)+V*AP
    2  AQ=Q2(I)+V*AQ
       HC=EXP(-V**2)*AP/AQ
       H=1-HC
      ELSE
       Y=1/V**2
       AP=P3(4)
       AQ=Q3(4)
       DO 3 I = 3,0,-1
       AP=P3(I)+Y*AP
    3  AQ=Q3(I)+Y*AQ
       HC=EXP(-V**2)*(C1+Y*AP/AQ)/V
       H=1-HC
      ENDIF
      IF(X .GT. 0) THEN
       DFREQ=HF+HF*H
      ELSE
       DFREQ=HF*HC
      ENDIF
      RETURN
      END
      FUNCTION DGAUSN(P)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C     Computes a "Normal Deviate"
C     Based on G.W. Hill & A.W. Davis, Algorithm 442 Normal Deviate
C     Collected Algorithms from CACM

C
      PARAMETER (C = 2.50662 82746 31000 50D0)
      PARAMETER (Z1 = 1, HF = Z1/2, C1 = 3*Z1/4, C2 = 7*Z1/8, C3 = Z1/3)

      IF(P .LE. 0 .OR. P .GE. 1) THEN
       H=0
       WRITE(6,101) P
      ELSEIF(P .EQ. HF) THEN
       H=0
      ELSE
       X=P
       IF(P .GT. HF) X=1-P
       X=SQRT(-2*LOG(X))
       X=X-((7.47395*X+494.877)*X+1637.720)/
     1     (((X+117.9407)*X+908.401)*X+659.935)
       IF(P .LT. HF) X=-X
       S=X**2
       Z=C*(P-DFREQ(X))*EXP(HF*S)
       H=(((((C1*S+C2)*Z+X)*X+HF)*C3*Z+HF*X)*Z+1)*Z+X
      ENDIF
      DGAUSN=H
      RETURN
  101 FORMAT('ARGUMENT P =',1P,D15.8,' NOT IN RANGE')
      END

      FUNCTION GAMDIS(X,A)
      IMPLICIT REAL (A-H,O-Z)
C     Calculates the gamma distribution function
C
C        G(x,a) = (1/gamma(a)) * int(0,x)[exp(-t) * t**(a-1)] dt.
C
C     Based on
C        W. Gautschi, ALGORITHM 542 Incomplete Gamma Functions,
C        ACM Trans. Math. Software 5 (1979) 482-489.

      DIMENSION C(14)

      PARAMETER (EPS = 1E-5, EPS1 = 5E-7)
      PARAMETER (ALH = -0.69314 72)
      PARAMETER (Z1 = 1, HALF = Z1/2, QUAR = Z1/4)
      PARAMETER (C1 = 3*Z1/2, KMAX = 300)

      DATA C
     1/ 0.5772157,-0.6558781,-0.0420026, 0.1665386,-0.0421977,
     2 -0.0096220, 0.0072189,-0.0011652,-0.0002152, 0.0001281,
     3 -0.0000201,-0.0000013, 0.0000011,-0.0000002/

      HST=0
      IF(X .EQ. 0) GO TO 99
      IF(X .LT. 0 .OR. A .LE. 0) THEN
       WRITE(6,101) X,A
       GO TO 99
      ELSE
       ALX=LOG(X)
      ENDIF
      IF(X .LT. QUAR) THEN
       ALFA=ALH/ALX
      ELSE
       ALFA=X+QUAR
      ENDIF
      IF(A .GT. ALFA) THEN
       TERM=1
       SUM=1
       DO 1 K = 1,KMAX
       TERM=X*TERM/(A+K)
       SUM=SUM+TERM
       IF(ABS(TERM) .LE. EPS*SUM) GO TO 2
    1  CONTINUE
       GO TO 98
    2  HST=SUM*EXP(A*ALX-X-ALOGAM(1+A))
      ELSEIF(X .GT. C1) THEN
       P=0
       S=1-A
       Q=(X+S)*(X-1-A)
       R=4*(X+S)
       TERM=1
       SUM=1
       RHO=0
       DO 3 K = 2,KMAX
       P=P+S
       Q=Q+R
       R=R+8
       S=S+2
       T=P*(1+RHO)
       RHO=T/(Q-T)
       TERM=RHO*TERM
       SUM=SUM+TERM
       IF(ABS(TERM) .LE. EPS*SUM) GO TO 4
    3  CONTINUE
       GO TO 98
    4  HST=1-(A*SUM/(X+1-A))*EXP(A*ALX-X-ALOGAM(1+A))
      ELSE
       IF(A .LT. HALF) THEN
        SUM=C(14)
        DO 12 K = 13,1,-1
   12   SUM=A*SUM+C(K)
        GA=-SUM/(1+A*SUM)
        Y=A*ALX
        IF(ABS(Y) .GE. 1) THEN
         U=GA-(EXP(Y)-1)/A
        ELSE
         SUM=1
         TERM=1
         DO 7 K = 2,KMAX
         TERM=Y*TERM/K
         SUM=SUM+TERM
         IF(ABS(TERM) .LE. EPS1*SUM) GO TO 8
    7    CONTINUE
         GO TO 98
    8    U=GA-SUM*ALX
        ENDIF
       ELSE
        U=GAMMA(A)-EXP(A*ALX)/A
       ENDIF
       P=A*X
       Q=A+1
       R=A+3
       TERM=1
       SUM=1
       DO 9 K = 2,KMAX
       P=P+X
       Q=Q+R
       R=R+2
       TERM=-P*TERM/Q
       SUM=SUM+TERM
       IF(ABS(TERM) .LE. EPS1*SUM) GO TO 10
    9  CONTINUE
       GO TO 98
   10  HST=1-A*(U+SUM*EXP((1+A)*ALX)/(1+A))/GAMMA(1+A)
      ENDIF
   99 GAMDIS=HST
      RETURN

   98 WRITE(6,102) X,A
      GO TO 99
  101 FORMAT('ILLEGAL ARGUMENT(S) X = ',E15.8,'  A = ',E15.8)
  102 FORMAT('PROBLEMS WITH CONVERGENCE, X = ',E15.8,'  A = ',E15.8)
      END
      SUBROUTINE RZERO(A,B,X0,R,EPS,MXF,F)
      IMPLICIT REAL (A-H,O-Z)

      EXTERNAL F

      PARAMETER (ONE = 1, HALF = ONE/2)

      XA=MIN(A,B)
      XB=MAX(A,B)
      FA=F(XA,1)
      FB=F(XB,2)
      IF(FA*FB .GT. 0) GO TO 5
      MC=0

    1 X0=HALF*(XA+XB)
      R=X0-XA
      EE=EPS*(ABS(X0)+1)
      IF(R .LE. EE) GO TO 4
      F1=FA
      X1=XA
      F2=FB
      X2=XB

    2 FX=F(X0,2)
      MC=MC+1
      IF(MC .GT. MXF) GO TO 6
      IF(FX*FA .GT. 0) THEN
       XA=X0
       FA=FX
      ELSE
       XB=X0
       FB=FX
      END IF

    3 U1=F1-F2
      U2=X1-X2
      U3=F2-FX
      U4=X2-X0
      IF(U2 .EQ. 0 .OR. U4 .EQ. 0) GO TO 1
      F3=FX
      X3=X0
      U1=U1/U2
      U2=U3/U4
      CA=U1-U2
      CB=(X1+X2)*U2-(X2+X0)*U1
      CC=(X1-X0)*F1-X1*(CA*X1+CB)
      IF(CA .EQ. 0) THEN
       IF(CB .EQ. 0) GO TO 1
       X0=-CC/CB
      ELSE
       U3=CB/(2*CA)
       U4=U3*U3-CC/CA
       IF(U4 .LT. 0) GO TO 1
       X0=-U3+SIGN(SQRT(U4),X0+U3)
      END IF
      IF(X0 .LT. XA .OR. X0 .GT. XB) GO TO 1

      R=MIN(ABS(X0-X3),ABS(X0-X2))
      EE=EPS*(ABS(X0)+1)
      IF(R .GT. EE) THEN
       F1=F2
       X1=X2
       F2=F3
       X2=X3
       GO TO 2
      END IF

      FX=F(X0,2)
      IF(FX .EQ. 0) GO TO 4
      IF(FX*FA .LT. 0) THEN
       XX=X0-EE
       IF(XX .LE. XA) GO TO 4
       FF=F(XX,2)
       FB=FF
       XB=XX
      ELSE
       XX=X0+EE
       IF(XX .GE. XB) GO TO 4
       FF=F(XX,2)
       FA=FF
       XA=XX
      END IF
      IF(FX*FF .GT. 0) THEN
       MC=MC+2
       IF(MC .GT. MXF) GO TO 6
       F1=F3
       X1=X3
       F2=FX
       X2=X0
       X0=XX
       FX=FF
       GO TO 3
      END IF

    4 R=EE
      FF=F(X0,3)
      RETURN
    5 Continue
      WRITE(*,100)
      R=-2*(XB-XA)
      X0=0
      RETURN
    6 Continue
      WRITE(*,101)
      R=-HALF*ABS(XB-XA)
      X0=0
      RETURN
  100 FORMAT(1X,'***** Modified RZERO ... F(A) AND F(B)',
     1          ' HAVE THE SAME SIGN')
  101 FORMAT(1X,'***** Modified RZERO ... TOO MANY FUNCTION CALLS')
      END

      FUNCTION ALGAMA(X)
C
C
      DIMENSION P1(7),Q1(7),P2(7),Q2(7),P3(7),Q3(7),C(5)

      PARAMETER (Z1 = 1, HF = Z1/2, HF1 = 1+HF)
      DATA P1
     1/+3.84287 36567 45991D+0, +5.27068 93753 00983D+1,
     2 +5.55840 45723 51531D+1, -2.15135 13573 72570D+2,
     3 -2.45872 61722 29242D+2, -5.75008 93603 04123D+1,
     4 -2.33590 98949 51284D+0/
      DATA Q1
     1/+1.00000 00000 00000D+0, +3.37330 47907 07074D+1,
     2 +1.93877 84034 37713D+2, +3.08829 54973 42428D+2,
     3 +1.50068 39064 89095D+2, +2.01068 51344 33395D+1,
     4 +4.57174 20282 50299D-1/
      DATA P2
     1/+4.87402 01396 83863 6D+0, +2.48845 25168 57407 6D+2,
     2 +2.17973 66058 89591 5D+3, +3.79751 24011 52511 8D+3,
     3 -1.97780 70769 84164 6D+3, -3.69298 34005 59128 2D+3,
     4 -5.60177 73537 80387 7D+2/
      DATA Q2
     1/+1.00000 00000 00000 0D+0, +9.50999 17418 20893 8D+1,
     2 +1.56120 45277 92863 5D+3, +7.23400 87928 94807 1D+3,
     3 +1.04595 76594 05895 9D+4, +4.16994 15153 20023 1D+3,
     4 +2.76785 83623 80410 1D+2/
      DATA P3
     1/-6.88062 40094 59425D+3, -4.30699 69819 57098D+5,
     2 -4.75045 94653 43956D+6, -2.94234 45930 32234D+6,
     3 +3.63218 04931 54257D+7, -3.35677 82814 54576D+6,
     4 -2.48043 69488 28593D+7/
      DATA Q3
     1/+1.00000 00000 00000D+0, -1.42168 29839 65146D+3,
     2 -1.55528 90280 85353D+5, -3.41525 17108 01107D+6,
     3 -2.09696 23255 80444D+7, -3.45441 75093 34395D+7,
     4 -9.16055 82863 71317D+6/
      DATA C
     1/ 1.12249 21356 561D-1,  7.95916 92961 204D-2,
     1 -1.70877 94611 020D-3,  9.18938 53320 467D-1,
     2  1.34699 05627 879D+0/

      ENTRY ALOGAM(X)

      IF(X .LE. 0) THEN
       H=0
       WRITE(6,101) X
      ELSE IF(X .EQ. 1 .OR. X .EQ. 2) THEN
       H=0
      ELSE IF(X .LE. HF) THEN
       Y=X+1
       AP=P1(1)
       AQ=Q1(1)
       DO 2 I = 2,7
       AP=P1(I)+Y*AP
    2  AQ=Q1(I)+Y*AQ
       H=-LOG(X)+X*AP/AQ
      ELSE IF(X .LE. HF1) THEN
       AP=P1(1)
       AQ=Q1(1)
       DO 3 I = 2,7
       AP=P1(I)+X*AP
    3  AQ=Q1(I)+X*AQ
       H=(X-1)*AP/AQ
      ELSE IF(X .LE. 4) THEN
       AP=P2(1)
       AQ=Q2(1)
       DO 4 I = 2,7
       AP=P2(I)+X*AP
    4  AQ=Q2(I)+X*AQ
       H=(X-2)*AP/AQ
      ELSE IF(X .LE. 12) THEN
       AP=P3(1)
       AQ=Q3(1)
       DO 5 I = 2,7
       AP=P3(I)+X*AP
    5  AQ=Q3(I)+X*AQ
       H=AP/AQ
      ELSE
       Y=1/X**2
       H=(X-HF)*LOG(X)-X+C(4)+(C(1)+Y*(C(2)+Y*C(3)))/
     1                                        ((C(5)+Y)*X)
      ENDIF
      ALGAMA=H
      RETURN
  101 FORMAT('NON-POSITIVE ARGUMENT  X = ',1P,E15.6)
      END
      FUNCTION GAMMA(X)
      IMPLICIT REAL (A-H,O-Z)
C
      DIMENSION C(0:15)

      DATA C( 0) /3.65738 77250 83382 44D0/
      DATA C( 1) /1.95754 34566 61268 27D0/
      DATA C( 2) /0.33829 71138 26160 39D0/
      DATA C( 3) /0.04208 95127 65575 49D0/
      DATA C( 4) /0.00428 76504 82129 09D0/
      DATA C( 5) /0.00036 52121 69294 62D0/
      DATA C( 6) /0.00002 74006 42226 42D0/
      DATA C( 7) /0.00000 18124 02333 65D0/
      DATA C( 8) /0.00000 01096 57758 66D0/
      DATA C( 9) /0.00000 00059 87184 05D0/
      DATA C(10) /0.00000 00003 07690 81D0/
      DATA C(11) /0.00000 00000 14317 93D0/
      DATA C(12) /0.00000 00000 00651 09D0/
      DATA C(13) /0.00000 00000 00025 96D0/
      DATA C(14) /0.00000 00000 00001 11D0/
      DATA C(15) /0.00000 00000 00000 04D0/

      U=X
      IF(U .LE. 0) THEN
       WRITE(6,101) U
       H=0
       B0=0
       B2=0
       F=0
       GO TO 9
      ENDIF
    8 F=1
      IF(U .LT. 3) THEN
       DO 1 I = 1,INT(4-U)
       F=F/U
    1  U=U+1
      ELSE
       DO 2 I = 1,INT(U-3)
       U=U-1
    2  F=F*U
      END IF
      H=U+U-7
      ALFA=H+H
      B1=0
      B2=0
      DO 3 I = 15,0,-1
      B0=C(I)+ALFA*B1-B2
      B2=B1
    3 B1=B0
    9 GAMMA =F*(B0-H*B2)
      RETURN
  101 FORMAT('ARGUMENT IS NEGATIVE = ',1P,E15.1)
      END
      SUBROUTINE RANLUX(RVEC,LENV)
C         Subtract-and-borrow random number generator proposed by
C         Marsaglia and Zaman, implemented by F. James with the name
C         RCARRY in 1991, and later improved by Martin Luescher
C         in 1993 to produce "Luxury Pseudorandom Numbers".
C     Fortran 77 coded by F. James, 1993
C
C   LUXURY LEVELS.
C   ------ ------      The available luxury levels are:
C
C  level 0  (p=24): equivalent to the original RCARRY of Marsaglia
C           and Zaman, very long period, but fails many tests.
C  level 1  (p=48): considerable improvement in quality over level 0,
C           now passes the gap test, but still fails spectral test.
C  level 2  (p=97): passes all known tests, but theoretically still
C           defective.
C  level 3  (p=223): DEFAULT VALUE.  Any theoretically possible
C           correlations have very small chance of being observed.
C  level 4  (p=389): highest possible luxury, all 24 bits chaotic.
C
C!!! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C!!!  Calling sequences for RANLUX:                                  ++
C!!!      CALL RANLUX (RVEC, LEN)   returns a vector RVEC of LEN     ++
C!!!                   32-bit random floating point numbers between  ++
C!!!                   zero (not included) and one (also not incl.). ++
C!!!      CALL RLUXGO(LUX,INT,K1,K2) initializes the generator from  ++
C!!!               one 32-bit integer INT and sets Luxury Level LUX  ++
C!!!               which is integer between zero and MAXLEV, or if   ++
C!!!               LUX .GT. 24, it sets p=LUX directly.  K1 and K2   ++
C!!!               should be set to zero unless restarting at a break++
C!!!               point given by output of RLUXAT (see RLUXAT).     ++
C!!!      CALL RLUXAT(LUX,INT,K1,K2) gets the values of four integers++
C!!!               which can be used to restart the RANLUX generator ++
C!!!               at the current point by calling RLUXGO.  K1 and K2++
C!!!               specify how many numbers were generated since the ++
C!!!               initialization with LUX and INT.  The restarting  ++
C!!!               skips over  K1+K2*E9   numbers, so it can be long.++
C!!!   A more efficient but less convenient way of restarting is by: ++
C!!!      CALL RLUXIN(ISVEC)    restarts the generator from vector   ++
C!!!                   ISVEC of 25 32-bit integers (see RLUXUT)      ++
C!!!      CALL RLUXUT(ISVEC)    outputs the current values of the 25 ++
C!!!                 32-bit integer seeds, to be used for restarting ++
C!!!      ISVEC must be dimensioned 25 in the calling program        ++
C!!! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
      DIMENSION RVEC(LENV)
      DIMENSION SEEDS(24), ISEEDS(24), ISDEXT(25)
      PARAMETER (MAXLEV=4, LXDFLT=3)
      DIMENSION NDSKIP(0:MAXLEV)
      DIMENSION NEXT(24)
      PARAMETER (TWOP12=4096., IGIGA=1000000000,JSDFLT=314159265)
      PARAMETER (ITWO24=2**24, ICONS=2147483563)
      SAVE NOTYET, I24, J24, CARRY, SEEDS, TWOM24, TWOM12, LUXLEV
      SAVE NSKIP, NDSKIP, IN24, NEXT, KOUNT, MKOUNT, INSEED
      INTEGER LUXLEV
      LOGICAL NOTYET
      DATA NOTYET, LUXLEV, IN24, KOUNT, MKOUNT /.TRUE., LXDFLT, 0,0,0/
      DATA I24,J24,CARRY/24,10,0./
C                               default
C  Luxury Level   0     1     2   *3*    4
      DATA NDSKIP/0,   24,   73,  199,  365 /
Corresponds to p=24    48    97   223   389
C     time factor 1     2     3     6    10   on slow workstation
C                 1    1.5    2     3     5   on fast mainframe
C
C  NOTYET is .TRUE. if no initialization has been performed yet.
C              Default Initialization by Multiplicative Congruential
      IF (NOTYET) THEN
         NOTYET = .FALSE.
         JSEED = JSDFLT
         INSEED = JSEED
         WRITE(6,'(A,I12)') ' RANLUX DEFAULT INITIALIZATION: ',JSEED
         LUXLEV = LXDFLT
         NSKIP = NDSKIP(LUXLEV)
         LP = NSKIP + 24
         IN24 = 0
         KOUNT = 0
         MKOUNT = 0
         WRITE(6,'(A,I2,A,I4)')  ' RANLUX DEFAULT LUXURY LEVEL =  ',
     +        LUXLEV,'      p =',LP
            TWOM24 = 1.
         DO 25 I= 1, 24
            TWOM24 = TWOM24 * 0.5
         K = JSEED/53668
         JSEED = 40014*(JSEED-K*53668) -K*12211
         IF (JSEED .LT. 0)  JSEED = JSEED+ICONS
         ISEEDS(I) = MOD(JSEED,ITWO24)
   25    CONTINUE
         TWOM12 = TWOM24 * 4096.
         DO 50 I= 1,24
         SEEDS(I) = REAL(ISEEDS(I))*TWOM24
         NEXT(I) = I-1
   50    CONTINUE
         NEXT(1) = 24
         I24 = 24
         J24 = 10
         CARRY = 0.
         IF (SEEDS(24) .EQ. 0.) CARRY = TWOM24
      ENDIF
C
C          The Generator proper: "Subtract-with-borrow",
C          as proposed by Marsaglia and Zaman,
C          Florida State University, March, 1989
C
      DO 100 IVEC= 1, LENV
      UNI = SEEDS(J24) - SEEDS(I24) - CARRY
      IF (UNI .LT. 0.)  THEN
         UNI = UNI + 1.0
         CARRY = TWOM24
      ELSE
         CARRY = 0.
      ENDIF
      SEEDS(I24) = UNI
      I24 = NEXT(I24)
      J24 = NEXT(J24)
      RVEC(IVEC) = UNI
C  small numbers (with less than 12 "significant" bits) are "padded".
      IF (UNI .LT. TWOM12)  THEN
         RVEC(IVEC) = RVEC(IVEC) + TWOM24*SEEDS(J24)
C        and zero is forbidden in case someone takes a logarithm
         IF (RVEC(IVEC) .EQ. 0.)  RVEC(IVEC) = TWOM24*TWOM24
      ENDIF
C        Skipping to luxury.  As proposed by Martin Luscher.
      IN24 = IN24 + 1
      IF (IN24 .EQ. 24)  THEN
         IN24 = 0
         KOUNT = KOUNT + NSKIP
         DO 90 ISK= 1, NSKIP
         UNI = SEEDS(J24) - SEEDS(I24) - CARRY
         IF (UNI .LT. 0.)  THEN
            UNI = UNI + 1.0
            CARRY = TWOM24
         ELSE
            CARRY = 0.
         ENDIF
         SEEDS(I24) = UNI
         I24 = NEXT(I24)
         J24 = NEXT(J24)
   90    CONTINUE
      ENDIF
  100 CONTINUE
      KOUNT = KOUNT + LENV
      IF (KOUNT .GE. IGIGA)  THEN
         MKOUNT = MKOUNT + 1
         KOUNT = KOUNT - IGIGA
      ENDIF
      RETURN
C
C           Entry to input and float integer seeds from previous run
      ENTRY RLUXIN(ISDEXT)
         NOTYET = .FALSE.
         TWOM24 = 1.
         DO 195 I= 1, 24
         NEXT(I) = I-1
  195    TWOM24 = TWOM24 * 0.5
         NEXT(1) = 24
         TWOM12 = TWOM24 * 4096.
      WRITE(6,'(A)') ' FULL INITIALIZATION OF RANLUX WITH 25 INTEGERS:'
      WRITE(6,'(5X,5I12)') ISDEXT
      DO 200 I= 1, 24
      SEEDS(I) = REAL(ISDEXT(I))*TWOM24
  200 CONTINUE
      CARRY = 0.
      IF (ISDEXT(25) .LT. 0)  CARRY = TWOM24
      ISD = IABS(ISDEXT(25))
      I24 = MOD(ISD,100)
      ISD = ISD/100
      J24 = MOD(ISD,100)
      ISD = ISD/100
      IN24 = MOD(ISD,100)
      ISD = ISD/100
      LUXLEV = ISD
        IF (LUXLEV .LE. MAXLEV) THEN
          NSKIP = NDSKIP(LUXLEV)
          WRITE (6,'(A,I2)') ' RANLUX LUXURY LEVEL SET BY RLUXIN TO: ',
     +                         LUXLEV
        ELSE  IF (LUXLEV .GE. 24) THEN
          NSKIP = LUXLEV - 24
          WRITE (6,'(A,I5)') ' RANLUX P-VALUE SET BY RLUXIN TO:',LUXLEV
        ELSE
          NSKIP = NDSKIP(MAXLEV)
          WRITE (6,'(A,I5)') ' RANLUX ILLEGAL LUXURY RLUXIN: ',LUXLEV
          LUXLEV = MAXLEV
        ENDIF
      INSEED = -1
      RETURN
C
C                    Entry to ouput seeds as integers
      ENTRY RLUXUT(ISDEXT)
      DO 300 I= 1, 24
         ISDEXT(I) = INT(SEEDS(I)*TWOP12*TWOP12)
  300 CONTINUE
      ISDEXT(25) = I24 + 100*J24 + 10000*IN24 + 1000000*LUXLEV
      IF (CARRY .GT. 0.)  ISDEXT(25) = -ISDEXT(25)
      RETURN
C
C                    Entry to output the "convenient" restart point
      ENTRY RLUXAT(LOUT,INOUT,K1,K2)
      LOUT = LUXLEV
      INOUT = INSEED
      K1 = KOUNT
      K2 = MKOUNT
      RETURN
C
C                    Entry to initialize from one or three integers
      ENTRY RLUXGO(LUX,INS,K1,K2)
         IF (LUX .LT. 0) THEN
            LUXLEV = LXDFLT
         ELSE IF (LUX .LE. MAXLEV) THEN
            LUXLEV = LUX
         ELSE IF (LUX .LT. 24 .OR. LUX .GT. 2000) THEN
            LUXLEV = MAXLEV
            WRITE (6,'(A,I7)') ' RANLUX ILLEGAL LUXURY RLUXGO: ',LUX
         ELSE
            LUXLEV = LUX
            DO 310 ILX= 0, MAXLEV
              IF (LUX .EQ. NDSKIP(ILX)+24)  LUXLEV = ILX
  310       CONTINUE
         ENDIF
      IF (LUXLEV .LE. MAXLEV)  THEN
         NSKIP = NDSKIP(LUXLEV)
         WRITE(6,'(A,I2,A,I4)') ' RANLUX LUXURY LEVEL SET BY RLUXGO :',
     +        LUXLEV,'     P=', NSKIP+24
      ELSE
          NSKIP = LUXLEV - 24
          WRITE (6,'(A,I5)') ' RANLUX P-VALUE SET BY RLUXGO TO:',LUXLEV
      ENDIF
      IN24 = 0
      IF (INS .LT. 0)  WRITE (6,'(A)')
     +   ' Illegal initialization by RLUXGO, negative input seed'
      IF (INS .GT. 0)  THEN
        JSEED = INS
        WRITE(6,'(A,3I12)') ' RANLUX INITIALIZED BY RLUXGO FROM SEEDS',
     +      JSEED, K1,K2
      ELSE
        JSEED = JSDFLT
        WRITE(6,'(A)')' RANLUX INITIALIZED BY RLUXGO FROM DEFAULT SEED'
      ENDIF
      INSEED = JSEED
      NOTYET = .FALSE.
      TWOM24 = 1.
         DO 325 I= 1, 24
           TWOM24 = TWOM24 * 0.5
         K = JSEED/53668
         JSEED = 40014*(JSEED-K*53668) -K*12211
         IF (JSEED .LT. 0)  JSEED = JSEED+ICONS
         ISEEDS(I) = MOD(JSEED,ITWO24)
  325    CONTINUE
      TWOM12 = TWOM24 * 4096.
         DO 350 I= 1,24
         SEEDS(I) = REAL(ISEEDS(I))*TWOM24
         NEXT(I) = I-1
  350    CONTINUE
      NEXT(1) = 24
      I24 = 24
      J24 = 10
      CARRY = 0.
      IF (SEEDS(24) .EQ. 0.) CARRY = TWOM24
C        If restarting at a break point, skip K1 + IGIGA*K2
C        Note that this is the number of numbers delivered to
C        the user PLUS the number skipped (if luxury .GT. 0).
      KOUNT = K1
      MKOUNT = K2
      IF (K1+K2 .NE. 0)  THEN
        DO 500 IOUTER= 1, K2+1
          INNER = IGIGA
          IF (IOUTER .EQ. K2+1)  INNER = K1
          DO 450 ISK= 1, INNER
            UNI = SEEDS(J24) - SEEDS(I24) - CARRY
            IF (UNI .LT. 0.)  THEN
               UNI = UNI + 1.0
               CARRY = TWOM24
            ELSE
               CARRY = 0.
            ENDIF
            SEEDS(I24) = UNI
            I24 = NEXT(I24)
            J24 = NEXT(J24)
  450     CONTINUE
  500   CONTINUE
C         Get the right value of IN24 by direct calculation
        IN24 = MOD(KOUNT, NSKIP+24)
        IF (MKOUNT .GT. 0)  THEN
           IZIP = MOD(IGIGA, NSKIP+24)
           IZIP2 = MKOUNT*IZIP + IN24
           IN24 = MOD(IZIP2, NSKIP+24)
        ENDIF
C       Now IN24 had better be between zero and 23 inclusive
        IF (IN24 .GT. 23) THEN
           WRITE (6,'(A/A,3I11,A,I5)')
     +    '  Error in RESTARTING with RLUXGO:','  The values', INS,
     +     K1, K2, ' cannot occur at luxury level', LUXLEV
           IN24 = 0
        ENDIF
      ENDIF
      RETURN
      END
      SUBROUTINE SORTR(A,NC,NR,NS)
C
C CERN PROGLIB# M107    SORTR           .VERSION KERNFOR  4.21  890323
C ORIG. 15/11/88 FCA.  Modified from CERNLIB source by S. Yellin
C
C A has NC columns and NR rows.  If NS > 0 the column number NS is put
C in increasing order, and all other row entries, if any, move with the
C row entry of column NS.  If NS < 0 the same thing happens but with decreasing
C order.  The convention for storage of A is all NC columns of row 1, then of
C row 2, ..., row NR, which is the Fortran convention for A(NC,NR).  So if
C NC=NS=1 with a 1-dimensional array A(NR), that array will end up in
C increasing order.
C
      DIMENSION A(NC,NR)
C
      NCS=ABS(NS)
      IF(NCS.EQ.0)  GO TO 999
      IF(NCS.GT.NC) GO TO 999
      IF(NR.LE.1)   GO TO 999
      IF (NS.LE.0)           GO TO 31

C----              Ascending order

      DO 30 J=1,NR-1
        LMIN = J
        HMIN = A(NCS,J)
        DO 10 K=J+1,NR
          IF(HMIN.GT.A(NCS,K)) THEN
            HMIN = A(NCS,K)
            LMIN = K
          ENDIF
   10   CONTINUE
        IF(LMIN.NE.J) THEN
          DO 25 L=LMIN, J, -1
            IF(A(NCS,L).EQ.A(NCS,J)) THEN
              DO 20 K=1,NC
                TEMP      = A(K,LMIN)
                A(K,LMIN) = A(K,L)
                A(K,L)    = TEMP
   20         CONTINUE
              LMIN = L
            ENDIF
   25     CONTINUE
        ENDIF
   30 CONTINUE
      GO TO 999

C----              Descending order

   31 DO 60 J=1,NR-1
        LMAX = J
        HMAX = A(NCS,J)
        DO 40 K=J+1,NR
          IF(HMAX.LT.A(NCS,K)) THEN
            HMAX = A(NCS,K)
            LMAX = K
          ENDIF
   40   CONTINUE
        IF(LMAX.NE.J) THEN
          DO 55 L=LMAX, J, -1
            IF(A(NCS,L).EQ.A(NCS,J)) THEN
              DO 50 K=1,NC
                TEMP      = A(K,LMAX)
                A(K,LMAX) = A(K,L)
                A(K,L)    = TEMP
   50         CONTINUE
              LMAX = L
            ENDIF
   55     CONTINUE
        ENDIF
   60 CONTINUE
  999 RETURN
      END

max_stars_repo_name: slwatkins/DarkLim

max_issues_repo_name: slwatkins/DarkLim

max_forks_repo_name: slwatkins/DarkLim


text:
      SUBROUTINE gratio(a,x,ans,qans,ind)
C ----------------------------------------------------------------------
C        EVALUATION OF THE INCOMPLETE GAMMA RATIO FUNCTIONS
C                      P(A,X) AND Q(A,X)
C
C                        ----------
C
C     IT IS ASSUMED THAT A AND X ARE NONNEGATIVE, WHERE A AND X
C     ARE NOT BOTH 0.
C
C     ANS AND QANS ARE VARIABLES. GRATIO ASSIGNS ANS THE VALUE
C     P(A,X) AND QANS THE VALUE Q(A,X). IND MAY BE ANY INTEGER.
C     IF IND = 0 THEN THE USER IS REQUESTING AS MUCH ACCURACY AS
C     POSSIBLE (UP TO 14 SIGNIFICANT DIGITS). OTHERWISE, IF
C     IND = 1 THEN ACCURACY IS REQUESTED TO WITHIN 1 UNIT OF THE
C     6-TH SIGNIFICANT DIGIT, AND IF IND .NE. 0,1 THEN ACCURACY
C     IS REQUESTED TO WITHIN 1 UNIT OF THE 3RD SIGNIFICANT DIGIT.
C
C     ERROR RETURN ...
C        ANS IS ASSIGNED THE VALUE 2 WHEN A OR X IS NEGATIVE,
C     WHEN A*X = 0, OR WHEN P(A,X) AND Q(A,X) ARE INDETERMINANT.
C     P(A,X) AND Q(A,X) ARE COMPUTATIONALLY INDETERMINANT WHEN
C     X IS EXCEEDINGLY CLOSE TO A AND A IS EXTREMELY LARGE.
C ----------------------------------------------------------------------
C     WRITTEN BY ALFRED H. MORRIS, JR.
C        NAVAL SURFACE WEAPONS CENTER
C        DAHLGREN, VIRGINIA
C     --------------------
C     .. Scalar Arguments ..
      DOUBLE PRECISION a,ans,qans,x
      INTEGER ind
C     ..
C     .. Local Scalars ..
      DOUBLE PRECISION a2n,a2nm1,acc,alog10,am0,amn,an,an0,apn,b2n,
     +                 b2nm1,c,c0,c1,c2,c3,c4,c5,c6,cma,d10,d20,d30,d40,
     +                 d50,d60,d70,e,e0,g,h,j,l,r,rt2pin,rta,rtpi,rtx,s,
     +                 sum,t,t1,third,tol,twoa,u,w,x0,y,z
      INTEGER i,iop,m,max,n
C     ..
C     .. Local Arrays ..
      DOUBLE PRECISION acc0(3),big(3),d0(13),d1(12),d2(10),d3(8),d4(6),
     +                 d5(4),d6(2),e00(3),wk(20),x00(3)
C     ..
C     .. External Functions ..
      DOUBLE PRECISION erf,erfc1,gam1,gamma,rexp,rlog,spmpar
      EXTERNAL erf,erfc1,gam1,gamma,rexp,rlog,spmpar
C     ..
C     .. Intrinsic Functions ..
      INTRINSIC abs,dble,dlog,dmax1,exp,int,sqrt
C     ..
C     .. Data statements ..
C     --------------------
C     --------------------
C     ALOG10 = LN(10)
C     RT2PIN = 1/SQRT(2*PI)
C     RTPI   = SQRT(PI)
C     --------------------
C     --------------------
C     --------------------
C     --------------------
C     --------------------
C     --------------------
C     --------------------
C     --------------------
C     --------------------
      DATA acc0(1)/5.D-15/,acc0(2)/5.D-7/,acc0(3)/5.D-4/
      DATA big(1)/20.0D0/,big(2)/14.0D0/,big(3)/10.0D0/
      DATA e00(1)/.25D-3/,e00(2)/.25D-1/,e00(3)/.14D0/
      DATA x00(1)/31.0D0/,x00(2)/17.0D0/,x00(3)/9.7D0/
      DATA alog10/2.30258509299405D0/
      DATA rt2pin/.398942280401433D0/
      DATA rtpi/1.77245385090552D0/
      DATA third/.333333333333333D0/
      DATA d0(1)/.833333333333333D-01/,d0(2)/-.148148148148148D-01/,
     +     d0(3)/.115740740740741D-02/,d0(4)/.352733686067019D-03/,
     +     d0(5)/-.178755144032922D-03/,d0(6)/.391926317852244D-04/,
     +     d0(7)/-.218544851067999D-05/,d0(8)/-.185406221071516D-05/,
     +     d0(9)/.829671134095309D-06/,d0(10)/-.176659527368261D-06/,
     +     d0(11)/.670785354340150D-08/,d0(12)/.102618097842403D-07/,
     +     d0(13)/-.438203601845335D-08/
      DATA d10/-.185185185185185D-02/,d1(1)/-.347222222222222D-02/,
     +     d1(2)/.264550264550265D-02/,d1(3)/-.990226337448560D-03/,
     +     d1(4)/.205761316872428D-03/,d1(5)/-.401877572016461D-06/,
     +     d1(6)/-.180985503344900D-04/,d1(7)/.764916091608111D-05/,
     +     d1(8)/-.161209008945634D-05/,d1(9)/.464712780280743D-08/,
     +     d1(10)/.137863344691572D-06/,d1(11)/-.575254560351770D-07/,
     +     d1(12)/.119516285997781D-07/
      DATA d20/.413359788359788D-02/,d2(1)/-.268132716049383D-02/,
     +     d2(2)/.771604938271605D-03/,d2(3)/.200938786008230D-05/,
     +     d2(4)/-.107366532263652D-03/,d2(5)/.529234488291201D-04/,
     +     d2(6)/-.127606351886187D-04/,d2(7)/.342357873409614D-07/,
     +     d2(8)/.137219573090629D-05/,d2(9)/-.629899213838006D-06/,
     +     d2(10)/.142806142060642D-06/
      DATA d30/.649434156378601D-03/,d3(1)/.229472093621399D-03/,
     +     d3(2)/-.469189494395256D-03/,d3(3)/.267720632062839D-03/,
     +     d3(4)/-.756180167188398D-04/,d3(5)/-.239650511386730D-06/,
     +     d3(6)/.110826541153473D-04/,d3(7)/-.567495282699160D-05/,
     +     d3(8)/.142309007324359D-05/
      DATA d40/-.861888290916712D-03/,d4(1)/.784039221720067D-03/,
     +     d4(2)/-.299072480303190D-03/,d4(3)/-.146384525788434D-05/,
     +     d4(4)/.664149821546512D-04/,d4(5)/-.396836504717943D-04/,
     +     d4(6)/.113757269706784D-04/
      DATA d50/-.336798553366358D-03/,d5(1)/-.697281375836586D-04/,
     +     d5(2)/.277275324495939D-03/,d5(3)/-.199325705161888D-03/,
     +     d5(4)/.679778047793721D-04/
      DATA d60/.531307936463992D-03/,d6(1)/-.592166437353694D-03/,
     +     d6(2)/.270878209671804D-03/
      DATA d70/.344367606892378D-03/
C     ..
C     .. Executable Statements ..
C     --------------------
C     ****** E IS A MACHINE DEPENDENT CONSTANT. E IS THE SMALLEST
C            FLOATING POINT NUMBER FOR WHICH 1.0 + E .GT. 1.0 .
C
      e = spmpar(1)
C
C     --------------------
      IF (a.LT.0.0D0 .OR. x.LT.0.0D0) GO TO 430
      IF (a.EQ.0.0D0 .AND. x.EQ.0.0D0) GO TO 430
      IF (a*x.EQ.0.0D0) GO TO 420
C
      iop = ind + 1
      IF (iop.NE.1 .AND. iop.NE.2) iop = 3
      acc = dmax1(acc0(iop),e)
      e0 = e00(iop)
      x0 = x00(iop)
C
C            SELECT THE APPROPRIATE ALGORITHM
C
      IF (a.GE.1.0D0) GO TO 10
      IF (a.EQ.0.5D0) GO TO 390
      IF (x.LT.1.1D0) GO TO 160
      t1 = a*dlog(x) - x
      u = a*exp(t1)
      IF (u.EQ.0.0D0) GO TO 380
      r = u* (1.0D0+gam1(a))
      GO TO 250
C
   10 IF (a.GE.big(iop)) GO TO 30
      IF (a.GT.x .OR. x.GE.x0) GO TO 20
      twoa = a + a
      m = int(twoa)
      IF (twoa.NE.dble(m)) GO TO 20
      i = m/2
      IF (a.EQ.dble(i)) GO TO 210
      GO TO 220

   20 t1 = a*dlog(x) - x
      r = exp(t1)/gamma(a)
      GO TO 40
C
   30 l = x/a
      IF (l.EQ.0.0D0) GO TO 370
      s = 0.5D0 + (0.5D0-l)
      z = rlog(l)
      IF (z.GE.700.0D0/a) GO TO 410
      y = a*z
      rta = sqrt(a)
      IF (abs(s).LE.e0/rta) GO TO 330
      IF (abs(s).LE.0.4D0) GO TO 270
C
      t = (1.0D0/a)**2
      t1 = (((0.75D0*t-1.0D0)*t+3.5D0)*t-105.0D0)/ (a*1260.0D0)
      t1 = t1 - y
      r = rt2pin*rta*exp(t1)
C
   40 IF (r.EQ.0.0D0) GO TO 420
      IF (x.LE.dmax1(a,alog10)) GO TO 50
      IF (x.LT.x0) GO TO 250
      GO TO 100
C
C                 TAYLOR SERIES FOR P/R
C
   50 apn = a + 1.0D0
      t = x/apn
      wk(1) = t
      DO 60 n = 2,20
          apn = apn + 1.0D0
          t = t* (x/apn)
          IF (t.LE.1.D-3) GO TO 70
          wk(n) = t
   60 CONTINUE
      n = 20
C
   70 sum = t
      tol = 0.5D0*acc
   80 apn = apn + 1.0D0
      t = t* (x/apn)
      sum = sum + t
      IF (t.GT.tol) GO TO 80
C
      max = n - 1
      DO 90 m = 1,max
          n = n - 1
          sum = sum + wk(n)
   90 CONTINUE
      ans = (r/a)* (1.0D0+sum)
      qans = 0.5D0 + (0.5D0-ans)
      RETURN
C
C                 ASYMPTOTIC EXPANSION
C
  100 amn = a - 1.0D0
      t = amn/x
      wk(1) = t
      DO 110 n = 2,20
          amn = amn - 1.0D0
          t = t* (amn/x)
          IF (abs(t).LE.1.D-3) GO TO 120
          wk(n) = t
  110 CONTINUE
      n = 20
C
  120 sum = t
  130 IF (.NOT.(abs(t).GT.acc)) GO TO 140
      amn = amn - 1.0D0
      t = t* (amn/x)
      sum = sum + t
      GO TO 130
C
  140 max = n - 1
      DO 150 m = 1,max
          n = n - 1
          sum = sum + wk(n)
  150 CONTINUE
      qans = (r/x)* (1.0D0+sum)
      ans = 0.5D0 + (0.5D0-qans)
      RETURN
C
C             TAYLOR SERIES FOR P(A,X)/X**A
C
  160 an = 3.0D0
      c = x
      sum = x/ (a+3.0D0)
      tol = 3.0D0*acc/ (a+1.0D0)
  170 an = an + 1.0D0
      c = -c* (x/an)
      t = c/ (a+an)
      sum = sum + t
      IF (abs(t).GT.tol) GO TO 170
      j = a*x* ((sum/6.0D0-0.5D0/ (a+2.0D0))*x+1.0D0/ (a+1.0D0))
C
      z = a*dlog(x)
      h = gam1(a)
      g = 1.0D0 + h
      IF (x.LT.0.25D0) GO TO 180
      IF (a.LT.x/2.59D0) GO TO 200
      GO TO 190

  180 IF (z.GT.-.13394D0) GO TO 200
C
  190 w = exp(z)
      ans = w*g* (0.5D0+ (0.5D0-j))
      qans = 0.5D0 + (0.5D0-ans)
      RETURN
C
  200 l = rexp(z)
      w = 0.5D0 + (0.5D0+l)
      qans = (w*j-l)*g - h
      IF (qans.LT.0.0D0) GO TO 380
      ans = 0.5D0 + (0.5D0-qans)
      RETURN
C
C             FINITE SUMS FOR Q WHEN A .GE. 1
C                 AND 2*A IS AN INTEGER
C
  210 sum = exp(-x)
      t = sum
      n = 1
      c = 0.0D0
      GO TO 230
C
  220 rtx = sqrt(x)
      sum = erfc1(0,rtx)
      t = exp(-x)/ (rtpi*rtx)
      n = 0
      c = -0.5D0
C
  230 IF (n.EQ.i) GO TO 240
      n = n + 1
      c = c + 1.0D0
      t = (x*t)/c
      sum = sum + t
      GO TO 230

  240 qans = sum
      ans = 0.5D0 + (0.5D0-qans)
      RETURN
C
C              CONTINUED FRACTION EXPANSION
C
  250 tol = dmax1(5.0D0*e,acc)
      a2nm1 = 1.0D0
      a2n = 1.0D0
      b2nm1 = x
      b2n = x + (1.0D0-a)
      c = 1.0D0
  260 a2nm1 = x*a2n + c*a2nm1
      b2nm1 = x*b2n + c*b2nm1
      am0 = a2nm1/b2nm1
      c = c + 1.0D0
      cma = c - a
      a2n = a2nm1 + cma*a2n
      b2n = b2nm1 + cma*b2n
      an0 = a2n/b2n
      IF (abs(an0-am0).GE.tol*an0) GO TO 260
C
      qans = r*an0
      ans = 0.5D0 + (0.5D0-qans)
      RETURN
C
C                GENERAL TEMME EXPANSION
C
  270 IF (abs(s).LE.2.0D0*e .AND. a*e*e.GT.3.28D-3) GO TO 430
      c = exp(-y)
      w = 0.5D0*erfc1(1,sqrt(y))
      u = 1.0D0/a
      z = sqrt(z+z)
      IF (l.LT.1.0D0) z = -z
      IF (iop.lt.2) GO TO 280
      IF (iop.eq.2) GO TO 290
      GO TO 300
C
  280 IF (abs(s).LE.1.D-3) GO TO 340
      c0 = ((((((((((((d0(13)*z+d0(12))*z+d0(11))*z+d0(10))*z+d0(9))*z+
     +     d0(8))*z+d0(7))*z+d0(6))*z+d0(5))*z+d0(4))*z+d0(3))*z+d0(2))*
     +     z+d0(1))*z - third
      c1 = (((((((((((d1(12)*z+d1(11))*z+d1(10))*z+d1(9))*z+d1(8))*z+
     +     d1(7))*z+d1(6))*z+d1(5))*z+d1(4))*z+d1(3))*z+d1(2))*z+d1(1))*
     +     z + d10
      c2 = (((((((((d2(10)*z+d2(9))*z+d2(8))*z+d2(7))*z+d2(6))*z+
     +     d2(5))*z+d2(4))*z+d2(3))*z+d2(2))*z+d2(1))*z + d20
      c3 = (((((((d3(8)*z+d3(7))*z+d3(6))*z+d3(5))*z+d3(4))*z+d3(3))*z+
     +     d3(2))*z+d3(1))*z + d30
      c4 = (((((d4(6)*z+d4(5))*z+d4(4))*z+d4(3))*z+d4(2))*z+d4(1))*z +
     +     d40
      c5 = (((d5(4)*z+d5(3))*z+d5(2))*z+d5(1))*z + d50
      c6 = (d6(2)*z+d6(1))*z + d60
      t = ((((((d70*u+c6)*u+c5)*u+c4)*u+c3)*u+c2)*u+c1)*u + c0
      GO TO 310
C
  290 c0 = (((((d0(6)*z+d0(5))*z+d0(4))*z+d0(3))*z+d0(2))*z+d0(1))*z -
     +     third
      c1 = (((d1(4)*z+d1(3))*z+d1(2))*z+d1(1))*z + d10
      c2 = d2(1)*z + d20
      t = (c2*u+c1)*u + c0
      GO TO 310
C
  300 t = ((d0(3)*z+d0(2))*z+d0(1))*z - third
C
  310 IF (l.LT.1.0D0) GO TO 320
      qans = c* (w+rt2pin*t/rta)
      ans = 0.5D0 + (0.5D0-qans)
      RETURN

  320 ans = c* (w-rt2pin*t/rta)
      qans = 0.5D0 + (0.5D0-ans)
      RETURN
C
C               TEMME EXPANSION FOR L = 1
C
  330 IF (a*e*e.GT.3.28D-3) GO TO 430
      c = 0.5D0 + (0.5D0-y)
      w = (0.5D0-sqrt(y)* (0.5D0+ (0.5D0-y/3.0D0))/rtpi)/c
      u = 1.0D0/a
      z = sqrt(z+z)
      IF (l.LT.1.0D0) z = -z
      IF (iop.lt.2) GO TO 340
      IF (iop.eq.2) GO TO 350
      GO TO 360
C
  340 c0 = ((((((d0(7)*z+d0(6))*z+d0(5))*z+d0(4))*z+d0(3))*z+d0(2))*z+
     +     d0(1))*z - third
      c1 = (((((d1(6)*z+d1(5))*z+d1(4))*z+d1(3))*z+d1(2))*z+d1(1))*z +
     +     d10
      c2 = ((((d2(5)*z+d2(4))*z+d2(3))*z+d2(2))*z+d2(1))*z + d20
      c3 = (((d3(4)*z+d3(3))*z+d3(2))*z+d3(1))*z + d30
      c4 = (d4(2)*z+d4(1))*z + d40
      c5 = (d5(2)*z+d5(1))*z + d50
      c6 = d6(1)*z + d60
      t = ((((((d70*u+c6)*u+c5)*u+c4)*u+c3)*u+c2)*u+c1)*u + c0
      GO TO 310
C
  350 c0 = (d0(2)*z+d0(1))*z - third
      c1 = d1(1)*z + d10
      t = (d20*u+c1)*u + c0
      GO TO 310
C
  360 t = d0(1)*z - third
      GO TO 310
C
C                     SPECIAL CASES
C
  370 ans = 0.0D0
      qans = 1.0D0
      RETURN
C
  380 ans = 1.0D0
      qans = 0.0D0
      RETURN
C
  390 IF (x.GE.0.25D0) GO TO 400
      ans = erf(sqrt(x))
      qans = 0.5D0 + (0.5D0-ans)
      RETURN

  400 qans = erfc1(0,sqrt(x))
      ans = 0.5D0 + (0.5D0-qans)
      RETURN
C
  410 IF (abs(s).LE.2.0D0*e) GO TO 430
  420 IF (x.LE.a) GO TO 370
      GO TO 380
C
C                     ERROR RETURN
C
  430 ans = 2.0D0
      RETURN

      END

max_stars_repo_name: Ennosigaeon/scipy

max_issues_repo_name: Ennosigaeon/scipy

max_forks_repo_name: Ennosigaeon/scipy