text:
module kind_mod
implicit none
private
public :: dp,bad_real,long_int,i32,i64,bad_int
integer, parameter :: dp = selected_real_kind(15, 307)
! integer, parameter :: dp = selected_real_kind(32) ! to get quadruple precision
real(kind=dp), parameter :: bad_real = -999.0_dp
integer, parameter :: i32=selected_int_kind(9), long_int = selected_int_kind(15), i64 = long_int, &
                      bad_int = -999
end module kind_mod

max_stars_repo_name: Beliavsky/RetirementSpending

max_issues_repo_name: Beliavsky/RetirementSpending

max_forks_repo_name: Beliavsky/RetirementSpending


!
!     * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
!     *                                                               *
!     *                  copyright (c) 1998 by UCAR                   *
!     *                                                               *
!     *       University Corporation for Atmospheric Research         *
!     *                                                               *
!     *                      all rights reserved                      *
!     *                                                               *
!     *                          Spherepack                           *
!     *                                                               *
!     *       A Package of Fortran Subroutines and Programs           *
!     *                                                               *
!     *              for Modeling Geophysical Processes               *
!     *                                                               *
!     *                             by                                *
!     *                                                               *
!     *                  John Adams and Paul Swarztrauber             *
!     *                                                               *
!     *                             of                                *
!     *                                                               *
!     *         the National Center for Atmospheric Research          *
!     *                                                               *
!     *                Boulder, Colorado  (80307)  U.S.A.             *
!     *                                                               *
!     *                   which is sponsored by                       *
!     *                                                               *
!     *              the National Science Foundation                  *
!     *                                                               *
!     * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
!
submodule(scalar_laplacian_routines) invert_scalar_laplacian_regular_grid

contains
    !     subroutine islapec(nlat, nlon, isym, nt, xlmbda, sf, ids, jds, a, b, &
    !                        mdab, ndab, wshsec, pertrb, ierror)
    !
    !     islapec inverts the laplace or helmholz operator on an equally
    !     spaced latitudinal grid using o(n**2) storage. given the
    !     spherical harmonic coefficients a(m, n) and b(m, n) of the right
    !     hand side slap(i, j), islapec computes a solution sf(i, j) to
    !     the following helmhotz equation :
    !
    !           2                2
    !     [d(sf(i, j))/dlambda /sint + d(sint*d(sf(i, j))/dtheta)/dtheta]/sint
    !
    !                   - xlmbda * sf(i, j) = slap(i, j)
    !
    !      where sf(i, j) is computed at colatitude
    !
    !                 theta(i) = (i-1)*pi/(nlat-1)
    !
    !            and longitude
    !
    !                 lambda(j) = (j-1)*2*pi/nlon
    !
    !            for i=1, ..., nlat and j=1, ..., nlon.
    !
    !
    !     input parameters
    !
    !     nlat   the number of colatitudes on the full sphere including the
    !            poles. for example, nlat = 37 for a five degree grid.
    !            nlat determines the grid increment in colatitude as
    !            pi/(nlat-1).  if nlat is odd the equator is located at
    !            grid point i=(nlat + 1)/2. if nlat is even the equator is
    !            located half way between points i=nlat/2 and i=nlat/2+1.
    !            nlat must be at least 3. note: on the half sphere, the
    !            number of grid points in the colatitudinal direction is
    !            nlat/2 if nlat is even or (nlat + 1)/2 if nlat is odd.
    !
    !     nlon   the number of distinct longitude points.  nlon determines
    !            the grid increment in longitude as 2*pi/nlon. for example
    !            nlon = 72 for a five degree grid. nlon must be greater
    !            than zero. the axisymmetric case corresponds to nlon=1.
    !            the efficiency of the computation is improved when nlon
    !            is a product of small prime numbers.
    !
    !     isym   this parameter should have the same value input to subroutine
    !            shaec to compute the coefficients a and b for the scalar field
    !            slap.  isym is set as follows:
    !
    !            = 0  no symmetries exist in slap about the equator. scalar
    !                 synthesis is used to compute sf on the entire sphere.
    !                 i.e., in the array sf(i, j) for i=1, ..., nlat and
    !                 j=1, ..., nlon.
    !
    !           = 1  sf and slap are antisymmetric about the equator. the
    !                synthesis used to compute sf is performed on the
    !                northern hemisphere only.  if nlat is odd, sf(i, j) is
    !                computed for i=1, ..., (nlat + 1)/2 and j=1, ..., nlon.  if
    !                nlat is even, sf(i, j) is computed for i=1, ..., nlat/2
    !                and j=1, ..., nlon.
    !
    !
    !           = 2  sf and slap are symmetric about the equator. the
    !                synthesis used to compute sf is performed on the
    !                northern hemisphere only.  if nlat is odd, sf(i, j) is
    !                computed for i=1, ..., (nlat + 1)/2 and j=1, ..., nlon.  if
    !                nlat is even, sf(i, j) is computed for i=1, ..., nlat/2
    !                and j=1, ..., nlon.
    !
    !
    !   nt       the number of solutions. in the program that calls islapec
    !            the arrays sf, a, and b can be three dimensional in which
    !            case multiple solutions are computed. the third index
    !            is the solution index with values k=1, ..., nt.
    !            for a single solution set nt=1. the description of the
    !            remaining parameters is simplified by assuming that nt=1
    !            and sf, a, b are two dimensional.
    !
    !   xlmbda   a one dimensional array with nt elements. if xlmbda is
    !            is identically zero islapec solves poisson's equation.
    !            if xlmbda > 0.0 islapec solves the helmholtz equation.
    !            if xlmbda < 0.0 the nonfatal error flag ierror=-1 is
    !            returned. negative xlambda could result in a division
    !            by zero.
    !
    !   ids      the first dimension of the array sf as it appears in the
    !            program that calls islapec.  if isym = 0 then ids must be at
    !            least nlat.  if isym > 0 and nlat is even then ids must be
    !            at least nlat/2. if isym > 0 and nlat is odd then ids must
    !            be at least (nlat + 1)/2.
    !
    !   jds      the second dimension of the array sf as it appears in the
    !            program that calls islapec. jds must be at least nlon.
    !
    !
    !   a, b      two or three dimensional arrays (see input parameter nt)
    !            that contain scalar spherical harmonic coefficients
    !            of the scalar field slap. a, b must be computed by shaec
    !            prior to calling islapec.
    !
    !
    !   mdab     the first dimension of the arrays a and b as it appears
    !            in the program that calls islapec.  mdab must be at
    !            least min(nlat, (nlon+2)/2) if nlon is even or at least
    !            min(nlat, (nlon + 1)/2) if nlon is odd.
    !
    !   ndab     the second dimension of the arrays a and b as it appears
    !            in the program that calls islapec. ndab must be at least
    !            least nlat.
    !
    !            mdab, ndab should have the same values input to shaec to
    !            compute the coefficients a and b.
    !
    !
    !   wshsec   an array which must be initialized by subroutine shseci.
    !            once initialized, wshsec can be used repeatedly by
    !            islapec as long as nlat and nlon  remain unchanged.
    !            wshsec must not be altered between calls of islapec.
    !
    !   lshsec   the dimension of the array wshsec as it appears in the
    !            program that calls islapec.  let
    !
    !               l1 = min(nlat, (nlon+2)/2) if nlon is even or
    !               l1 = min(nlat, (nlon + 1)/2) if nlon is odd
    !
    !            and
    !
    !               l2 = nlat/2        if nlat is even or
    !               l2 = (nlat + 1)/2    if nlat is odd
    !
    !            then lshsec must be greater than or equal to
    !
    !               2*nlat*l2+3*((l1-2)*(2*nlat-l1-1))/2+nlon+15
    !
    !     output parameters
    !
    !
    !    sf      two or three dimensional arrays (see input parameter nt)
    !            that contain the solution to either the helmholtz
    !            (xlmbda>0.0) or poisson's equation. sf(i, j) is computed
    !            at colatitude
    !
    !                 theta(i) = (i-1)*pi/(nlat-1)
    !
    !            and longitude
    !
    !                 lambda(j) = (j-1)*2*pi/nlon
    !
    !            for i=1, ..., nlat and j=1, ..., nlon.
    !
    !   pertrb  a one dimensional array with nt elements (see input
    !           parameter nt). in the discription that follows we assume
    !           that nt=1. if xlmbda > 0.0 then pertrb=0.0 is always
    !           returned because the helmholtz operator is invertible.
    !           if xlmbda = 0.0 then a solution exists only if a(1, 1)
    !           is zero. islapec sets a(1, 1) to zero. the resulting
    !           solution sf(i, j) solves poisson's equation with
    !           pertrb = a(1, 1)/(2.*sqrt(2.)) subtracted from the
    !           right side slap(i, j).
    !
    !
    !  ierror   a parameter which flags errors in input parameters as follows:
    !
    !            =-1  xlmbda is input negative (nonfatal error)
    !            = 0  no errors detected
    !            = 1  error in the specification of nlat
    !            = 2  error in the specification of nlon
    !            = 3  error in the specification of isym
    !            = 4  error in the specification of nt
    !            = 5  error in the specification of ids
    !            = 6  error in the specification of jds
    !            = 7  error in the specification of mdab
    !            = 8  error in the specification of ndab
    !            = 9  error in the specification of lshsec
    !
    module subroutine islapec(nlat, nlon, isym, nt, xlmbda, sf, ids, jds, &
        a, b, mdab, ndab, wshsec, pertrb, ierror)

        ! Dummy arguments
        integer(ip), intent(in)  :: nlat
        integer(ip), intent(in)  :: nlon
        integer(ip), intent(in)  :: isym
        integer(ip), intent(in)  :: nt
        real(wp),    intent(in)  :: xlmbda(:)
        real(wp),    intent(out) :: sf(ids, jds, nt)
        integer(ip), intent(in)  :: ids
        integer(ip), intent(in)  :: jds
        real(wp),    intent(in)  :: a(mdab, ndab, nt)
        real(wp),    intent(in)  :: b(mdab, ndab, nt)
        integer(ip), intent(in)  :: mdab
        integer(ip), intent(in)  :: ndab
        real(wp),    intent(in)  :: wshsec(:)
        real(wp),    intent(out) :: pertrb(:)
        integer(ip), intent(out) :: ierror

        ! Local variables
        integer(ip)                  :: required_wavetable_size
        type(ScalarSynthesisUtility) :: util

        ! Check input arguments
        required_wavetable_size = util%get_lshsec(nlat, nlon)

        call util%check_scalar_transform_inputs(isym, ids, jds, &
            mdab, ndab, nlat, nlon, nt, required_wavetable_size, &
            wshsec, ierror)

        ! Check sign of xlmbda
        if (any(xlmbda < ZERO)) then
            ierror = -1
        end if

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

        call invert_scalar_laplacian_lower_utility_routine(nlat, nlon, isym, &
            nt, xlmbda, pertrb, sf, a, b, wshsec, shsec, ierror)

    end subroutine islapec

end submodule invert_scalar_laplacian_regular_grid


text:
      real*8 function tlum(np,x,y,xm,ym,sigx,sigy,nmesh)
      use tfstk
      implicit none
      integer*8 ibin
      integer*4 np,nmesh
      real*8 x(2,np),y(2,np),xm,ym,tluma,sigx,sigy
      nmesh=max(int(sqrt(real(np))),6)
      ibin=ktaloc(2*(nmesh+1)**2)
      if(ibin .le. 0)then
        write(*,*)'Insufficient memory in tlum, mesh size =',nmesh
        tlum=0.d0
        return
      endif
      tlum=tluma(np,x,y,xm,ym,sigx,sigy,rlist(ibin),nmesh)
      call tfree(ibin)
      return
      end

ext: f

max_stars_repo_name: noboruatkek/SAD

max_issues_repo_name: noboruatkek/SAD

max_forks_repo_name: noboruatkek/SAD


text:
module da_gen_be

   !----------------------------------------------------------------------------
   ! Purpose: Collection of routines required by gen_be BE stats calculation
   ! code.
   !
   !  Update: Multivariate BE option (cv_options=6)
   !          Syed RH Rizvi (MMM/NESL/NCAR)   Date: 02/01/2010
   !
   !  Note: Please acknowledge author/institute in work that uses this code.
   !------------------------------------------------------------------------

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

   use da_control, only : stdout,vertical_ip, t0,es_beta,es_alpha, &
      es_gamma,kappa,rd_over_rv,rd_over_rv1,t_kelvin, gravity, &
      filename_len,vertical_ip_0, trace_use, trace_use_dull, cv_options, use_rf, do_normalize
   use da_reporting, only : da_error, da_warning, da_message, message
   use da_tools_serial, only : da_get_unit, da_free_unit, da_array_print
   use da_lapack, only : dsyev
   use da_wavelet, only: lf,namw,nb,nij,ws,wsd

   implicit none

   real, parameter :: base_pres = 100000.0 ! Hardwired - link to namelist later.

contains

! Stubs to avoid picking up vast ammounts of unnecessary stuff from wrfvar

subroutine da_trace_entry(name, message, messages, maxnocalls)
   implicit none
   character (len=*),           intent(in) :: name
   character (len=*), optional, intent(in) :: message
   character (len=*), optional, intent(in) :: messages(:)
   integer, optional,           intent(in) :: maxnocalls
end subroutine da_trace_entry

subroutine da_trace_exit(name, message, messages, maxnocalls)
   implicit none
   character (len=*), intent(in)           :: name
   character (len=*), optional, intent(in) :: message
   character (len=*), optional, intent(in) :: messages(:)
   integer, optional, intent(in)           :: maxnocalls
end subroutine da_trace_exit

#include "da_create_bins.inc"
#include "da_filter_regcoeffs.inc"
#include "da_get_field.inc"
#include "da_get_height.inc"
#include "da_get_trh.inc"
#include "da_print_be_stats_h_global.inc"
#include "da_print_be_stats_h_regional.inc"
#include "da_print_be_stats_p.inc"
#include "da_print_be_stats_v.inc"
#include "da_readwrite_be_stage1.inc"
#include "da_readwrite_be_stage2.inc"
#include "da_readwrite_be_stage3.inc"
#include "da_readwrite_be_stage4.inc"
#include "da_stage0_initialize.inc"

   ! Files from other modules:
#include "da_transform_vptovv.inc"
#include "da_eof_decomposition.inc"
#include "da_eof_decomposition_test.inc"
#include "da_perform_2drf.inc"
#include "da_recursive_filter_1d.inc"

end module da_gen_be

subroutine wrf_abort
   stop
end subroutine wrf_abort

   LOGICAL FUNCTION wrf_dm_on_monitor()
      IMPLICIT NONE
      wrf_dm_on_monitor = .TRUE.
      RETURN
   END FUNCTION wrf_dm_on_monitor

max_stars_repo_name: matzegoebel/WRF-fluxavg

max_issues_repo_name: teb-model/wrf-teb

max_forks_repo_name: teb-model/wrf-teb


text:
      SUBROUTINE getiota(phipog, bsupu, bsupv)
      USE vmec_main
      USE realspace
      IMPLICIT NONE
C-----------------------------------------------
C   D u m m y   A r g u m e n t s
C-----------------------------------------------
      REAL(rprec), DIMENSION(nrzt), INTENT(in) :: phipog, bsupv
      REAL(rprec), DIMENSION(nrzt), INTENT(inout) :: bsupu
C-----------------------------------------------
C   L o c a l   P a r a m e t e r s
C-----------------------------------------------
      REAL(rprec), PARAMETER :: p5=0.5_dp, c1p5=1.5_dp
C-----------------------------------------------
C   L o c a l   V a r i a b l e s
C-----------------------------------------------
      INTEGER :: js, l
      REAL(rprec) :: top, bot
C-----------------------------------------------

      IF (ncurr .eq. 0) RETURN
      DO js = 2, ns
         top = icurv(js)                                !!Integrated toroidal current
         bot = 0
         DO l = js, nrzt, ns
            top = top - wint(l)*(guu(l)*bsupu(l) + guv(l)*bsupv(l))
            bot = bot + wint(l)*phipog(l)*guu(l)
         END DO
         iotas(js) = top/bot
      END DO

!     Do not compute iota too near origin
      iotaf(1)  = c1p5*iotas(2) - p5*iotas(3)           !!zero gradient near axis
!      iotaf(ns) = c1p5*iotas(ns) - p5*iotas(ns-1)
!      DO js = 2, ns-1
!         iotaf(js) = p5*(iotas(js) + iotas(js+1))
!      END DO
      DO js = 2, ns
         iotaf(js) = 2*iotas(js) - iotaf(js-1)
      END DO

      END SUBROUTINE getiota


! ====================================================================
!     Subroutine module
!
!     Purpose: store the variables that will be used in the model
!              frequently.
! ====================================================================

    MODULE parm
    IMPLICIT NONE

    INTEGER, PUBLIC, PARAMETER :: KR = SELECTED_REAL_KIND(r=20, p=20) !Accuracy of Real Number.
    INTEGER, PUBLIC, PARAMETER :: KI = SELECTED_INT_KIND(4)           !Accuracy of Int Number.
    INTEGER, PUBLIC, PARAMETER :: MMPLD=50          ! Max print time and boundary condition time.
    INTEGER, PUBLIC, PARAMETER :: NMATD=8           ! Max numbers of soil materials (Vertical)
    INTEGER, PUBLIC, PARAMETER :: NREGD=8           ! Max numbers of sub-regions.
    INTEGER, PUBLIC, PARAMETER :: NPARD=10          ! Max numbers of hydraulic parameters.
    INTEGER, PUBLIC, PARAMETER :: NOBSD=10          ! Max numbers of observation nodes.
    INTEGER, PUBLIC, PARAMETER :: NLAYERD=130       ! Max numbers of soil layers.
    INTEGER, PUBLIC, PARAMETER :: numc=1            ! Max numbers of plant species.

    LOGICAL (KIND=4) :: lCheck         ! If check the input data and write the check.out file.
    LOGICAL (KIND=4) :: lWat           ! True, transient water flow; False, steady water flow.
    LOGICAL (KIND=4) :: lChem          ! If calculate the solute transport.
    LOGICAL (KIND=4) :: lTemp          ! If calculate the temperature.
    LOGICAL (KIND=4) :: lCrop          ! If calculate the crop root uptake.
    LOGICAL (KIND=4) :: IfPM           ! True, calculate the ET0 with P-M equation; False, direct input.

    INTEGER (KIND=KI) :: Terr          ! Error check index.
    INTEGER (KIND=KI) :: Nlayer        ! Numbers of real soil layers.
    INTEGER (KIND=KI) :: Wlayer
    INTEGER (KIND=KI) :: NObs          ! Numbers of the observation points.
    INTEGER (KIND=KI) :: Nmat          ! Numbers of soil materials.
    INTEGER (KIND=KI) :: NReg          ! Numbers of the sub-region, for mass balance calculate.
    INTEGER (KIND=KI) :: MPL           ! Numbers of print times.
    INTEGER (KIND=KI) :: Bup           ! Upper boundary condition. 1, atmospheric boundary, open 01.wea;
                                       ! 0, given flux (irrigation), open met.in;
                                       ! <0, precipitation and irrigation, open both files.
    INTEGER (KIND=KI) :: Bdn           ! Lower boundary condition. 0, 0 flux boundary; -1, given water content,
                                       ! open met.in; -2, given groundwater depth, open met in.
    INTEGER (KIND=KI) :: Npar          ! Numbers of soil hydraulic parameters.
    INTEGER (KIND=KI) :: Drng          ! Which Drainage Function 1/2/3/4/5/6/7.
    INTEGER (kind=KI) :: Dfit          ! The Empirical Formula of Diffusion -3/-2/-1/1/2/3.
    INTEGER (KIND=KI) :: ddn           ! Diffision process related.
    INTEGER (KIND=KI) :: Nup           ! Numbers of flux upper boundary node.
    INTEGER (KIND=KI) :: Ndn           ! Numbers of flux lower boundary node.
    INTEGER (KIND=KI) :: MaxAL         ! Numbers of atmospheric data-records.
    INTEGER (KIND=KI) :: Plevel        ! Print time counter.
    INTEGER (KIND=KI) :: Tlevel        ! Calculate time counter.
    INTEGER (KIND=KI) :: Alevel        ! Atmosphere time counter.
    INTEGER (KIND=KI), ALLOCATABLE :: Obs(:)    ! The observation points.
    INTEGER (KIND=KI), ALLOCATABLE :: MATuz(:)  ! The material serial number of each layer.
    INTEGER (KIND=KI), ALLOCATABLE :: REGuz(:)  ! The sub-region serial numer of each layer.

    CHARACTER (LEN=8) :: date          ! The date of the calculation.
    CHARACTER (LEN=5) :: LUnit         ! Length Unit.
    CHARACTER (LEN=5) :: TUnit         ! Time Unit.
    CHARACTER (LEN=5) :: MUnit         ! Mass Unit.

    REAL (KIND=KR) :: Tol=1E-10_KR     ! The tolerance while calculation.
    REAL (KIND=KR) :: CosAlf           ! The Cosine value of the angle.
    REAL (KIND=KR) :: dt               ! Time step.
    REAL (KIND=KR) :: dtOld            ! The initial input time step.
    REAL (KIND=KR) :: t                ! Time.
    REAL (KIND=KR) :: tinit            ! The initial time.
    REAL (KIND=KR) :: tEnd             ! The endding time.
    REAL (KIND=KR) :: tAtm             ! The time for the atmospheric condition.
    REAL (KIND=KR) :: qair             ! Precipitation.
    REAL (KIND=KR) :: qirr             ! Irrigation.
    REAL (KIND=KR) :: TotalPerco       ! Percolation flux.
    REAL (KIND=KR) :: qairt            ! Sum of Precipitation.
    REAL (KIND=KR) :: qirrt            ! Sum of Irrigation.
    REAL (KIND=KR) :: qbtmt            ! Sum of Bottom Flux.
    REAL (KIND=KR) :: CumE             ! Sum of Evaporation.
    REAL (KIND=KR) :: CumT             ! Sum of Transpiration.
    REAL (KIND=KR) :: sinkt            ! Sum of source/sink.
    REAL (KIND=KR) :: voluz            ! The initial water volumn.
    REAL (KIND=KR) :: voluz1           ! The new water volumn.
    REAL (KIND=KR) :: Epa	           ! Actual Evaporation.
    REAL (KIND=KR) :: Tra	           ! Actual Transpiration.
    REAL (KIND=KR) :: ptab
    REAL (KIND=KR) :: xConv            ! Length conversion coefficient.
    REAL (KIND=KR) :: tConv            ! Time conversion coefficient.
    REAL (KIND=KR) :: mConv            ! Mass conversion coefficient.
    REAL (KIND=KR) :: rTop             ! If lWat=false, upper flux
    REAL (KIND=KR) :: rBot             ! If lWat=false, lower flux
    REAL (KIND=KR) :: rET            ! If lWat=false, potential evapotranspiration

    REAL (KIND=KR), ALLOCATABLE :: TPrint(:)        ! The print time.
    REAL (KIND=KR), ALLOCATABLE :: ths(:)           ! ths for each materials.
    REAL (KIND=KR), ALLOCATABLE :: thF(:)           ! Field capacity.
    REAL (KIND=KR), ALLOCATABLE :: thw(:)           ! Wilting point.
    REAL (KIND=KR), ALLOCATABLE :: up(:,:)          ! Met.in file, Upper input.
    REAL (KIND=KR), ALLOCATABLE :: dn(:,:)          ! Met.in file, Lower input.
    REAL (KIND=KR), ALLOCATABLE :: sp(:,:)          ! The Evaporation cumulative distribution function.
    REAL (KIND=KR), ALLOCATABLE :: Par(:,:)         ! The Hydraulic parameters.
    REAL (KINd=KR), ALLOCATABLE :: evatra(:,:),precip(:,:) ! The input value.
    ! Nlayer dimension array
    REAL (KIND=KR), ALLOCATABLE :: dz(:)            ! Thickness of Each Layer.
    REAL (KIND=KR), ALLOCATABLE :: zx(:)            ! The z Coordinates of Every Node Descend from Up to Down.
    REAL (KIND=KR), ALLOCATABLE :: th(:)            ! Soil Water Content.
    REAL (KIND=KR), ALLOCATABLE :: SInk1d(:)        ! The Actural Evapotranspiration for Each Layer.
    REAL (KIND=KR), ALLOCATABLE :: Epi(:)           ! Evaporation for Each Layer.
    REAL (KIND=KR), ALLOCATABLE :: Tri(:)           ! Transpiration for Each Layer.
    REAL (KIND=KR), ALLOCATABLE :: Slope(:)         ! Empirical Formula of Diffusion.
    REAL (KIND=KR), ALLOCATABLE :: Intercept(:)     ! Empirical Formula of Diffusion.
    REAL (KIND=KR), ALLOCATABLE :: par_n(:)         ! Empirical Formula of n.

!   Solute transport
    LOGICAL (KIND=4) :: MIM             ! If mobile-immobile.
    LOGICAL (KIND=4) :: DiCr            ! If Dissolution and Crystallization process.

    REAL (KIND=KR) :: ct
    REAL (KIND=KR) :: mvoluz            ! The mobile water volumn.
    REAL (KIND=KR) :: imvoluz           ! The immobile water volumn.
    REAL (KIND=KR) :: gdep              ! Groundwater level (L).
    REAL (KIND=KR) :: dth               ! Given th at lower boundary.

    REAL (KIND=KR), ALLOCATABLE :: qflux(:,:) ! The water flux of each layer.
    REAL (KIND=KR), ALLOCATABLE :: thini(:)   ! The initial soil water content.
    REAL (KIND=KR), ALLOCATABLE :: tht(:)     !
    REAL (KIND=KR), ALLOCATABLE :: thm(:)     ! Mobile soil water content.
    REAL (KIND=KR), ALLOCATABLE :: thim(:)    ! Immobile Soil Water Content.
    REAL (KIND=KR), ALLOCATABLE :: thmini(:)  ! The initial mobile soil water content.
    REAL (KIND=KR), ALLOCATABLE :: thimini(:) ! The initial immobile soil water content.

    INTEGER (KIND=KI) :: CBup      ! Solute up boundary condition.
    INTEGER (KIND=KI) :: CBdn      ! Solute lower boundary condition.
    INTEGER (KIND=KI) :: CNup      ! Numbers of Solute Upper Boundary Node.
    INTEGER (KIND=KI) :: CNdn      ! Numbers of Solute Lower Boundary Node.
    INTEGER (KIND=KI) :: Cind      ! The factor of solute calculation time step.
    INTEGER (KIND=KI) :: index

    REAL (KIND=KR) :: w            ! The weight in solute convection equation.
    REAL (KIND=KR) :: Concup       ! The concentration of the upper boundary.
    REAL (KIND=KR) :: Concdn       ! The Concentration of the lower boundary.
    REAL (KIND=KR) :: Svoluz       ! The solute volumn.
    REAL (KIND=KR) :: Smvoluz      ! The mobile region solute volumn.
    REAL (KIND=KR) :: Simvoluz     ! The immobile region solute volumn.
    REAL (KIND=KR) :: Sup          ! The solution comes from upper boundary.
    REAL (KIND=KR) :: Sdn          ! The solution comes from lower boundary.
    REAL (KIND=KR) :: SSinkt       ! The solution source/sink term.
    REAL (KIND=KR) :: Csat         ! Solubility saturation

    REAL (KIND=KR), ALLOCATABLE :: rou(:)       ! density, M3/L.
    REAL (KIND=KR), ALLOCATABLE :: kx(:)        ! Adsorption Coefficient, L3/M-1.
    REAL (KIND=KR), ALLOCATABLE :: miuw(:)      ! Water First Order Reaction Rate Coefficient, T-1.
    REAL (KIND=KR), ALLOCATABLE :: mius(:)      ! Solid First Order Reaction Rate Coefficient, T-1.
    REAL (KIND=KR), ALLOCATABLE :: gamaw(:)     ! Water Zero Order Reaction Rate Coefficient, ML-3T-1.
    REAL (KIND=KR), ALLOCATABLE :: gamas(:)     ! Solid Zero Order Reaction Rate Coefficient, T-1.
    REAL (KIND=KR), ALLOCATABLE :: tao(:)       ! Dissolution Rate of Crystalline Salt, T-1.
    REAL (KIND=KR), ALLOCATABLE :: Cup(:,:)     ! Solute Upper Boundary.
    REAL (KIND=KR), ALLOCATABLE :: Cdn(:,:)     ! Solute Lower Boundary.
    REAL (KIND=KR), ALLOCATABLE :: ChPar(:,:)   ! The Hydraulic parameters.
    REAL (KIND=KR), ALLOCATABLE :: Conc(:)      ! Solute Concentration. MIM True, Mobile Solute Concentration.
    REAL (KIND=KR), ALLOCATABLE :: Concini(:)   ! The Initial Solute Concentration.
    REAL (KIND=KR), ALLOCATABLE :: Conc1(:)
    REAL (KIND=KR), ALLOCATABLE :: Conc2(:)
    REAL (KIND=KR), ALLOCATABLE :: Conc3(:)
    REAL (KIND=KR), ALLOCATABLE :: Conc4(:)
    REAL (KIND=KR), ALLOCATABLE :: Conim(:)     ! Immobile Solute Concentration.
    REAL (KIND=KR), ALLOCATABLE :: Conimini(:)  ! The Initial Immobile Solute Concentration.
    REAL (KIND=KR), ALLOCATABLE :: Ssalt(:)     ! The Initial Mobile Soil Water Content.
    REAL (KIND=KR), ALLOCATABLE :: Ssalim(:)    ! The Initial Immobile Soil Water Content.
    REAL (KIND=KR), ALLOCATABLE :: Ssink1d(:)
    REAL (KIND=KR), ALLOCATABLE :: Ssink1dm(:)
    REAL (KIND=KR), ALLOCATABLE :: Ssink1dim(:)

    END MODULE parm


text:
!!!!!! sub.f90 !!!!!!
! This program has some of the subroutines.
! subroutine list is written bellow:
!
! addab(a, b)
! addab2c(a, b, c)
! calc(n, Dict, puzzle, cover, mode)
!

!!! addab !!!
subroutine addab(a, b)
  integer a, b
  a = a + b
  b = a + b
  return
end subroutine addab

!!! addab2c !!!
subroutine addab2c(a, b, c)
  integer a, b, c
  c = a + b
  return
end subroutine addab2c

max_stars_repo_name: sushi-kimura/Crossword-LocalSearch

max_issues_repo_name: sushi-kimura/Crossword-LocalSearch

max_forks_repo_name: sushi-kimura/Crossword-LocalSearch


text:
! reference: R. J. Purser (J. of Clim. and Apld. Meteorology, 1987)

module filter_module

  ! 3-point filter: T=1, R=S=0
  use amrex_fort_module, only : rt => amrex_real
  real(rt)        , dimension(0:1), parameter :: ff1 = &
       (/ 0.5e0_rt, 0.25e0_rt /)
  ! For boundary cell
  real(rt)        , dimension(0:1), parameter :: ff1b = &
       (/ 0.75e0_rt, 0.25e0_rt /)

  ! 5-point filter: T=2, R+S+1=T
  ! S = 0, 1
  real(rt)        , dimension(0:2,0:1), parameter :: ff2 = reshape( &
       [ 0.625e0_rt, 0.25e0_rt, -0.0625e0_rt, &
       & 0.375e0_rt, 0.25e0_rt,  0.0625e0_rt ], &
       [3,2] )
  ! For first bondary cell
  real(rt)        , dimension(0:2), parameter :: ff2b0 = &
       (/ 17.e0_rt/16.e0_rt, -2.e0_rt/16.e0_rt, 1.e0_rt/16.e0_rt /)
  ! for second boundary cell
  real(rt)        , dimension(-1:2), parameter :: ff2b1 = &
       (/ -2.e0_rt/16.e0_rt, 21.e0_rt/16.e0_rt, -4.e0_rt/16.e0_rt, 1.e0_rt/16.e0_rt /)

  ! 7-point filter: T=3, R+S+1=T
  ! S = 0, 1, 2
  real(rt)        , dimension(0:3,0:2), parameter :: ff3 = reshape( &
       [ 44.e0_rt/64.e0_rt, 15.e0_rt/64.e0_rt, -6.e0_rt/64.e0_rt,  1.e0_rt/64.e0_rt, &
       & 32.e0_rt/64.e0_rt, 18.e0_rt/64.e0_rt,  0.e0_rt      , -2.e0_rt/64.e0_rt, &
       & 20.e0_rt/64.e0_rt, 15.e0_rt/64.e0_rt,  6.e0_rt/64.e0_rt,  1.e0_rt/64.e0_rt ], &
       [4,3] )
  ! for boundary cells
  real(rt)        , dimension(0:3), parameter :: ff3b0 = &
       (/ 63.e0_rt/64.e0_rt, 3.e0_rt/64.e0_rt, -3.e0_rt/64.e0_rt, 1.e0_rt/64.e0_rt /)
  real(rt)        , dimension(-1:3), parameter :: ff3b1 = &
       (/ 3.e0_rt/64.e0_rt, &
       54.e0_rt/64.e0_rt, 12.e0_rt/64.e0_rt, -6.e0_rt/64.e0_rt, 1.e0_rt/64.e0_rt /)
  real(rt)        , dimension(-2:3), parameter :: ff3b2 = &
       (/ -3.e0_rt/64.e0_rt, 12.e0_rt/64.e0_rt, &
       45.e0_rt/64.e0_rt, 15.e0_rt/64.e0_rt, -6.e0_rt/64.e0_rt, 1.e0_rt/64.e0_rt /)

  ! 9-point filter: T=4, R+S+1=T
  ! S = 0, 1, 2, 3
  real(rt)        , dimension(0:4,0:3), parameter :: ff4 = reshape( &
       [ 186.e0_rt/256.e0_rt, 56.e0_rt/256.e0_rt, -28.e0_rt/256.e0_rt,  8.e0_rt/256.e0_rt, -1.e0_rt/256.e0_rt, &
       & 146.e0_rt/256.e0_rt, 72.e0_rt/256.e0_rt, -12.e0_rt/256.e0_rt, -8.e0_rt/256.e0_rt,  3.e0_rt/256.e0_rt, &
       & 110.e0_rt/256.e0_rt, 72.e0_rt/256.e0_rt,  12.e0_rt/256.e0_rt, -8.e0_rt/256.e0_rt, -3.e0_rt/256.e0_rt, &
       &  70.e0_rt/256.e0_rt, 56.e0_rt/256.e0_rt,  28.e0_rt/256.e0_rt,  8.e0_rt/256.e0_rt,  1.e0_rt/256.e0_rt ], &
       [5,4] )
  ! for boundary cells
  real(rt)        , dimension(0:4), parameter :: ff4b0 = &
       (/ 257.e0_rt/256.e0_rt, -4.e0_rt/256.e0_rt, 6.e0_rt/256.e0_rt, -4.e0_rt/256.e0_rt, 1.e0_rt/256.e0_rt /)
  real(rt)        , dimension(-1:4), parameter :: ff4b1 = &
       (/ -4.e0_rt/256.e0_rt,  &
       273.e0_rt/256.e0_rt, -28.e0_rt/256.e0_rt, 22.e0_rt/256.e0_rt, -8.e0_rt/256.e0_rt, 1.e0_rt/256.e0_rt /)
  real(rt)        , dimension(-2:4), parameter :: ff4b2 = &
       (/ 6.e0_rt/256.e0_rt, -28.e0_rt/256.e0_rt, &
       309.e0_rt/256.e0_rt, -52.e0_rt/256.e0_rt, 28.e0_rt/256.e0_rt, -8.e0_rt/256.e0_rt, 1.e0_rt/256.e0_rt /)
  real(rt)        , dimension(-3:4), parameter :: ff4b3 = &
       (/ -4.e0_rt/256.e0_rt, 22.e0_rt/256.e0_rt, -52.e0_rt/256.e0_rt, &
       325.e0_rt/256.e0_rt, -56.e0_rt/256.e0_rt, 28.e0_rt/256.e0_rt, -8.e0_rt/256.e0_rt, 1.e0_rt/256.e0_rt /)

end module filter_module

max_stars_repo_name: yingtchen/Castro

max_issues_repo_name: yingtchen/Castro

max_forks_repo_name: yingtchen/Castro


text:
module irc_data

implicit none

logical :: GAMESS, Molpro
integer :: m
character(len=80) :: g_irc_output_reac, g_irc_output_prod, & ! GAMESS files
                     m_irc_output                            ! Molpro files

integer                       :: n1, n2, numbr
double precision, allocatable :: geom_reac(:,:,:), geom_prod(:,:,:), en_reac(:), en_prod(:), &
                                 irc_reac(:), irc_prod(:), ref(:,:)
double precision              :: l(5)=0.0d0,h(5)=0.0d0

contains
!-------------------------------------------------------------------
subroutine read_input_file(file_inp)

implicit none
character(len=*),intent(in)              :: file_inp
namelist /programs/ GAMESS, Molpro
namelist /gamess_input_files/ g_irc_output_reac, g_irc_output_prod
namelist /molpro_input_files/ m_irc_output
!-------------------------------------------------------------------
open (unit = 99, file = file_inp)
read(unit = 99, nml = programs)

if (GAMESS .and. (.not. Molpro)) then
 m = 1
 write(*,*) "Working with GAMESS input files"
 read(unit = 99, nml = gamess_input_files)
else if (Molpro .and. (.not. GAMESS)) then
 m = 2
 write(*,*) "Working with Molpro input files"
 read(unit = 99, nml = molpro_input_files)
else if ((GAMESS .and. Molpro)) then
  write(*,*) "Error! You chose both GAMESS and Molpro to be TRUE. &
             Check your input and return"
  call exit
else if ((.not. GAMESS) .and. (.not. Molpro)) then
  write(*,*) "Error! You chose both GAMESS and Molpro to be FALSE.  &
             Check your input and return"
  call exit
end if

end subroutine read_input_file
!----------------------------------------------------------------------
end module irc_data

max_stars_repo_name: svarganov/NAST

max_issues_repo_name: svarganov/NAST

max_forks_repo_name: svarganov/NAST


text:
      SUBROUTINE slAMPQ (RA, DA, AMPRMS, RM, DM)
*+
*     - - - - - -
*      A M P Q
*     - - - - - -
*
*  Convert star RA,Dec from geocentric apparent to mean place
*
*  The mean coordinate system is the post IAU 1976 system,
*  loosely called FK5.
*
*  Use of this routine is appropriate when efficiency is important
*  and where many star positions are all to be transformed for
*  one epoch and equinox.  The star-independent parameters can be
*  obtained by calling the slMAPA routine.
*
*  Given:
*     RA       d      apparent RA (radians)
*     DA       d      apparent Dec (radians)
*
*     AMPRMS   d(21)  star-independent mean-to-apparent parameters:
*
*       (1)      time interval for proper motion (Julian years)
*       (2-4)    barycentric position of the Earth (AU)
*       (5-7)    heliocentric direction of the Earth (unit vector)
*       (8)      (grav rad Sun)*2/(Sun-Earth distance)
*       (9-11)   ABV: barycentric Earth velocity in units of c
*       (12)     sqrt(1-v**2) where v=modulus(ABV)
*       (13-21)  precession/nutation (3,3) matrix
*
*  Returned:
*     RM       d      mean RA (radians)
*     DM       d      mean Dec (radians)
*
*  References:
*     1984 Astronomical Almanac, pp B39-B41.
*     (also Lederle & Schwan, Astron. Astrophys. 134,
*      1-6, 1984)
*
*  Note:
*
*     Iterative techniques are used for the aberration and
*     light deflection corrections so that the routines
*     slAMP (or slAMPQ) and slMAP (or slMAPQ) are
*     accurate inverses;  even at the edge of the Sun's disc
*     the discrepancy is only about 1 nanoarcsecond.
*
*  Called:  slDS2C, slDIMV, slDVDV, slDVN, slDC2S,
*           slDA2P
*
*  P.T.Wallace   Starlink   7 May 2000
*
*  Copyright (C) 2000 Rutherford Appleton Laboratory
*
*  License:
*    This program is free software; you can redistribute it and/or modify
*    it under the terms of the GNU General Public License as published by
*    the Free Software Foundation; either version 2 of the License, or
*    (at your option) any later version.
*
*    This program is distributed in the hope that it will be useful,
*    but WITHOUT ANY WARRANTY; without even the implied warranty of
*    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
*    GNU General Public License for more details.
*
*    You should have received a copy of the GNU General Public License
*    along with this program (see SLA_CONDITIONS); if not, write to the
*    Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor,
*    Boston, MA  02110-1301  USA
*
*  Copyright (C) 1995 Association of Universities for Research in Astronomy Inc.
*-

      IMPLICIT NONE

      DOUBLE PRECISION RA,DA,AMPRMS(21),RM,DM

      INTEGER I,J

      DOUBLE PRECISION GR2E,AB1,EHN(3),ABV(3),P3(3),P2(3),
     :                 AB1P1,P1DV,P1DVP1,P1(3),W,PDE,PDEP1,P(3)

      DOUBLE PRECISION slDVDV,slDA2P

*  Unpack scalar and vector parameters
      GR2E = AMPRMS(8)
      AB1 = AMPRMS(12)
      DO I=1,3
         EHN(I) = AMPRMS(I+4)
         ABV(I) = AMPRMS(I+8)
      END DO

*  Apparent RA,Dec to Cartesian
      CALL slDS2C(RA,DA,P3)

*  Precession and nutation
      CALL slDIMV(AMPRMS(13),P3,P2)

*  Aberration
      AB1P1 = AB1+1D0
      DO I=1,3
         P1(I) = P2(I)
      END DO
      DO J=1,2
         P1DV = slDVDV(P1,ABV)
         P1DVP1 = 1D0+P1DV
         W = 1D0+P1DV/AB1P1
         DO I=1,3
            P1(I) = (P1DVP1*P2(I)-W*ABV(I))/AB1
         END DO
         CALL slDVN(P1,P3,W)
         DO I=1,3
            P1(I) = P3(I)
         END DO
      END DO

*  Light deflection
      DO I=1,3
         P(I) = P1(I)
      END DO
      DO J=1,5
         PDE = slDVDV(P,EHN)
         PDEP1 = 1D0+PDE
         W = PDEP1-GR2E*PDE
         DO I=1,3
            P(I) = (PDEP1*P1(I)-GR2E*EHN(I))/W
         END DO
         CALL slDVN(P,P2,W)
         DO I=1,3
            P(I) = P2(I)
         END DO
      END DO

*  Mean RA,Dec
      CALL slDC2S(P,RM,DM)
      RM = slDA2P(RM)

      END

max_stars_repo_name: ysBach/irafdocgen

max_issues_repo_name: kirxkirx/iraf

max_forks_repo_name: kirxkirx/iraf


text:
      subroutine dlk_comp_joint_surv(bev,lnu,xiv,phi,psiv,k,piv,Pi,si2,
     c typ,n,M,npX,XM,npW,W,YM,PM,btv,Pp1,dev,tv,dlk)

      integer k,n,M,npX,npW,i,mi,m1,dev(n),PM(n,M)
      double precision bev(npX),lnu,xiv(k),phi,psiv(npW),piv(k),Pi(k)
      double precision si2,XM(n,M,npX),YM(n,M),W(n,npW),btv(M)
      double precision Pp1(k,M,n),tv(n),btvnu(M),dbtvnu(M-1)
      double precision l2,l3,tmp00(k),Wp(n),tmp0(k),muv(k),etmp0(k),lk
      double precision nu,dl2bev(npX),dl3nu,dl2xiv(k),dl3xiv(k),dl3phi
      double precision dl3psiv(npW),tmp1(k),tmp(k),dlk(npX+1+k+1+npW)
      character (len=4) typ

c preliminaries
      nu = exp(lnu)
      btvnu = btv**nu
	  dbtvnu = btvnu(2:M)-btvnu(1:(M-1))
c compute functions
      dl2bev = 0
      dl3nu = 0
      dl2xiv = 0
      dl3xiv = 0
      dl3phi = 0
      dl3psiv = 0
      do i=1,n
		tmp0 = exp(xiv*phi+sum(W(i,:)*psiv))
		mi = count(btv.le.tv(i))
		do m1=1,mi
		  if(PM(i,m1)==1) then
			muv = xiv+sum(XM(i,m1,:)*bev)
			if (typ=="norm") then
              dl2bev = dl2bev+XM(i,m1,:)*sum(((YM(i,m1)-muv)/si2)*
     c                 Pp1(:,m1,i))
              dl2xiv = dl2xiv+(YM(i,m1)-muv)/si2*Pp1(:,m1,i)
			end if
			if (typ =="bino") then
			  muv = exp(muv)/(1+exp(muv))
			  tmp1 = (YM(i,m1)-muv)*Pp1(:,m1,i)
			  dl2bev = dl2bev+XM(i,m1,:)*sum(tmp1)
			  dl2xiv = dl2xiv+tmp1
			end if
		  end if
		end do
		if (mi>1) then
		  do m1=2,mi
            if (btv(m1-1)>0) then
		      tmp = log(btv(m1-1))
		    else
		      tmp = 0
		    end if
		    dl3nu = dl3nu-sum((tmp0*(log(btv(m1))*btvnu(m1)-
     c   		    tmp*btvnu(m1-1)))*Pp1(:,m1-1,i))
		    dl3xiv = dl3xiv-(tmp0*dbtvnu(m1-1))*phi*Pp1(:,m1-1,i)
		    dl3phi = dl3phi-sum((tmp0*dbtvnu(m1-1)*xiv)*Pp1(:,m1-1,i))
		    dl3psiv = dl3psiv-W(i,:)*sum((tmp0*dbtvnu(m1-1))*
     c                Pp1(:,m1-1,i))
		  end do
		end if
		if (btv(mi)>0) then
		  tmp = log(btv(mi))
		else
		  tmp = 0
		end if
		if (dev(i)==1) then
          dl3nu = dl3nu+sum(((1/nu+log(tv(i)))-tmp0*(log(tv(i))*
     c            tv(i)**nu-tmp*btvnu(mi)))*Pp1(:,mi,i))
		  dl3xiv = dl3xiv+(phi-tmp0*(tv(i)**nu-btvnu(mi))*phi)*
     c             Pp1(:,mi,i)
		  dl3phi = dl3phi+sum((xiv-tmp0*(tv(i)**nu-btvnu(mi))*
     c             xiv)*Pp1(:,mi,i))
		  dl3psiv = dl3psiv+W(i,:)*sum((1-tmp0*(tv(i)**nu-
     c              btvnu(mi)))*Pp1(:,mi,i))
		else
         dl3nu = dl3nu-sum((tmp0*(log(tv(i))*tv(i)**nu-tmp*btvnu(mi)))*
     c            Pp1(:,mi,i))
		  dl3xiv = dl3xiv-(tmp0*(tv(i)**nu-btvnu(mi))*phi)*Pp1(:,mi,i)
		  dl3phi = dl3phi-sum((tmp0*(tv(i)**nu-btvnu(mi))*xiv)*
     c             Pp1(:,mi,i))
		  dl3psiv = dl3psiv-W(i,:)*sum((tmp0*(tv(i)**nu-btvnu(mi)))*
     c              Pp1(:,mi,i))
		end if
      end do
c output
      dlk(1:npX) = dl2bev
      dlk(npX+1) = dl3nu*nu
      dlk((npX+2):(k+npX+1)) = dl2xiv+dl3xiv
      dlk(k+npX+2) = dl3phi
      dlk((k+npX+3):(k+npX+2+npW)) = dl3psiv

      end

max_stars_repo_name: afarcome/lmjm

max_issues_repo_name: afarcome/lmjm

max_forks_repo_name: afarcome/lmjm


text:
*DECK DEWSET
      SUBROUTINE DEWSET (N, ITOL, RTOL, ATOL, YCUR, EWT)
C***BEGIN PROLOGUE  DEWSET
C***SUBSIDIARY
C***PURPOSE  Set error weight vector.
C***TYPE      DOUBLE PRECISION (SEWSET-S, DEWSET-D)
C***AUTHOR  Hindmarsh, Alan C., (LLNL)
C***DESCRIPTION
C
C  This subroutine sets the error weight vector EWT according to
C      EWT(i) = RTOL(i)*ABS(YCUR(i)) + ATOL(i),  i = 1,...,N,
C  with the subscript on RTOL and/or ATOL possibly replaced by 1 above,
C  depending on the value of ITOL.
C
C***SEE ALSO  DLSODE
C***ROUTINES CALLED  (NONE)
C***REVISION HISTORY  (YYMMDD)
C   791129  DATE WRITTEN
C   890501  Modified prologue to SLATEC/LDOC format.  (FNF)
C   890503  Minor cosmetic changes.  (FNF)
C   930809  Renamed to allow single/double precision versions. (ACH)
C***END PROLOGUE  DEWSET
C**End
      INTEGER N, ITOL
      INTEGER I
      DOUBLE PRECISION RTOL, ATOL, YCUR, EWT
      DIMENSION RTOL(*), ATOL(*), YCUR(N), EWT(N)
C
C***FIRST EXECUTABLE STATEMENT  DEWSET
      GO TO (10, 20, 30, 40), ITOL
 10   CONTINUE
      DO 15 I = 1,N
 15     EWT(I) = RTOL(1)*ABS(YCUR(I)) + ATOL(1)
      RETURN
 20   CONTINUE
      DO 25 I = 1,N
 25     EWT(I) = RTOL(1)*ABS(YCUR(I)) + ATOL(I)
      RETURN
 30   CONTINUE
      DO 35 I = 1,N
 35     EWT(I) = RTOL(I)*ABS(YCUR(I)) + ATOL(1)
      RETURN
 40   CONTINUE
      DO 45 I = 1,N
 45     EWT(I) = RTOL(I)*ABS(YCUR(I)) + ATOL(I)
      RETURN
C----------------------- END OF SUBROUTINE DEWSET ----------------------
      END


text:
      PROGRAM EX_TSORT
      IMPLICIT NONE
      INTEGER NL,ND,NC,NO,IDEP,IORD,IPOS,ICODE,I,J,IL,IR
      PARAMETER(NL=15,ND=44,NC=69)
      CHARACTER*(20) LABEL
      DIMENSION IDEP(ND,2),LABEL(NL),IORD(NL),IPOS(NL),ICODE(NC)
      DATA LABEL/'DES_SYSTEM_LIB','DW01','DW02','DW03','DW04','DW05',
     1 'DW06','DW07','DWARE','GTECH','RAMLIB','STD_CELL_LIB','SYNOPSYS',
     2 'STD','IEEE'/
      DATA ICODE/1,14,13,12,1,3,2,11,15,0,2,15,2,9,10,0,3,15,3,9,0,4,14,
     213,9,4,3,2,15,10,0,5,5,15,2,9,10,0,6,6,15,9,0,7,7,15,9,0,8,15,9,0,
     39,15,9,0,10,15,10,0,11,14,15,0,12,15,12,0,0/

C DECODE DEPENDENCIES AND BUILD IDEP ARRAY
      I=0
      J=0
   10 I=I+1
      IL=ICODE(I)
      IF(IL.EQ.0) GO TO 30
   20 I=I+1
      IR=ICODE(I)
      IF(IR.EQ.0) GO TO 10
      J=J+1
      IDEP(J,1)=IL
      IDEP(J,2)=IR
      GO TO 20
   30 CONTINUE

C SORT LIBRARIES ACCORDING TO DEPENDENCIES (TOPOLOGICAL SORT)
      CALL TSORT(NL,ND,IDEP,IORD,IPOS,NO)

      PRINT*,'COMPILE ORDER'
      DO 40 I=1,NO
   40 PRINT*,LABEL(IORD(I))
      PRINT*,'UNORDERED LIBRARIES'
      DO 50 I=NO+1,NL
   50 PRINT*,LABEL(IORD(I))
      END

max_issues_repo_name: seanwallawalla-forks/RosettaCodeData

max_forks_repo_name: seanwallawalla-forks/RosettaCodeData


text:
PROGRAM create_dependent_set_tape5

  USE type_kinds
  USE error_handler
  USE lblrtm_utility
  USE dependent_profile_set
  USE meteorology

  IMPLICIT NONE

  CHARACTER( LEN = 24 ) :: database_file = 'dependent_profile_set.nc'

  INTEGER( Long ) :: i
  INTEGER( Long ) :: error_status
  INTEGER( Long ) :: n_profiles
  INTEGER( Long ) :: n_levels
  INTEGER( Long ) :: n_absorbers
  INTEGER( Long ) :: profile
  INTEGER( Long ) :: absorber

  REAL( Single )  :: surface_altitude

! -- Profile data arrays
  REAL( Double ),       ALLOCATABLE, DIMENSION( : )       :: pressure
  REAL( Double ),       ALLOCATABLE, DIMENSION( :, : )    :: temperature
  REAL( Double ),       ALLOCATABLE, DIMENSION( :, :, : ) :: absorber_amount

  REAL( Single ),       ALLOCATABLE, DIMENSION( : )       :: profile_t
  REAL( Single ),       ALLOCATABLE, DIMENSION( : )       :: profile_mr
  REAL( Single ),       ALLOCATABLE, DIMENSION( : )       :: altitude
  REAL( Single ),       ALLOCATABLE, DIMENSION( :, : )    :: profile_a

  INTEGER( Long ),      ALLOCATABLE, DIMENSION( : )       :: absorber_id
  INTEGER( Long ),      ALLOCATABLE, DIMENSION( : )       :: profile_number

  CHARACTER( LEN = 1 ), ALLOCATABLE, DIMENSION( : )       :: absorber_units

! -- Spectral band parameters for LBLRTM runs
  REAL( Single ), PARAMETER :: v_begin   = 1000.0
  REAL( Single ), PARAMETER :: v_end     = 2000.0

  INTEGER( Long ) :: profile_begin
  INTEGER( Long ) :: profile_end
  INTEGER( Long ) :: tmp

  REAL( Single ) :: v1
  REAL( Single ) :: v2
  REAL( Single ) :: za

  CHARACTER( LEN = 78 ) :: tape5_filename
  CHARACTER( LEN = 78 ) :: tape5_header

! -- RCS Id

  CHARACTER( LEN = 128 ), PARAMETER :: rcs_Id = &

!-------------------------------------------------------------------------------
!                         -- Create the TAPE5 files! --
!-------------------------------------------------------------------------------

! --------------------
! Inquire the database
! --------------------

  error_status = inquire_profile_set( database_file, &
                                      n_profiles  = n_profiles, &
                                      n_absorbers = n_absorbers, &
                                      n_levels    = n_levels )

  IF ( error_status /= 1 ) THEN
    WRITE( *, '( /5x, "Error inquiring file." )' )
    STOP
  END IF

! -- Output results
  WRITE( *, '( /5x, "n_profiles           : ", i6, &
              &/5x, "n_absorbers          : ", i6, &
              &/5x, "n_levels             : ", i6 )' ) &
            n_profiles, &
            n_absorbers, &
            n_levels

! ---------------
! Allocate arrays
! ---------------

  ALLOCATE( pressure( n_levels ), &
            temperature( n_levels, n_profiles ), &
            profile_t( n_levels ), &
            absorber_amount( n_levels, n_absorbers, n_profiles), &
            profile_a( n_levels, n_absorbers ), &
            profile_mr( n_levels ), &
            altitude( n_levels ), &
            absorber_id( n_absorbers ), &
            absorber_units( n_absorbers ), &
            profile_number( n_profiles ) )

! -------------
! Read the data
! -------------

  error_status = read_profile_set( database_file, &
                                   pressure        = pressure, &
                                   temperature     = temperature, &
                                   absorber_amount = absorber_amount, &
                                   absorber_id     = absorber_id, &
                                   profile_number  = profile_number )

  IF ( error_status /= 1 ) THEN
    WRITE( *, '( /5x, "Error reading file." )' )
    STOP
  END IF

! --------------------------------
! Assign the absorber units (ppmv)
! --------------------------------

  DO i = 1, n_absorbers
    absorber_units( i ) = 'A'
  END DO

! ----------------------------------------------------
! Determine the number of profiles to begin and end at
! ----------------------------------------------------

  WRITE( *, '( /5x, "Enter profile # limits [0,",i2,"] : " )', &
            ADVANCE = 'NO' ) n_profiles - 1
  READ( *, * ) profile_begin, profile_end

  profile_begin = MAX( 0,              profile_begin )
  profile_end   = MIN( n_profiles - 1, profile_end   )

  IF ( profile_begin > profile_end ) THEN
    tmp = profile_begin
    profile_begin = profile_end
    profile_end   = tmp
  END IF

! --------------------------------
! Assign zenith angle - in DEGREES
! --------------------------------

  za = 0.0

! ------------------
! Loop over profiles
! ------------------

  profile_loop: DO profile = profile_begin, profile_end

!   -- Separate out the current profiles
    profile_t = REAL( temperature( :, profile + 1 ), Single )
    profile_a = REAL( absorber_amount( :, :, profile + 1 ), Single )

!   -- Convert h2o ppmv to g/kg
    profile_mr = 1.0e-03 * profile_a( :, 1 ) * 18.0153 / 28.9658

!   -- Calculate the geopotential altitudes
    surface_altitude = 0.0
    altitude = geopotential_height( REAL( pressure, Single ), &
                                    profile_t,  &
                                    profile_mr,  &
                                    surface_altitude )

!   -- Create a TAPE5 file directory and name
    WRITE( tape5_filename, '( "tape5.profile", i2.2, ".rdk" )' ) &
                           profile

!   -- Create a header string
    WRITE( tape5_header, '( "X" )' )

!   -- Create the TAPE5 file
    error_status = lblrtm_create_tape5( REAL( pressure, Single ), &
                                        profile_t, &
                                        profile_a, &
                                        absorber_units, &
                                        absorber_id, &
                                        altitude, &
                                        surface_altitude, &
                                        v_begin, v_end, &
                                        zenith_angle = za, &
                                        calculation_levels = altitude, &
                                        continuum_option = 0, &
                                        merge_option = 3, &
                                        filename = './tape5_files/' // &
                                                   TRIM( tape5_filename ), &
                                        header   = TRIM( tape5_header ) )

    IF ( error_status /= success_state ) THEN
      PRINT *,'Error creating ', tape5_filename
      EXIT profile_loop
    END IF

  END DO profile_loop

END PROGRAM create_dependent_set_tape5

max_stars_repo_name: hsbadr/crtm

max_issues_repo_name: hsbadr/crtm

max_forks_repo_name: hsbadr/crtm


text:
use IOClass
use ArrayClass
implicit none

type(IO_) :: f
type(Math_) :: math
real(real64),allocatable :: t(:),ut(:),vt(:),at(:)
real(real64) :: w

! Frequency
w = 1.0d0

t = linspace([0.0d0,1.0d0], 1024)

vt = 2.0d0*t

ut = I_dx(x=t,fx=vt,f0=0.0d0)
at = d_dx(x=t,fx=vt)

call f%open("u.txt","w")
call f%write(t,ut)
call f%close()
call f%plot(option="with lines")

call f%open("a.txt","w")
call f%write(t,at)
call f%close()
call f%plot(option="with lines")

call f%open("v.txt","w")
call f%write(t,vt)
call f%close()
call f%plot(option="with lines")

end

max_stars_repo_name: kazulagi/plantFEM


SUBROUTINE formlump(diag,emm,g)
!
! This subroutine forms the lumped global mass matrix as a vector diag.
!
 IMPLICIT NONE
 INTEGER,PARAMETER::iwp=SELECTED_REAL_KIND(15)
 REAL(iwp),INTENT(IN)::emm(:,:)
 REAL(iwp),INTENT(OUT)::diag(0:)
 INTEGER,INTENT(IN)::g(:)
 INTEGER::i,ndof
 ndof=UBOUND(emm,1)
 DO i=1,ndof
   diag(g(i))=diag(g(i))+emm(i,i)
 END DO
RETURN
END SUBROUTINE formlump


text:
      PROGRAM RUNLEN
C Run Length Encoding, February 9, 2021
C
C Copyright 2021 Dennis Decker Jensen
C
C https://programmingpraxis.com/2021/02/09/run-length-encoding-2/2/
C
C Tectonics: gfortran -o runlen runlen.f
C Execute: ./runlen
C
      CHARACTER*(*) STR
      PARAMETER (STR='aaaabbbcca')
      EXTERNAL RUNLEN1, RUNLEN2
      CALL RUNLEN1(STR)
      CALL RUNLEN2(STR)
      END

      SUBROUTINE RUNLEN1(STR)
C Uses a small state machine and a direct counter
      CHARACTER*(*) STR
      CHARACTER C

      IF (LEN(STR).EQ.0) THEN
         WRITE (*, *) '[]'
      ELSE
         WRITE (*, *) '['
         N=1
         C=STR(1:1)
         DO 10 I=2,LEN(STR)
            IF (C.EQ.STR(I:I)) THEN
               N=N+1
            ELSE
               WRITE (*,5) C, N
    5          FORMAT (' ''',A1, ''': ',I0)
               N=1
               C=STR(I:I)
            ENDIF
   10    CONTINUE
         WRITE (*,5) C, N
         WRITE (*,*) ']'
      END IF
      END

      SUBROUTINE RUNLEN2(STR)
C Uses straight structured programming and an indirect counter
      CHARACTER*(*) STR

      J=1
      WRITE (*,*) '['
   10 I=J
      J=I+1
      IF (I.GT.LEN(STR)) GO TO 40
   20 IF (J.GT.LEN(STR)) GO TO 30
      IF (STR(I:I).EQ.STR(J:J)) THEN
         J=J+1
         GO TO 20
      ELSE
         WRITE (*,100) STR(I:I), J-I
         GO TO 10
      END IF
   30 WRITE (*,100) STR(I:I), J-I
   40 WRITE (*,*) ']'
  100 FORMAT (' ''',A1, ''': ',I0)
      END

max_stars_repo_name: dennisdjensen/Sketchbook

max_issues_repo_name: dennisdjensen/sketchbook

max_forks_repo_name: dennisdjensen/sketchbook


text:
+incdir+./verif/
+incdir+./verif/agents/apb_mstr_agent/
+incdir+./verif/env/apb_mstr_env/
+incdir+./verif/env/dma_perf/
+incdir+./verif/sequence_lib/apb_mstr_sequence_lib/
+incdir+./verif/tests/
verif/agents/apb_mstr_agent/apb_interface.sv
verif/env/dma_perf/dma_perf_interface.sv
verif/agents/apb_mstr_agent/apb_agent_pkg.sv
verif/env/dma_perf/dma_perf_pkg.sv
verif/env/apb_mstr_env/apb_env_pkg.sv
verif/sequence_lib/apb_mstr_sequence_lib/apb_seq_lib_pkg.sv
verif/tests/apb_test_pkg.sv


text:
The Samurai Jogger Lady can be seen practicing her wicked beautiful moves early in the morning at John Barovetto Park or jogging along the greenbelt between Arroyo Avenue and 2nd Street.

20100731 10:53:50 nbsp Not a samurai.  Bet you its http://en.wikipedia.org/wiki/Wushu Wushu.  UCD has a competitive team.  Users/EdWins

20100801 03:17:56 nbsp Im kind of surprised she didnt get into trouble.  Id figure someone would have pegged her for Brandishing a weapon or something stupid.   Users/MasonMurray
  Out here in Mace Ranch, the outback of Davis, we live a more rugged and independent lifestyle.  Besides, I reckon shell be at the head of the pack when we have to drive the zombies back downtown. users/jimstewart

20110811 12:15:17 nbsp Is this also the lady who walks her cat past central park in the morning? Users/MikeyCrews

20110811 14:54:33 nbsp Hasnt she already signed on for Kill Bill Part 3? Users/OldDavis73

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


MODULE module_sf_gfdl

CONTAINS

   SUBROUTINE SF_GFDL(U3D,V3D,T3D,QV3D,P3D, &
                     CP,ROVCP,R,XLV,PSFC,CHS,CHS2,CQS2, CPM, &
                     DT, SMOIS,num_soil_layers,ISLTYP,ZNT, &

                     UST,PSIM,PSIH, &
                     XLAND,HFX,QFX,TAUX,TAUY,LH,GSW,GLW,TSK,FLHC,FLQC, &
                     QGH,QSFC,U10,V10, &
                     GZ1OZ0,WSPD,BR,ISFFLX, &
                     EP1,EP2,KARMAN,NTSFLG,SFENTH, &
                     ids,ide, jds,jde, kds,kde, &
                     ims,ime, jms,jme, kms,kme, &
                     its,ite, jts,jte, kts,kte )

      USE MODULE_GFS_MACHINE, ONLY : kind_phys
      USE MODULE_GFS_FUNCPHYS , ONLY : gfuncphys,fpvs
      USE MODULE_GFS_PHYSCONS, grav => con_g

      IMPLICIT NONE
      character*255 :: message
      INTEGER, INTENT(IN) :: ids,ide, jds,jde, kds,kde, &
                                        ims,ime, jms,jme, kms,kme, &
                                        its,ite, jts,jte, kts,kte, &
                                        ISFFLX,NUM_SOIL_LAYERS,NTSFLG
      REAL, INTENT(IN) :: &
                                        CP, &
                                        EP1, &
                                        EP2, &
                                        KARMAN, &
                                        R, &
                                        ROVCP, &
                                        DT, &
                                        SFENTH, &
                                        XLV
      REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(IN) :: &
                                        P3D, &
                                        QV3D, &
                                        T3D, &
                                        U3D, &
                                        V3D
      INTEGER, DIMENSION( ims:ime , jms:jme ), INTENT(IN ):: ISLTYP
      REAL, DIMENSION( ims:ime , 1:num_soil_layers, jms:jme ), INTENT(INOUT):: SMOIS
      REAL, DIMENSION(ims:ime, jms:jme), INTENT(IN) :: &
                                        PSFC, &
                                        GLW, &
                                        GSW, &
                                        XLAND
      REAL, DIMENSION(ims:ime, jms:jme), INTENT(INOUT) :: &
                                        TSK, &
                                        BR, &
                                        CHS, &
                                        CHS2, &
                                        CPM, &
                                        CQS2, &
                                        FLHC, &
                                        FLQC, &
                                        GZ1OZ0, &
                                        HFX, &
                                        LH, &
                                        PSIM, &
                                        PSIH, &
                                        QFX, &
                                        QGH, &
                                        QSFC, &
                                        UST, &
                                        ZNT, &
                                        WSPD, &
                                        TAUX, &
                                        TAUY
      REAL, DIMENSION(ims:ime, jms:jme), INTENT(OUT) :: &
                                        U10, &
                                        V10
      REAL :: ESAT, &
                                        cpcgs, &
                                        smc, &
                                        smcdry, &
                                        smcmax
      REAL (kind=kind_phys) :: &
                                        RHOX
      REAL, DIMENSION(1:30) :: MAXSMC, &
                                        DRYSMC
      REAL (kind=kind_phys), DIMENSION(its:ite) :: &
                                        CH, &
                                        CM, &
                                        DDVEL, &
                                        DRAIN, &
                                        EP, &
                                        EVAP, &
                                        FH, &
                                        FH2, &
                                        FM, &
                                        HFLX, &
                                        PH, &
                                        PM, &
                                        PRSL1, &
                                        PRSLKI, &
                                        PS, &
                                        Q1, &
                                        Q2M, &
                                        QSS, &
                                        QSURF, &
                                        RB, &
                                        RCL, &
                                        RHO1, &
                                        SLIMSK, &
                                        STRESS, &
                                        T1, &
                                        T2M, &
                                        THGB, &
                                        THX, &
                                        TSKIN, &
                                        SHELEG, &
                                        U1, &
                                        U10M, &
                                        USTAR, &
                                        V1, &
                                        V10M, &
                                        WIND, &
                                        Z0RL, &
                                        Z1
      REAL, DIMENSION(kms:kme, ims:ime) :: &
                                        rpc, &
                                        tpc, &
                                        upc, &
                                        vpc
     REAL, DIMENSION(ims:ime) :: &
                                        pspc, &
                                        pkmax, &
                                        tstrc, &
                                        zoc, &
                                        mzoc, &
                                        wetc, &
                                        slwdc, &
                                        rib, &
                                        zkmax, &
                                        tkmax, &
                                        fxmx, &
                                        fxmy, &
                                        cdm, &
                                        fxh, &
                                        fxe, &
                                        xxfh, &
                                        xxfh2, &
                                        wind10, &
                                        tjloc
      INTEGER :: &
                                        I, &
                                        II, &
                                        IGPVS, &
                                        IM, &
                                        J, &
                                        K, &
                                        KM
        DATA MAXSMC/0.339, 0.421, 0.434, 0.476, 0.476, 0.439, &
                    0.404, 0.464, 0.465, 0.406, 0.468, 0.468, &
                    0.439, 1.000, 0.200, 0.421, 0.000, 0.000, &
                    0.000, 0.000, 0.000, 0.000, 0.000, 0.000, &
                    0.000, 0.000, 0.000, 0.000, 0.000, 0.000/
        DATA DRYSMC/0.010, 0.028, 0.047, 0.084, 0.084, 0.066, &
                    0.067, 0.120, 0.103, 0.100, 0.126, 0.138, &
                    0.066, 0.000, 0.006, 0.028, 0.000, 0.000, &
                    0.000, 0.000, 0.000, 0.000, 0.000, 0.000, &
                    0.000, 0.000, 0.000, 0.000, 0.000, 0.000/
      DATA IGPVS/0/
      save igpvs
   if(igpvs.eq.0) then
   endif
   igpvs=1
   IM=ITE-ITS+1
   KM=KTE-KTS+1
   WRITE(message,*)'WITHIN THE GFDL SCHEME, NTSFLG=1 FOR GFDL SLAB  2010 UPGRADS',NTSFLG
   call wrf_debug(1,message)
   DO J=jts,jte
      DO i=its,ite
        DDVEL(I)=0.
        RCL(i)=1.
        PRSL1(i)=P3D(i,kts,j)*.001
         wetc(i)=1.0
         if(xland(i,j).lt.1.99) then
         smc=smois(i,1,j)
         smcdry=drysmc(isltyp(i,j))
         smcmax=maxsmc(isltyp(i,j))
         wetc(i)=(smc-smcdry)/(smcmax-smcdry)
         wetc(i)=amin1(1.,amax1(wetc(i),0.))
         endif
        pspc(i)=PSFC(i,j)*10.
        pkmax(i)=P3D(i,kts,j)*10.
        PS(i)=PSFC(i,j)*.001
        Q1(I) = QV3D(i,kts,j)
        rpc(kts,i)=QV3D(i,kts,j)
        QSURF(I)=0.
        SHELEG(I)=0.
        SLIMSK(i)=ABS(XLAND(i,j)-2.)
        TSKIN(i)=TSK(i,j)
        tstrc(i)=TSK(i,j)
        T1(I) = T3D(i,kts,j)
        tpc(kts,i)=T3D(i,kts,j)
        U1(I) = U3D(i,kts,j)
        upc(kts,i)=U3D(i,kts,j) * 100.
        USTAR(I) = UST(i,j)
        V1(I) = V3D(i,kts,j)
        vpc(kts,i)=v3D(i,kts,j) * 100.
        Z0RL(I) = ZNT(i,j)*100.
        zoc(i)=ZNT(i,j)*100.
         if(XLAND(i,j).gt.1.99) zoc(i)=- zoc(i)
        slwdc(i)=gsw(i,j)+glw(i,j)
        slwdc(i)=0.239*60.*slwdc(i)*1.e-4
         tjloc(i)=float(j)
      ENDDO
      DO i=its,ite
         PRSLKI(i)=(PS(I)/PRSL1(I))**ROVCP
         THGB(I)=TSKIN(i)*(100./PS(I))**ROVCP
         THX(I)=T1(i)*(100./PRSL1(I))**ROVCP
         RHO1(I)=PRSL1(I)*1000./(R*T1(I)*(1.+EP1*Q1(I)))
         Q1(I)=Q1(I)/(1.+Q1(I))
      ENDDO
     CALL MFLUX2( fxh,fxe,fxmx,fxmy,cdm,rib,xxfh,zoc,mzoc,tstrc, &
                   pspc,pkmax,wetc,slwdc,tjloc, &
                   upc,vpc,tpc,rpc,dt,J,wind10,xxfh2,ntsflg,SFENTH, &
                   ids,ide, jds,jde, kds,kde, &
                   ims,ime, jms,jme, kms,kme, &
                   its,ite, jts,jte, kts,kte )
1010 format(2I4,9F11.6)
      DO i=its,ite
        IF(NTSFLG==1) then
        tsk(i,j) = tstrc(i)
        if(j.eq.jde) then
         tsk(i,j) = tsk(i,j-1)
         endif
        if(j.eq.jds) then
         tsk(i,j) = tsk(i,j+1)
         endif
        if(i.eq.ide) tsk(i,j) = tsk(i-1,j)
        if(i.eq.ids) tsk(i,j) = tsk(i+1,j)
        endif
        znt(i,j)= 0.01*abs(zoc(i))
        wspd(i,j) = SQRT(upc(kts,i)*upc(kts,i) + vpc(kts,i)*vpc(kts,i))
        wspd(i,j) = amax1(wspd(i,j) ,100.)/100.
        u10m(i) = u1(i)*(wind10(i)/wspd(i,j))/100.
        v10m(i) = v1(i)*(wind10(i)/wspd(i,j))/100.
        zkmax(i) = -R*tpc(kts,i)*alog(pkmax(i)/pspc(i))/grav
        gz1oz0(i,j)=alog(zkmax(i)/znt(i,j))
        ustar (i)= 0.01*sqrt(cdm(i)* &
                   (upc(kts,i)*upc(kts,i) + vpc(kts,i)*vpc(kts,i)))
        qfx (i,j)=-10.*fxe(i)
        hfx (i,j)= -10.*CP*fxh(i)
        taux (i,j)= fxmx(i)/10.
        tauy (i,j)= fxmy(i)/10.
        fm(i) = karman/sqrt(cdm(i))
        fh(i) = karman*xxfh(i)
        PSIM(i,j)=GZ1OZ0(i,j)-FM(i)
        PSIH(i,j)=GZ1OZ0(i,j)-FH(i)
        fh2(i) = karman*xxfh2(i)
        ch(i) = karman*karman/(fm(i) * fh(i))
        cm(i) = cdm(i)
        U10(i,j)=U10M(i)
        V10(i,j)=V10M(i)
        BR(i,j)=rib(i)
        CHS(I,J)=CH(I)*wspd (i,j)
        CHS2(I,J)=USTAR(I)*KARMAN/FH2(I)
        CPM(I,J)=CP*(1.+0.8*QV3D(i,kts,j))
        esat = fpvs(t1(i))
        QGH(I,J)=ep2*esat/(1000.*ps(i)-esat)
        esat = fpvs(tskin(i))
        qss(i) = ep2*esat/(1000.*ps(i)-esat)
        QSFC(I,J)=qss(i)
        UST(i,j)=ustar(i)
      ENDDO
      DO i=its,ite
        FLHC(i,j)=CPM(I,J)*RHO1(I)*CHS(I,J)
        FLQC(i,j)=RHO1(I)*CHS(I,J)
        CQS2(i,j)=CHS2(I,J)
      ENDDO
      IF (ISFFLX.EQ.0) THEN
        DO i=its,ite
          HFX(i,j)=0.
          LH(i,j)=0.
          QFX(i,j)=0.
        ENDDO
      ELSE
        DO i=its,ite
          IF(XLAND(I,J)-1.5.GT.0.)THEN
            hfx (i,j)= -10.*CP*fxh(i)
          ELSEIF(XLAND(I,J)-1.5.LT.0.)THEN
            hfx (i,j)= -10.*CP*fxh(i)
            HFX(I,J)=AMAX1(HFX(I,J),-250.)
          ENDIF
          qfx(i,j)=-10.*fxe(i)
          QFX(I,J)=AMAX1(QFX(I,J),0.)
          LH(I,J)=XLV*QFX(I,J)
        ENDDO
      ENDIF
    if(j.eq.jds.or.j.eq.jde) then
    write(message,*) "TSFC in gfdl sf mod,dt, its,ite,jts,jts", dt,its,ite,jts,jte,ids,ide,jds,jde
    call wrf_debug(1,message)
    write(message,*) "TSFC", (TSK(i,j),i=its,ite)
    call wrf_debug(1,message)
    endif
   ENDDO
   END SUBROUTINE SF_GFDL
       SUBROUTINE MFLUX2( fxh,fxe,fxmx,fxmy,cdm,rib,xxfh,zoc,mzoc,tstrc, &
                          pspc,pkmax,wetc,slwdc,tjloc, &
                          upc,vpc,tpc,rpc,dt,jfix,wind10,xxfh2,ntsflg,sfenth, &
                          ids,ide, jds,jde, kds,kde, &
                          ims,ime, jms,jme, kms,kme, &
                          its,ite, jts,jte, kts,kte )
      IMPLICIT NONE
      integer,intent(in) :: ids,ide, jds,jde, kds,kde
      integer,intent(in) :: ims,ime, jms,jme, kms,kme
      integer,intent(in) :: its,ite, jts,jte, kts,kte
      integer,intent(in) :: jfix,ntsflg
      real, intent (out), dimension (ims :ime ) :: fxh
      real, intent (out), dimension (ims :ime ) :: fxe
      real, intent (out), dimension (ims :ime ) :: fxmx
      real, intent (out), dimension (ims :ime ) :: fxmy
      real, intent (out), dimension (ims :ime ) :: cdm
      real, intent (out), dimension (ims :ime ) :: rib
      real, intent (out), dimension (ims :ime ) :: xxfh
      real, intent (out), dimension (ims :ime ) :: xxfh2
      real, intent (out), dimension (ims :ime ) :: wind10
      real, intent ( inout), dimension (ims :ime ) :: zoc,mzoc
      real, intent ( inout), dimension (ims :ime ) :: tstrc
      real, intent ( in) :: dt
      real, intent ( in) :: sfenth
      real, intent ( in), dimension (ims :ime ) :: pspc
      real, intent ( in), dimension (ims :ime ) :: pkmax
      real, intent ( in), dimension (ims :ime ) :: wetc
      real, intent ( in), dimension (ims :ime ) :: slwdc
      real, intent ( in), dimension (ims :ime ) :: tjloc
      real, intent ( in), dimension (kms:kme, ims :ime ) :: upc
      real, intent ( in), dimension (kms:kme, ims :ime ) :: vpc
      real, intent ( in), dimension (kms:kme, ims :ime ) :: tpc
      real, intent ( in), dimension (kms:kme, ims :ime ) :: rpc
      integer, parameter :: icntx = 30
      integer, dimension(1 :ime) :: ifz
      integer, dimension(1 :ime) :: indx
      integer, dimension(1 :ime) :: istb
      integer, dimension(1 :ime) :: it
      integer, dimension(1 :ime) :: iutb
      real, dimension(1 :ime) :: aap
      real, dimension(1 :ime) :: bq1
      real, dimension(1 :ime) :: bq1p
      real, dimension(1 :ime) :: delsrad
      real, dimension(1 :ime) :: ecof
      real, dimension(1 :ime) :: ecofp
      real, dimension(1 :ime) :: estso
      real, dimension(1 :ime) :: estsop
      real, dimension(1 :ime) :: fmz1
      real, dimension(1 :ime) :: fmz10
      real, dimension(1 :ime) :: fmz2
      real, dimension(1 :ime) :: fmzo1
      real, dimension(1 :ime) :: foft
      real, dimension(1 :ime) :: foftm
      real, dimension(1 :ime) :: frac
      real, dimension(1 :ime) :: land
      real, dimension(1 :ime) :: pssp
      real, dimension(1 :ime) :: qf
      real, dimension(1 :ime) :: rdiff
      real, dimension(1 :ime) :: rho
      real, dimension(1 :ime) :: rkmaxp
      real, dimension(1 :ime) :: rstso
      real, dimension(1 :ime) :: rstsop
      real, dimension(1 :ime) :: sf10
      real, dimension(1 :ime) :: sf2
      real, dimension(1 :ime) :: sfm
      real, dimension(1 :ime) :: sfzo
      real, dimension(1 :ime) :: sgzm
      real, dimension(1 :ime) :: slwa
      real, dimension(1 :ime) :: szeta
      real, dimension(1 :ime) :: szetam
      real, dimension(1 :ime) :: t1
      real, dimension(1 :ime) :: t2
      real, dimension(1 :ime) :: tab1
      real, dimension(1 :ime) :: tab2
      real, dimension(1 :ime) :: tempa1
      real, dimension(1 :ime) :: tempa2
      real, dimension(1 :ime) :: theta
      real, dimension(1 :ime) :: thetap
      real, dimension(1 :ime) :: tsg
      real, dimension(1 :ime) :: tsm
      real, dimension(1 :ime) :: tsp
      real, dimension(1 :ime) :: tss
      real, dimension(1 :ime) :: ucom
      real, dimension(1 :ime) :: uf10
      real, dimension(1 :ime) :: uf2
      real, dimension(1 :ime) :: ufh
      real, dimension(1 :ime) :: ufm
      real, dimension(1 :ime) :: ufzo
      real, dimension(1 :ime) :: ugzm
      real, dimension(1 :ime) :: uzeta
      real, dimension(1 :ime) :: uzetam
      real, dimension(1 :ime) :: vcom
      real, dimension(1 :ime) :: vrtkx
      real, dimension(1 :ime) :: vrts
      real, dimension(1 :ime) :: wind
      real, dimension(1 :ime) :: windp
      real, dimension(1 :ime) :: xxfm
      real, dimension(1 :ime) :: xxsh
      real, dimension(1 :ime) :: z10
      real, dimension(1 :ime) :: z2
      real, dimension(1 :ime) :: zeta
      real, dimension(1 :ime) :: zkmax
      real, dimension(1 :ime) :: pss
      real, dimension(1 :ime) :: tstar
      real, dimension(1 :ime) :: ukmax
      real, dimension(1 :ime) :: vkmax
      real, dimension(1 :ime) :: tkmax
      real, dimension(1 :ime) :: rkmax
      real, dimension(1 :ime) :: zot
      real, dimension(1 :ime) :: fhzo1
      real, dimension(1 :ime) :: sfh
      real :: ux13, yo, y,xo,x,ux21,ugzzo,ux11,ux12,uzetao,xnum,alll
      real :: ux1,ugz,x10,uzo,uq,ux2,ux3,xtan,xden,y10,uzet1o,ugz10
      real :: szet2, zal2,ugz2
      real :: rovcp,boycon,cmo2,psps1,zog,enrca,rca,cmo1,amask,en,ca,a,c
      real :: sgz,zal10,szet10,fmz,szo,sq,fmzo,rzeta1,zal1g,szetao,rzeta2,zal2g
      real :: hcap,xks,pith,teps,diffot,delten,alevp,psps2,alfus,nstep
      real :: shfx,sigt4,reflect
      real :: cor1,cor2,szetho,zal2gh,cons_p000001,cons_7,vis,ustar,restar,rat
      real :: wndm,ckg
      real :: yz,y1,y2,y3,y4,windmks,znott,znotm
      integer:: i,j,ii,iq,nnest,icnt,ngd,ip
      real, dimension (223) :: tab
      real, dimension (223) :: table
      real, dimension (101) :: tab11
      real, dimension (41) :: table4
      real, dimension (42) :: tab3
      real, dimension (54) :: table2
      real, dimension (54) :: table3
      real, dimension (74) :: table1
      real, dimension (80) :: tab22
      equivalence (tab(1),tab11(1))
      equivalence (tab(102),tab22(1))
      equivalence (tab(182),tab3(1))
      equivalence (table(1),table1(1))
      equivalence (table(75),table2(1))
      equivalence (table(129),table3(1))
      equivalence (table(183),table4(1))
      data amask/ -98.0/
      data tab11/21*0.01403,0.01719,0.02101,0.02561,0.03117,0.03784, &
     &.04584,.05542,.06685,.08049,.09672,.1160,.1388,.1658,.1977,.2353, &
     &.2796,.3316,.3925,.4638,.5472,.6444,.7577,.8894,1.042,1.220,1.425, &
     &1.662,1.936,2.252,2.615,3.032,3.511,4.060,4.688,5.406,6.225,7.159, &
     &8.223,9.432,10.80,12.36,14.13,16.12,18.38,20.92,23.80,27.03,30.67, &
     &34.76,39.35,44.49,50.26,56.71,63.93,71.98,80.97,90.98,102.1,114.5, &
     &128.3,143.6,160.6,179.4,200.2,223.3,248.8,276.9,307.9,342.1,379.8, &
     &421.3,466.9,517.0,572.0,632.3,698.5,770.9,850.2,937.0,1032./
      data tab22/1146.6,1272.0,1408.1,1556.7,1716.9,1890.3,2077.6,2279.6 &
     &,2496.7,2729.8,2980.0,3247.8,3534.1,3839.8,4164.8,4510.5,4876.9, &
     &5265.1,5675.2,6107.8,6566.2,7054.7,7575.3,8129.4,8719.2,9346.5, &
     &10013.,10722.,11474.,12272.,13119.,14017.,14969.,15977.,17044., &
     &18173.,19367.,20630.,21964.,23373.,24861.,26430.,28086.,29831., &
     &31671.,33608.,35649.,37796.,40055.,42430.,44927.,47551.,50307., &
     &53200.,56236.,59422.,62762.,66264.,69934.,73777.,77802.,82015., &
     &86423.,91034.,95855.,100890.,106160.,111660.,117400.,123400., &
     &129650.,136170.,142980.,150070.,157460.,165160.,173180.,181530., &
     &190220.,199260./
      data tab3/208670.,218450.,228610.,239180.,250160.,261560.,273400., &
     &285700.,298450.,311690.,325420.,339650.,354410.,369710.,385560., &
     &401980.,418980.,436590.,454810.,473670.,493170.,513350.,534220., &
     &555800.,578090.,601130.,624940.,649530.,674920.,701130.,728190., &
     &756110.,784920.,814630.,845280.,876880.,909450.,943020.,977610., &
     &1013250.,1049940.,1087740./
      data table1/20*0.0,.3160e-02,.3820e-02,.4600e-02,.5560e-02,.6670e-02, &
     & .8000e-02,.9580e-02,.1143e-01,.1364e-01,.1623e-01,.1928e-01, &
     &.2280e-01,.2700e-01,.3190e-01,.3760e-01,.4430e-01,.5200e-01, &
     &.6090e-01,.7130e-01,.8340e-01,.9720e-01,.1133e+00,.1317e-00, &
     &.1526e-00,.1780e-00,.2050e-00,.2370e-00,.2740e-00,.3160e-00, &
     &.3630e-00,.4170e-00,.4790e-00,.5490e-00,.6280e-00,.7180e-00, &
     &.8190e-00,.9340e-00,.1064e+01,.1209e+01,.1368e+01,.1560e+01, &
     &.1770e+01,.1990e+01,.2260e+01,.2540e+01,.2880e+01,.3230e+01, &
     &.3640e+01,.4090e+01,.4590e+01,.5140e+01,.5770e+01,.6450e+01, &
     &.7220e+01/
      data table2/.8050e+01,.8990e+01,.1001e+02,.1112e+02,.1240e+02, &
     &.1380e+02,.1530e+02,.1700e+02,.1880e+02,.2080e+02,.2310e+02, &
     &.2550e+02,.2810e+02,.3100e+02,.3420e+02,.3770e+02,.4150e+02, &
     &.4560e+02,.5010e+02,.5500e+02,.6030e+02,.6620e+02,.7240e+02, &
     &.7930e+02,.8680e+02,.9500e+02,.1146e+03,.1254e+03,.1361e+03, &
     &.1486e+03,.1602e+03,.1734e+03,.1873e+03,.2020e+03,.2171e+03, &
     &.2331e+03,.2502e+03,.2678e+03,.2863e+03,.3057e+03,.3250e+03, &
     &.3457e+03,.3664e+03,.3882e+03,.4101e+03,.4326e+03,.4584e+03, &
     &.4885e+03,.5206e+03,.5541e+03,.5898e+03,.6273e+03,.6665e+03, &
     &.7090e+03/
      data table3/.7520e+03,.7980e+03,.8470e+03,.8980e+03,.9520e+03, &
     &.1008e+04,.1067e+04,.1129e+04,.1194e+04,.1263e+04,.1334e+04, &
     &.1409e+04,.1488e+04,.1569e+04,.1656e+04,.1745e+04,.1840e+04, &
     &.1937e+04,.2041e+04,.2147e+04,.2259e+04,.2375e+04,.2497e+04, &
     &.2624e+04,.2756e+04,.2893e+04,.3036e+04,.3186e+04,.3340e+04, &
     &.3502e+04,.3670e+04,.3843e+04,.4025e+04,.4213e+04,.4408e+04, &
     &.4611e+04,.4821e+04,.5035e+04,.5270e+04,.5500e+04,.5740e+04, &
     &.6000e+04,.6250e+04,.6520e+04,.6810e+04,.7090e+04,.7390e+04, &
     &.7700e+04,.8020e+04,.8350e+04,.8690e+04,.9040e+04,.9410e+04, &
     &.9780e+04/
      data table4/.1016e+05,.1057e+05,.1098e+05,.1140e+05,.1184e+05, &
     &.1230e+05,.1275e+05,.1324e+05,.1373e+05,.1423e+05,.1476e+05, &
     &.1530e+05,.1585e+05,.1642e+05,.1700e+05,.1761e+05,.1822e+05, &
     &.1886e+05,.1950e+05,.2018e+05,.2087e+05,.2158e+05,.2229e+05, &
     &.2304e+05,.2381e+05,.2459e+05,.2539e+05,.2621e+05,.2706e+05, &
     &.2792e+05,.2881e+05,.2971e+05,.3065e+05,.3160e+05,.3257e+05, &
     &.3357e+05,.3459e+05,.3564e+05,.3669e+05,.3780e+05,.0000e+00/
      real,parameter :: cp = 1.00464e7
      real,parameter :: g = 980.6
      real,parameter :: rgas = 2.87e6
      real,parameter :: og = 1./g
      character*255 :: message
      j = IFIX(tjloc(its))
       cor1 = .120
       cor2 = 720.
       yz= -0.0001344
       y1= 3.015e-05
       y2= 1.517e-06
       y3= -3.567e-08
       y4= 2.046e-10
      do i = its,ite
        z10(i) = 1000.
        z2 (i) = 200.
        pss(i) = pspc(i)
        tstar(i) = tstrc(i)
        ukmax(i) = upc(1,i)
        vkmax(i) = vpc(1,i)
        tkmax(i) = tpc(1,i)
        rkmax(i) = rpc(1,i)
      enddo
      do i = its,ite
        windp(i) = SQRT(ukmax(i)*ukmax(i) + vkmax(i)*vkmax(i))
        wind (i) = amax1(windp(i),100.)
        if (zoc(i) .LT. amask) zoc(i) = -0.0185*0.001*wind(i)*wind(i)*og
        if (zoc(i) .GT. 0.0) then
          ecof(i) = wetc(i)
          land(i) = 1.0
          zot (i) = zoc(i)
        else
          ecof(i) = wetc(i)
          land(i) = 0.0
          zot (i) = zoc(i)
          windmks=wind(i)*.01
          znott=0.2375*exp(-0.5250*windmks) + 0.0025*exp(-0.0211*windmks)
          znott=0.01*znott
          if(windmks.le. 7.) &
          znott = (0.0185/9.8*(7.59e-8*wind(i)**2+ &
                                  2.46e-4*wind(i))**2)
          zot (i) = 100.*znott
              zoc(i)=-(0.0185/9.8*(7.59e-8*wind(i)**2+ &
                                  2.46e-4*wind(i))**2)*100.
            if(wind(i).ge.1250.0) &
             zoc(i)=-(.000739793 * wind(i) -0.58)/10
        if(wind(i).ge.3000.) then
        windmks=wind(i)*.01
        znotm = yz +windmks*y1 +windmks**2*y2 +windmks**3*y3 +windmks**4*y4
        zoc(i) = 100.*znotm
        endif
        endif
      mzoc(i) = zoc(i)
      enddo
      en = 2.
      c = .76
      a = 5.
      ca = .4
      cmo1 = .5*a - 1.648
      cmo2 = 17.193 + .5*a - 10.*c
      boycon = .61
        rovcp=rgas/cp
      do i = its,ite
        theta(i) = tkmax(i)/((pkmax(i)/pspc(i))**rovcp)
        vrtkx(i) = 1.0 + boycon*rkmax(i)
        zkmax(i) = -rgas*tkmax(i)*alog(pkmax(i)/pspc(i))*og
      enddo
      do i = its,ite
        tsg (i) = tstar(i)
        tab1 (i) = tstar(i) - 153.16
        it (i) = IFIX(tab1(i))
        tab2 (i) = tab1(i) - FLOAT(it(i))
        t1 (i) = tab(it(i) + 1)
        t2 (i) = table(it(i) + 1)
        estso(i) = t1(i) + tab2(i)*t2(i)
         psps1 = (pss(i) - estso(i))
          if(psps1 .EQ. 0.0)then
           psps1 = .1
          endif
        rstso(i) = 0.622*estso(i)/psps1
        vrts (i) = 1. + boycon*ecof(i)*rstso(i)
      enddo
      do i = its,ite
        tempa1(i) = theta(i)*vrtkx(i) - tstar(i)*vrts(i)
        tempa2(i) = tempa1(i)*(theta(i) - tstar(i))
        if (tempa2(i) .LT. 0.) tempa1(i) = 1.0e-4
        tab1(i) = ABS(tempa1(i))
        if (tab1(i) .LT. 1.0e-4) tempa1(i) = 1.0e-4
        rib (i) = g*zkmax(i)*tempa1(i)/ &
                                    (tkmax(i)*vrtkx(i)*wind(i)*wind(i))
        tab2(i) = ABS(zoc(i))
        tab1(i) = 0.95/(c*(1. - tab2(i)/zkmax(i)))
        if (rib(i) .GT. tab1(i)) rib(i) = tab1(i)
      enddo
      do i = its,ite
        zeta(i) = ca*rib(i)/0.03
      enddo
      rca = 1./ca
      enrca = en*rca
      enrca = 0.0
      zog = .0185*og
      ip = 0
      do i = its,ite
        if (zeta(i) .GE. 0.0) then
          ip = ip + 1
          istb(ip) = i
        endif
      enddo
      if (ip .EQ. 0) go to 170
      do i = 1,ip
        szetam(i) = 1.0e+30
        sgzm(i) = 0.0e+00
        szeta(i) = zeta(istb(i))
        ifz(i) = 1
      enddo
      do icnt = 1,icntx
        do i = 1,ip
          if (ifz(i) .EQ. 0) go to 80
          zal1g = ALOG(szeta(i))
          if (szeta(i) .LE. 0.5) then
            fmz1(i) = (zal1g + a*szeta(i))*rca
          else if (szeta(i) .GT. 0.5 .AND. szeta(i) .LE. 10.) then
            rzeta1 = 1./szeta(i)
            fmz1(i) = (8.*zal1g + 4.25*rzeta1 - &
                                          0.5*rzeta1*rzeta1 + cmo1)*rca
          else if (szeta(i) .GT. 10.) then
            fmz1(i) = (c*szeta(i) + cmo2)*rca
          endif
          szetao = ABS(zoc(istb(i)))/zkmax(istb(i))*szeta(i)
          zal2g = ALOG(szetao)
          if (szetao .LE. 0.5) then
            fmzo1(i) = (zal2g + a*szetao)*rca
            sfzo (i) = 1. + a*szetao
          else if (szetao .GT. 0.5 .AND. szetao .LE. 10.) then
            rzeta2 = 1./szetao
            fmzo1(i) = (8.*zal2g + 4.25*rzeta2 - &
                                          0.5*rzeta2*rzeta2 + cmo1)*rca
            sfzo (i) = 8.0 - 4.25*rzeta2 + rzeta2*rzeta2
          else if (szetao .GT. 10.) then
            fmzo1(i) = (c*szetao + cmo2)*rca
            sfzo (i) = c*szetao
          endif
          szetho = ABS(zot(istb(i)))/zkmax(istb(i))*szeta(i)
          zal2gh = ALOG(szetho)
          if (szetho .LE. 0.5) then
            fhzo1(i) = (zal2gh + a*szetho)*rca
            sfzo (i) = 1. + a*szetho
          else if (szetho .GT. 0.5 .AND. szetho .LE. 10.) then
            rzeta2 = 1./szetho
            fhzo1(i) = (8.*zal2gh + 4.25*rzeta2 - &
                                          0.5*rzeta2*rzeta2 + cmo1)*rca
            sfzo (i) = 8.0 - 4.25*rzeta2 + rzeta2*rzeta2
          else if (szetho .GT. 10.) then
            fhzo1(i) = (c*szetho + cmo2)*rca
            sfzo (i) = c*szetho
          endif
            szet10 = ABS(z10(istb(i)))/zkmax(istb(i))*szeta(i)
            zal10 = ALOG(szet10)
            if (szet10 .LE. 0.5) then
              fmz10(i) = (zal10 + a*szet10)*rca
            else if (szet10 .GT. 0.5 .AND. szet10 .LE. 10.) then
              rzeta2 = 1./szet10
              fmz10(i) = (8.*zal10 + 4.25*rzeta2 - &
                                          0.5*rzeta2*rzeta2 + cmo1)*rca
            else if (szet10 .GT. 10.) then
              fmz10(i) = (c*szet10 + cmo2)*rca
            endif
            sf10(i) = fmz10(i) - fmzo1(i)
            szet2 = ABS(z2 (istb(i)))/zkmax(istb(i))*szeta(i)
            zal2 = ALOG(szet2 )
            if (szet2 .LE. 0.5) then
              fmz2 (i) = (zal2 + a*szet2 )*rca
            else if (szet2 .GT. 0.5 .AND. szet2 .LE. 2.) then
              rzeta2 = 1./szet2
              fmz2 (i) = (8.*zal2 + 4.25*rzeta2 - &
                                          0.5*rzeta2*rzeta2 + cmo1)*rca
            else if (szet2 .GT. 2.) then
              fmz2 (i) = (c*szet2 + cmo2)*rca
            endif
            sf2 (i) = fmz2 (i) - fmzo1(i)
          sfm(i) = fmz1(i) - fmzo1(i)
          sfh(i) = fmz1(i) - fhzo1(i)
          sgz = ca*rib(istb(i))*sfm(i)*sfm(i)/ &
                                               (sfh(i) + enrca*sfzo(i))
          fmz = (sgz - szeta(i))/szeta(i)
          fmzo = ABS(fmz)
          if (fmzo .GE. 5.0e-5) then
            sq = (sgz - sgzm(i))/(szeta(i) - szetam(i))
            if(sq .EQ. 1) then
             write(message,*)'NCO ERROR DIVIDE BY ZERO IN MFLUX2 (STABLE CASE)'
             call wrf_message(message)
             write(message,*)'sq is 1 ',fmzo,sgz,sgzm(i),szeta(i),szetam(i)
             call wrf_message(message)
            endif
            szetam(i) = szeta(i)
            szeta (i) = (sgz - szeta(i)*sq)/(1.0 - sq)
            sgzm (i) = sgz
          else
            ifz(i) = 0
          endif
80 continue
        enddo
      enddo
      do i = 1,ip
        if (ifz(i) .GE. 1) go to 110
      enddo
      go to 130
110 continue
      write(6,120)
120 format(2X, ' NON-CONVERGENCE FOR STABLE ZETA IN ROW ')
130 continue
      do i = 1,ip
        szo = zoc(istb(i))
        if (szo .LT. 0.0) then
        wndm=wind(istb(i))*0.01
          if(wndm.lt.15.0) then
          ckg=0.0185*og
          else
           ckg=(sfenth*(4*0.000308*wndm) + (1.-sfenth)*0.0185 )*og
          endif
       szo = - ckg*wind(istb(i))*wind(istb(i))/ &
                             (sfm(i)*sfm(i))
        cons_p000001 = .000001
        cons_7 = 7.
        vis = 1.4E-1
        ustar = sqrt( -szo / zog)
        restar = -ustar * szo / vis
        restar = max(restar,cons_p000001)
        rat = 2.67 * restar ** .25 - 2.57
        rat = min(rat ,cons_7)
 rat=0.
        zot(istb(i)) = szo * exp(-rat)
         else
         zot(istb(i)) = zoc(istb(i))
        endif
             zoc(istb(i)) = szo
      enddo
      do i = 1,ip
        xxfm(istb(i)) = sfm(i)
        xxfh(istb(i)) = sfh(i)
        xxfh2(istb(i)) = sf2 (i)
        xxsh(istb(i)) = sfzo(i)
      enddo
        do i = 1,ip
         wind10(istb(i)) = sf10(i)*wind(istb(i))/sfm(i)
         wind10(istb(i)) = wind10(istb(i)) * 1.944
           if(wind10(istb(i)) .GT. 6000.0) then
       wind10(istb(i))=wind10(istb(i))+wind10(istb(i))*cor1 &
                        - cor2
           endif
         wind10(istb(i)) = wind10(istb(i)) / 1.944
        enddo
170 continue
      iq = 0
      do i = its,ite
        if (zeta(i) .LT. 0.0) then
          iq = iq + 1
          iutb(iq) = i
        endif
      enddo
      if (iq .EQ. 0) go to 290
      do i = 1,iq
        uzeta (i) = zeta(iutb(i))
        ifz (i) = 1
        uzetam(i) = 1.0e+30
        ugzm (i) = 0.0e+00
      enddo
      do icnt = 1,icntx
        do i = 1,iq
          if (ifz(i) .EQ. 0) go to 200
          ugzzo = ALOG(zkmax(iutb(i))/ABS(zot(iutb(i))))
          uzetao = ABS(zot(iutb(i)))/zkmax(iutb(i))*uzeta(i)
          ux11 = 1. - 16.*uzeta(i)
          ux12 = 1. - 16.*uzetao
          y = SQRT(ux11)
          yo = SQRT(ux12)
          ufzo(i) = 1./yo
          ux13 = (1. + y)/(1. + yo)
          ux21 = ALOG(ux13)
          ufh(i) = (ugzzo - 2.*ux21)*rca
          ugzzo = ALOG(zkmax(iutb(i))/ABS(zoc(iutb(i))))
          uzetao = ABS(zoc(iutb(i)))/zkmax(iutb(i))*uzeta(i)
          ux11 = 1. - 16.*uzeta(i)
          ux12 = 1. - 16.*uzetao
          y = SQRT(ux11)
          yo = SQRT(ux12)
          ux13 = (1. + y)/(1. + yo)
          ux21 = ALOG(ux13)
          x = SQRT(y)
          xo = SQRT(yo)
          xnum = (x**2 + 1.)*((x + 1.)**2)
          xden = (xo**2 + 1.)*((xo + 1.)**2)
          xtan = ATAN(x) - ATAN(xo)
          ux3 = ALOG(xnum/xden)
          ufm(i) = (ugzzo - ux3 + 2.*xtan)*rca
            ugz10 = ALOG(z10(iutb(i))/ABS(zoc(iutb(i))))
            uzet1o = ABS(z10(iutb(i)))/zkmax(iutb(i))*uzeta(i)
            uzetao = ABS(zoc(iutb(i)))/zkmax(iutb(i))*uzeta(i)
            ux11 = 1. - 16.*uzet1o
            ux12 = 1. - 16.*uzetao
            y = SQRT(ux11)
            y10 = SQRT(ux12)
            ux13 = (1. + y)/(1. + y10)
            ux21 = ALOG(ux13)
            x = SQRT(y)
            x10 = SQRT(y10)
            xnum = (x**2 + 1.)*((x + 1.)**2)
            xden = (x10**2 + 1.)*((x10 + 1.)**2)
            xtan = ATAN(x) - ATAN(x10)
            ux3 = ALOG(xnum/xden)
            uf10(i) = (ugz10 - ux3 + 2.*xtan)*rca
          ugz2 = ALOG(z2 (iutb(i))/ABS(zoc(iutb(i))))
          uzet1o = ABS(z2 (iutb(i)))/zkmax(iutb(i))*uzeta(i)
          uzetao = ABS(zoc(iutb(i)))/zkmax(iutb(i))*uzeta(i)
          ux11 = 1. - 16.*uzet1o
          ux12 = 1. - 16.*uzetao
          y = SQRT(ux11)
          yo = SQRT(ux12)
          ux13 = (1. + y)/(1. + yo)
          ux21 = ALOG(ux13)
          uf2 (i) = (ugzzo - 2.*ux21)*rca
          ugz = ca*rib(iutb(i))*ufm(i)*ufm(i)/(ufh(i) + enrca*ufzo(i))
          ux1 = (ugz - uzeta(i))/uzeta(i)
          ux2 = ABS(ux1)
          if (ux2 .GE. 5.0e-5) then
            uq = (ugz - ugzm(i))/(uzeta(i) - uzetam(i))
            uzetam(i) = uzeta(i)
            if(uq .EQ. 1) then
             write(0,*)'NCO ERROR DIVIDE BY ZERO IN MFLUX2 (UNSTABLE CASE)'
             write(0,*)'uq is 1 ',ux2,ugz,ugzm(i),uzeta(i),uzetam(i)
            endif
            uzeta (i) = (ugz - uzeta(i)*uq)/(1.0 - uq)
            ugzm (i) = ugz
          else
            ifz(i) = 0
          endif
200 continue
        enddo
      enddo
      do i = 1,iq
        if (ifz(i) .GE. 1) go to 230
      enddo
      go to 250
230 continue
      write(6,240)
240 format(2X, ' NON-CONVERGENCE FOR UNSTABLE ZETA IN ROW ')
250 continue
      do i = 1,iq
        uzo = zoc(iutb(i))
        if (zoc(iutb(i)) .LT. 0.0) then
        wndm=wind(iutb(i))*0.01
         if(wndm.lt.15.0) then
         ckg=0.0185*og
         else
            ckg=(4*0.000308*wndm)*og
            ckg=(sfenth*(4*0.000308*wndm) + (1.-sfenth)*0.0185 )*og
         endif
        uzo =-ckg*wind(iutb(i))*wind(iutb(i))/(ufm(i)*ufm(i))
        cons_p000001 = .000001
        cons_7 = 7.
        vis = 1.4E-1
        ustar = sqrt( -uzo / zog)
        restar = -ustar * uzo / vis
        restar = max(restar,cons_p000001)
        rat = 2.67 * restar ** .25 - 2.57
        rat = min(rat ,cons_7)
 rat=0.0
        zot(iutb(i)) = uzo * exp(-rat)
         else
         zot(iutb(i)) = zoc(iutb(i))
        endif
            zoc(iutb(i)) = uzo
      enddo
        do i = 1,iq
          wind10(iutb(i)) = uf10(i)*wind(iutb(i))/ufm(i)
          wind10(iutb(i)) = wind10(iutb(i)) * 1.944
           if(wind10(iutb(i)) .GT. 6000.0) then
         wind10(iutb(i))=wind10(iutb(i))+wind10(iutb(i))*cor1 &
                         - cor2
           endif
          wind10(iutb(i)) = wind10(iutb(i)) / 1.944
        enddo
      do i = 1,iq
        xxfm(iutb(i)) = ufm(i)
        xxfh(iutb(i)) = ufh(i)
        xxfh2(iutb(i)) = uf2 (i)
        xxsh(iutb(i)) = ufzo(i)
      enddo
290 continue
      do i = its,ite
        ucom(i) = ukmax(i)
        vcom(i) = vkmax(i)
        if (windp(i) .EQ. 0.0) then
          windp(i) = 100.0
          ucom (i) = 100.0/SQRT(2.0)
          vcom (i) = 100.0/SQRT(2.0)
        endif
        rho(i) = pss(i)/(rgas*(tsg(i) + enrca*(theta(i) - &
                tsg(i))*xxsh(i)/(xxfh(i) + enrca*xxsh(i))))
        bq1(i) = wind(i)*rho(i)/(xxfm(i)*(xxfh(i) + enrca*xxsh(i)))
      enddo
      if (ntsflg .EQ. 0) go to 370
      alll = 600.
      xks = 0.01
      hcap = .5/2.39e-8
      pith = SQRT(4.*ATAN(1.0))
      alfus = alll/2.39e-8
      teps = 0.1
      ip = 0
      do i = its,ite
        if (land(i) .EQ. 1) then
          ip = ip + 1
          indx (ip) = i
          slwa (ip) = slwdc(i)/(2.39e-8*60.)
          tss (ip) = tstar(i)
          thetap (ip) = theta(i)
          rkmaxp (ip) = rkmax(i)
          aap (ip) = 5.673e-5
          pssp (ip) = pss(i)
          ecofp (ip) = ecof(i)
          estsop (ip) = estso(i)
          rstsop (ip) = rstso(i)
          bq1p (ip) = bq1(i)
          bq1p (ip) = amax1(bq1p(ip),0.1e-3)
          delsrad(ip) = dt *pith/(hcap*SQRT(3600.*24.*xks))
        endif
      enddo
      do i = 1,ip
        ifz (i) = 1
        tsm (i) = tss(i)
        rdiff(i) = amin1(0.0,(rkmaxp(i) - rstsop(i)))
300 format(2X, ' SURFACE EQUILIBRIUM CALCULATION ')
          foftm(i) = tss(i) + delsrad(i)*(slwa(i) - aap(i)*tsm(i)**4 - &
           cp*bq1p(i)*(tsm(i) - thetap(i)) + ecofp(i)*alfus*bq1p(i)* &
           rdiff(i))
        tsp(i) = foftm(i)
      enddo
      do icnt = 1,icntx
        do i = 1,ip
          if (ifz(i) .EQ. 0) go to 330
          tab1 (i) = tsp(i) - 153.16
          it (i) = IFIX(tab1(i))
          tab2 (i) = tab1(i) - FLOAT(it(i))
          t1 (i) = tab(it(i) + 1)
          t2 (i) = table(it(i) + 1)
          estsop(i) = t1(i) + tab2(i)*t2(i)
             psps2 = (pssp(i) - estsop(i))
             if(psps2 .EQ. 0.0)then
               psps2 = .1
             endif
          rstsop(i) = 0.622*estsop(i)/psps2
          rdiff (i) = amin1(0.0,(rkmaxp(i) - rstsop(i)))
            foft(i) = tss(i) + delsrad(i)*(slwa(i) - aap(i)*tsp(i)**4 - &
             cp*bq1p(i)*(tsp(i) - thetap(i)) + ecofp(i)*alfus*bq1p(i)* &
             rdiff(i))
          frac(i) = ABS((foft(i) - tsp(i))/tsp(i))
          if (frac(i) .GE. teps) then
            qf (i) = (foft(i) - foftm(i))/(tsp(i) - tsm(i))
            tsm (i) = tsp(i)
            tsp (i) = (foft(i) - tsp(i)*qf(i))/(1. - qf(i))
            foftm(i) = foft(i)
          else
            ifz(i) = 0
          endif
330 continue
        enddo
      enddo
      do i = 1,ip
        if (ifz(i) .EQ. 1) then
        write(6, 340) tsp(i), i, j
340 format(2X, ' NON-CONVERGENCE OF T* PREDICTED (T*,I,J) = ', E14.8, &
            2I5)
        write(6,345) indx(i), j, tstar(indx(i)), tsp(i), ip
345 format(2X, ' I, J, OLD T*, NEW T*, NPTS ', 2I5, 2E14.8, I5)
350 format(2X, ' REFLECT, SIGT4, SHFX, ALEVP, DELTEN, DIFFOT ', &
            6E14.8)
        endif
      enddo
      do i = 1,ip
        ii = indx(i)
        tstrc(ii) = tsp (i)
      enddo
370 continue
      do i = its,ite
        fxh(i) = bq1(i)*(theta(i) - tsg(i))
        fxe(i) = ecof(i)*bq1(i)*(rkmax(i) - rstso(i))
        if (fxe(i) .GT. 0.0) fxe(i) = 0.0
        fxmx(i) = rho(i)/(xxfm(i)*xxfm(i))*wind(i)*wind(i)*ucom(i)/ &
                 windp(i)
        fxmy(i) = rho(i)/(xxfm(i)*xxfm(i))*wind(i)*wind(i)*vcom(i)/ &
        windp(i)
        cdm(i) = 1./(xxfm(i)*xxfm(i))
      enddo
      return
      end subroutine MFLUX2
  SUBROUTINE hwrfsfcinit(isn,XICE,VEGFRA,SNOW,SNOWC,CANWAT,SMSTAV, &
                        SMSTOT, SFCRUNOFF,UDRUNOFF,GRDFLX,ACSNOW, &
                        ACSNOM,IVGTYP,ISLTYP,TSLB,SMOIS,DZS,SFCEVP, &
                        TMN, &
                        num_soil_layers, &
                        allowed_to_read, &
                        ids,ide, jds,jde, kds,kde, &
                        ims,ime, jms,jme, kms,kme, &
                        its,ite, jts,jte, kts,kte )
   IMPLICIT NONE
   INTEGER, INTENT(IN ) :: ids,ide, jds,jde, kds,kde, &
                                    ims,ime, jms,jme, kms,kme, &
                                    its,ite, jts,jte, kts,kte
   INTEGER, INTENT(IN) :: num_soil_layers
   REAL, DIMENSION( num_soil_layers), INTENT(IN) :: DZS
   REAL, DIMENSION( ims:ime, num_soil_layers, jms:jme ) , &
            INTENT(INOUT) :: SMOIS, &
                                                         TSLB
   REAL, DIMENSION( ims:ime, jms:jme ) , &
            INTENT(INOUT) :: SNOW, &
                                                         SNOWC, &
                                                        CANWAT, &
                                                        SMSTAV, &
                                                        SMSTOT, &
                                                     SFCRUNOFF, &
                                                      UDRUNOFF, &
                                                        SFCEVP, &
                                                        GRDFLX, &
                                                        ACSNOW, &
                                                          XICE, &
                                                        VEGFRA, &
                                                        TMN, &
                                                        ACSNOM
   INTEGER, DIMENSION( ims:ime, jms:jme ) , &
            INTENT(INOUT) :: IVGTYP, &
                                                        ISLTYP
  INTEGER, INTENT(IN) :: isn
  LOGICAL, INTENT(IN) :: allowed_to_read
  INTEGER :: iseason
  INTEGER :: icm,jcm,itf,jtf
  INTEGER :: I,J,L
   itf=min0(ite,ide-1)
   jtf=min0(jte,jde-1)
   icm = ide/2
   jcm = jde/2
   iseason=isn
   DO J=jts,jtf
       DO I=its,itf
       SNOWC(i,j)=0.
    ENDDO
   ENDDO
  END SUBROUTINE hwrfsfcinit
END MODULE module_sf_gfdl


text:
great(aa1,p1).
great(a1,j1).
great(c1,bb1).
great(jj1,ll1).
great(q1,bb1).
great(hh1,m1).
great(hh1,j1).
great(hh1,x1).
great(d1,g1).
great(ff1,s1).
great(c1,j1).
great(l1,aa1).
great(jj1,y1).
great(i1,x1).
great(c1,a1).
great(a1,o1).
great(j1,aa1).
great(h1,l1).
great(ii1,ff1).
great(bb1,x1).
great(h1,t1).
great(cc1,dd1).
great(ff1,u1).
great(ii1,s1).
great(v1,ee1).
great(g1,b1).

max_stars_repo_name: manoelfranca/cilppp

max_issues_repo_name: manoelfranca/cilppp

max_forks_repo_name: manoelfranca/cilppp