MODULE SOCPBL 10,1
!@sum  SOCPBL deals with boundary layer physics
!@auth Ye Cheng/G. Hartke (modifications by G. Schmidt)
!@ver  1.0 (from PBLB336E)
!@cont pbl,advanc,stars,getl,dflux,simil,griddr,tfix
!@cont ccoeff0,getk,e_eqn,t_eqn,q_eqn,uv_eqn
!@cont t_eqn_sta,q_eqn_sta,uv_eqn_sta
!@cont inits,tcheck,ucheck,check1,output,rtsafe

      USE CONSTANT, only : grav,omega,pi,radian,bygrav,teeny,deltx,tf
     &                     ,by3,shv,sha,shw,lhe,stbo,rhow,rgas

      IMPLICIT NONE

      integer, parameter :: n=8  !@param n  no of pbl. layers

!@var ZS1    = height of the first model layer (m)
!@var TGV    = virtual potential temperature of the ground (K)
!@var TKV    = virtual potential temperature of first model layer (K)
!@var HEMI   = 1 for northern hemisphere, -1 for southern hemisphere
!@var POLE   = .TRUE. if at the north or south pole, .FALSE. otherwise

!@var US     = x component of surface wind, positive eastward (m/s)
!@var VS     = y component of surface wind, positive northward (m/s)
!@var WS     = magnitude of the surface wind (m/s)
!@var WSM    = magnitude of the surface wind - ocean currents (m/s)
!@var WSH    = magnitude of surface wind modified by buoyancy flux(m/s)
!@var TSV    = virtual potential temperature of the surface (K)
!@var QS     = surface value of the specific moisture
!@var PSI    = angular diff. btw geostrophic and surface winds (rads)
!@var DBL    = boundary layer height (m)
!@var KMS    = momentum transport coefficient at ZGS (m**2/s)
!@var KHS    = heat transport coefficient at ZGS (m**2/s)
!@var KHQ    = moist transport coefficient at ZGS (m**2/s)
!@var PPBL   = pressure at DBL (mb)
!@var USTAR  = friction speed (square root of momentum flux) (m/s)
!@var CM     = drag coefficient (dimensionless surface momentum flux)
!@var CH     = Stanton number   (dimensionless surface heat flux)
!@var CQ     = Dalton number    (dimensionless surface moisture flux)
!@var z0m   = roughness length for momentum,
!@+           prescribed for itype=3,4 but computed for itype=1,2 (m)
!@var z0h   = roughness length for temperature (m)
!@var z0q   = roughness length for water vapor (m)
!@var UG     = eastward component of the geostrophic wind (m/s)
!@var VG     = northward component of the geostrophic wind (m/s)
!@var WG     = magnitude of the geostrophic wind (m/s)
!@var MDF    = downdraft mass flux (m/s)
!@var WINT   = integrated surface wind speed over sgs wind distribution
      real*8 :: zs1,tgv,tkv,qg_sat,qg_aver,hemi,dtsurf,w2_1,mdf,wint
      real*8 :: us,vs,ws,wsm,wsh,tsv,qsrf,psi,dbl,kms,khs,kqs,ppbl
     *         ,ustar,cm,ch,cq,z0m,z0h,z0q,ug,vg,wg,XCDpbl=1d0
      logical :: pole

      real*8 ::  dpdxr,dpdyr,dpdxr0,dpdyr0
      real*8 :: rimax,ghmin,ghmax,gmmax0,d1,d2,d3,d4,d5
     *     ,s0,s1,s2,s4,s5,s6,c1,c2,c3,c4,c5,b1,b123,b2,prt

      !for level 3 model only:
      real*8 :: g0,d1_3,d2_3,d3_3,d4_3,d5_3
     *         ,s0_3,s1_3,s2_3,s3_3,s4_3,s5_3,s6_3
     *         ,g1,g2,g3,g4,g5,g6,g7,g8

C**** boundary layer parameters
      real*8, parameter :: kappa=0.40d0 !@var kappa  Von Karman constant
      real*8, parameter :: zgs=10. !@var zgs height of surface layer (m)



C**** parameters for surface fluxes
      !Hogstrom 1988:
      real*8, parameter :: sigma=0.95d0,sigma1=1.-sigma
      real*8, parameter :: gamamu=19.3d0,gamahu=11.6d0,gamams=6.d0,
     *     gamahs=7.8d0/sigma

      ! Businger 1971:
ccc   real*8, parameter :: sigma=0.74d0,sigma1=1.-sigma
ccc   real*8, parameter :: gamamu=15.0d0,gamahu=9.d0,gamams=4.7d0,
ccc  *     gamahs=4.7d0/sigma

C***
C***  Thread-Private Common
C***
      COMMON /PBLTPC/ dpdxr,dpdyr,dpdxr0,dpdyr0

!$OMP  THREADPRIVATE (/PBLTPC/)

      COMMON /PBLPAR/ZS1,TGV,TKV,QG_SAT,qg_aver,HEMI,POLE
      COMMON /PBLOUT/US,VS,WS,WSM,WSH,TSV,QSRF,PSI,DBL,KMS,KHS,KQS,PPBL,
     *     USTAR,CM,CH,CQ,Z0M,Z0H,Z0Q,UG,VG,WG,W2_1,MDF,WINT
!$OMP  THREADPRIVATE (/PBLPAR/,/PBLOUT/)

CCC !@var bgrid log-linear gridding parameter
CCC      real*8 :: bgrid

!@var smax,smin,cmax,cmin limits on drag coeffs.
!@var emax limit on turbulent kinetic energy
      real*8, parameter :: smax=0.25d0,smin=0.005d0,cmax=smax*smax,
     *     cmin=smin*smin,emax=1.d5

!@param twoby3 2/3 constant
      real*8, parameter :: twoby3 = 2d0/3d0

      CONTAINS


      subroutine advanc( 1,12
     3     coriol,utop,vtop,ttop,qtop,tgrnd,
     &     qgrnd,qgrnd_sat,evap_max,fr_sat,

     4     psurf,trhr0,ztop,dtime,ufluxs,vfluxs,tfluxs,qfluxs,
     5     uocean,vocean,ts_guess,ilong,jlat,itype)
!@sum  advanc  time steps the solutions for the boundary layer variables
!@auth  Ye Cheng/G. Hartke
!@ver   1.0
c    input:
!@var  z0m  roughness height for momentum (if itype>2)
!@var  zgs  height of the surface layer (nominally 10 m)
!@var  coriol  2.*omega*sin(latitude), the coriolis factor
!@var  utop  x component of wind at the top of the layer
!@var  vtop  y component of wind at the top of the layer
!@var  ttop  virtual potential temperature at the top of the layer
!@var  qtop  moisture at the top of the layer
!@var  tgrnd  virt. pot. temp. of the ground, at the roughness height
!@var  qgrnd  moisture at the ground, at the roughness height
!@var  qgrnd_sat saturated moisture at the ground, at the roughness height
!@var  evap_max maximal evaporation from unsaturated soil
!@var  fr_sat fraction of saturated soil
!@var  ztop height of the first model layer, approx 200 m if lm=9
!@var  dtime  time step
!@var  ilong  longitude identifier
!@var  jlat   latitude identifier
!@var  itype  1, ocean; 2, ocean ice; 3, land ice; 4, land
c   output:
!@var  us  x component of the surface wind (i.e., due east)
!@var  vs  y component of the surface wind (i.e., due north)
!@var  tsv  virtual potential surface temperature
!@var  qsrf  surface specific moisture
!@var  kms  surface value of km
!@var  khs  surface value of kh
!@var  kqs  surface value of kq
!@var  wsh  magnitude of surface wind modified by buoyancy flux (m/s)
!@var  ustar  friction speed
!@var  cm  dimensionless momentum flux at surface (drag coeff.)
!@var  ch  dimensionless heat flux at surface (stanton number)
!@var  cq  dimensionless moisture flux at surface (dalton number)
!@var  z0m  roughness height for momentum (if itype=1 or 2)
!@var  z0h  roughness height for heat
!@var  z0q  roughness height for moisture
!@var  psurf surface pressure
!@var  trhr0 incident long wave radiation 


c  internals:
!@var  n     number of the local, vertical grid points
!@var  lscale turbulence length scale. computed on secondary grid.
!@var  z     altitude of primary vertical grid points
!@var  zhat  altitude of secondary vertical grid points
!@var  dzh   dz evaluated at zhat(i)
!@var  dz    dz evaluated at z(i)
!@var  dxi   (z(n)-z(1))/(n-1)
!@var  km    turbulent momentum tranport coefficient.
!@var  kh    turbulent thermometric conductivity. computed
!@var  ke    transport coefficient for the turbulent kinetic energy.
!@var  ipbl  stores bl properties of last time step
      implicit none

      real*8, intent(in) :: coriol,utop,vtop,ttop,qtop,uocean,vocean
     *                      ,ts_guess
      real*8, intent(in) :: tgrnd,qgrnd,qgrnd_sat,evap_max,fr_sat
     &     ,ztop,dtime
      real*8, intent(out) :: ufluxs,vfluxs,tfluxs,qfluxs
      integer, intent(in) :: ilong,jlat,itype
      real*8, intent(in) :: psurf,trhr0

      real*8 :: lmonin,tstar,qstar,ustar0,test,wstar3,wstar3fac,wstar2h
      real*8 :: bgrid,an2,as2,dudz,dvdz,tau
      real*8, parameter ::  tol=1d-3,w=.5d0
      integer, parameter ::  itmax=50
      integer, parameter :: iprint=0,jprint=41  ! set iprint>0 to debug
      real*8, dimension(n) :: z,dz,xi,usave,vsave,tsave,qsave
     *       ,usave1,vsave1,tsave1,qsave1             
      real*8, dimension(n-1) :: lscale,zhat,dzh,xihat,km,kh,kq,ke,gm,gh
     *     ,esave,esave1
      integer :: i,j,iter,ierr  !@var i,j,iter loop variable
C****
      real*8 :: sig0,delt,wt,wmin,wmax,sig
      integer :: icase

!@var  u  local due east component of wind
!@var  v  local due north component of wind
!@var  t  local virtual potential temperature
!@var  q  local specific humidity (a passive scalar)
!@var  e  local turbulent kinetic energy
c**** special threadprivate common block (compaq compiler stupidity)
      real*8, dimension(n) :: u,v,t,q
      real*8, dimension(n-1) :: e
      common/pbluvtq/u,v,t,q,e
!$OMP  THREADPRIVATE (/pbluvtq/)
C**** end special threadprivate common block

      call griddr(z,zhat,xi,xihat,dz,dzh,zgs,ztop,bgrid,n,ierr)
      if (ierr.gt.0) then
        print*,"In advanc: i,j,itype =",ilong,jlat,itype,us,vs,tsv,qsrf
        call abort
        call stop_model("PBL error in advanc",255)
      end if

      usave(:)=u(:)
      vsave(:)=v(:)
      tsave(:)=t(:)
      qsave(:)=q(:)
      esave(:)=e(:)

      do i=1,n-1
        usave1(i)=usave(i)
        vsave1(i)=vsave(i)
        tsave1(i)=tsave(i)
        qsave1(i)=qsave(i)
        esave1(i)=esave(i)
      end do
      ustar0=0.

      do iter=1,itmax

        if(iter.eq.1) then
          call getl1(e,zhat,dzh,lscale,n)
        else
          call getl(e,u,v,t,zhat,dzh,lmonin,ustar,lscale,dbl,n)
        endif

        call getk(km,kh,kq,ke,gm,gh,u,v,t,e,lscale,dzh,n)
        call stars(ustar,tstar,qstar,lmonin,tgrnd,qgrnd,
     2             u,v,t,q,z,z0m,z0h,z0q,cm,ch,cq,
     3             km,kh,kq,dzh,itype,n)


        call e_eqn(esave,e,u,v,t,km,kh,ke,lscale,dz,dzh,
     2                 ustar,dtime,n)

        !@var wstar the convection-induced wind according to
        !@+ M.J.Miller et al. 1992, J. Climate, 5(5), 418-434, Eqs(6-7),
        !@+ for heat and mositure
        if(          wstar3fac=-dbl*grav*2.*(          wstar2h = (wstar3fac*kh(1))**twoby3
        else
          wstar2h = 0.
        endif

        wsh = sqrt((u(1)-uocean)**2+(
        call q_eqn(qsave,q,kq,dz,dzh,cq,wsh,qgrnd_sat,qtop,dtime,n
     &       ,evap_max,fr_sat)

        call t_eqn(u,v,tsave,t,q,z,kh,kq,dz,dzh,ch,wsh,tgrnd
     &            ,ttop,dtime,n)

        call uv_eqn(usave,vsave,u,v,z,km,dz,dzh,
     2              ustar,cm,z0m,utop,vtop,dtime,coriol,
     3              ug,vg,uocean,vocean,n)

        if ((ttop.gt.tgrnd).and.(lmonin.lt.0.)) then
          call tfix(t,z,ttop,tgrnd,lmonin,tstar,ustar,kh(1),
     *              ts_guess,n)
        endif

        test=abs(2.*(ustar-ustar0)/(ustar+ustar0))
        if (test.lt.tol) exit

        if (iter.lt.itmax) then
        do i=1,n-1
          u(i)=w*usave1(i)+(1.-w)*u(i)
          v(i)=w*vsave1(i)+(1.-w)*v(i)
          t(i)=w*tsave1(i)+(1.-w)*t(i)
          q(i)=w*qsave1(i)+(1.-w)*q(i)
          e(i)=w*esave1(i)+(1.-w)*e(i)
          usave1(i)=u(i)
          vsave1(i)=v(i)
          tsave1(i)=t(i)
          qsave1(i)=q(i)
          esave1(i)=e(i)
        end do
        end if
        ustar0=ustar

      end do
c     write(97,*) "iter=",iter,ustar

c**** cannot update wsh without taking care that wsh used for tracers is
c**** the same as that used for q
c      wsh = sqrt((u(1)-uocean)**2+(v(1)-vocean)**2+wstar2h)
      wsm = wsh

C**** Preliminary coding for use of sub-gridscale wind distribution
C**** generic calculations for all tracers
C**** To use, uncomment next two lines and adapt the next chunk for
C**** your tracers. The integrated wind value is passed back to SURFACE
C**** and GHY_DRV. This may need to be tracer dependent?
csgs      delt = t(1)/(1.+q(1)*deltx) - tgrnd/(1.+qgrnd*deltx)
csgs      sig0 = sig(e(1),mdf,dbl,delt,ch,wsh,t(1))
csgsC**** possibly tracer specific coding
csgs      wt = 3.                 ! threshold velocity
csgs      wmin = wt               ! minimum wind velocity (usually wt?)
csgs      wmax = 50.              ! maxmimum wind velocity 
csgs      icase=3                 ! icase=3 ==> w^3 dependency  
csgs      call integrate_sgswind(sig0,wt,wmin,wmax,wsh,icase,wint)



      us  = u(1)
      vs  = v(1)
      tsv = t(1)
      qsrf  = q(1)
      kms = km(1)
      khs = kh(1)
      kqs = kq(1)

      ufluxs=km(1)*(      vfluxs=km(1)*(      tfluxs=kh(1)*(      qfluxs=kq(1)*(
      an2=2.*grav*(      dudz=(      dvdz=(      as2=dudz*dudz+dvdz*dvdz
      tau=B1*lscale(n-1)/max(sqrt(2.*e(n-1)),teeny)
      !@var w2_1 the vertical component of e at GCM layer 1
      w2_1=twoby3*e(n-1)-tau*by3
     &     *((3*g3-g2)*km(n-1)*as2+4.*g4*kh(n-1)*an2)
      w2_1=max(0.24d0*e(n-1),w2_1) ! Mellor-Yamada 1982

c     call check1(ustar,1,ilong,jlat,2)

c Diagnostics printed at a selected point:

      if ((ilong.eq.iprint).and.(jlat.eq.jprint)) then
        call output(u,v,t,q,e,lscale,z,zhat,dzh,
     2              km,kh,kq,ke,gm,gh,cm,ch,cq,z0m,z0h,z0q,
     3              ustar,tstar,qstar,lmonin,tgrnd,qgrnd,
     4              utop,vtop,ttop,qtop,
     5              dtime,bgrid,ilong,jlat,iter,itype,n)
      endif

      return
      end subroutine advanc


      subroutine stars(ustar,tstar,qstar,lmonin,tgrnd,qgrnd, 2,2
     2                 u,v,t,q,z,z0m,z0h,z0q,cm,ch,cq,
     3                 km,kh,kq,dzh,itype,n)
!@sum computes USTAR,TSTAR and QSTAR
!@+   Momentum flux = USTAR*USTAR
!@+   Heat flux     = USTAR*TSTAR
!@+   MOISTURE flux = USTAR*QSTAR
!@auth Ye Cheng/G. Hartke (modifications by G. Schmidt)
!@ver  1.0 (from PBLB336E)
!@var USTAR the friction speed
!@var TSTAR the virtual potential temperature scale
!@var QSTAR the moisture scale
!@var LMONIN the Monin-Obukhov length scale
      USE CONSTANT, only : teeny
      implicit none

      integer, intent(in) :: itype,n
      real*8, dimension(n), intent(in) :: u,v,t,q,z
      real*8, dimension(n-1), intent(in) :: km,kh,kq,dzh
      real*8, intent(in) :: tgrnd,qgrnd
      real*8, intent(inout) :: z0m
      real*8, intent(out) :: ustar,tstar,qstar,lmonin
      real*8, intent(out) :: z0h,z0q,cm,ch,cq

      real*8 dz,vel1,du1,dv1,dudz,dtdz,dqdz,zgs

      dz     = dzh(1)
      vel1   = sqrt(      du1=u(2)-u(1)
      dv1=v(2)-v(1)
      dudz=sqrt(du1*du1+dv1*dv1)/dz
      dtdz   = (      dqdz   = (      ustar  = sqrt(km(1)*dudz)
      ustar  = max(ustar,teeny)
      tstar  = kh(1)*dtdz/ustar
      qstar  = kq(1)*dqdz/ustar
      zgs    = z(1)
      if (ustar.gt.smax*vel1) ustar=smax*vel1
      if (abs(tstar).gt.smax*abs(      if (abs(qstar).gt.smax*abs(      if (ustar.lt.smin*vel1) ustar=smin*vel1
      if (abs(tstar).lt.smin*abs(      if (abs(qstar).lt.smin*abs(
      lmonin = ustar*ustar*tgrnd/(kappa*grav*tstar)
      if(abs(lmonin).lt.teeny) lmonin=sign(teeny,lmonin)
c     To compute the drag coefficient,Stanton number and Dalton number
      call dflux(lmonin,ustar,vel1,z0m,z0h,z0q,zgs,cm,ch,cq,itype)

      return
      end subroutine stars


      subroutine getl1(e,zhat,dzh,lscale,n) 2
!@sum getl1 estimates the master length scale of the turbulence model
!@+   on the secondary grid
!@auth  Ye Cheng/G. Hartke
      implicit none
  
      integer, intent(in) :: n     !@var n  array dimension
      real*8, dimension(n-1), intent(in) :: e,zhat,dzh
      real*8, dimension(n-1), intent(out) :: lscale
      real*8, parameter :: alpha=0.2d0
  
      real*8 :: sum1,sum2,l0,l1
      integer :: j
  
      sum1=0.
      sum2=0.
      do j=1,n-1
        sum1=sum1+sqrt(        sum2=sum2+sqrt(      end do
      l0=alpha*sum1/sum2
  
      do j=1,n-1
        l1=kappa*zhat(j)
        lscale(j)=l0*l1/(l0+l1)
      end do
  
      return
      end subroutine getl1


      subroutine getl(e,u,v,t,zhat,dzh,lmonin,ustar,lscale,dbl,n) 2,1
!@sum   getl computes the master length scale of the turbulence model
!@+     on the secondary grid. l0 in this routine is 0.16*(pbl height)
!@+     according to the LES data (Moeng and Sullivan 1992)
!@auth  Ye Cheng/G. Hartke
!@ver   1.0
!@var e z-profle of turbulent kinetic energy
!@var u z-profle of west-east   velocity component
!@var v z-profle of south-north velocity component
!@var t z-profle of virtual potential temperature
!@var lscale z-profile of the turbulent dissipation length scale
!@var z vertical grids (main, meter)
!@var zhat vertical grids (secondary, meter)
!@var dzh(j)  z(j+1)-z(j)
!@var dbl PBL height (m)
!@var n number of vertical subgrid main layers

      USE CONSTANT, only : by3
      implicit none

      integer, intent(in) :: n   !@var n  array dimension
      real*8, dimension(n-1), intent(in) :: e,zhat,dzh
      real*8, dimension(n), intent(in) :: u,v,t
      real*8, dimension(n-1), intent(out) :: lscale
      real*8, intent(in) :: dbl,lmonin,ustar

      integer :: i   !@var i  array dimension
      real*8 kz,l0,ls,lb,an2,an,dudz,dvdz,as2,qty,qturb,zeta

      l0=.16d0*dbl ! Moeng and Sullivan 1994

      if (l0.lt.zhat(1)) l0=zhat(1)

      kz=kappa*zhat(1)
      lscale(1)=l0*kz/(l0+kz)

      do i=2,n-1
        kz=kappa*zhat(i)
        zeta=zhat(i)/lmonin
        ! Nakanishi (2001)
        if(zeta.ge.1.) then
          ls=kz/3.7d0
        elseif(zeta.ge.0.) then
          ls=kz/(1.+2.7d0*zeta)
        else
          ls=kz*(1.-100.*zeta)**0.2d0
        endif
        if (          an2=2.*grav*(          an=sqrt(an2)
          qturb=sqrt(2*e(i))
          if(zeta.ge.0.) then
             lb=qturb/an
          else
             qty=(ustar/((-kappa*lmonin*l0*l0)**by3*an))**0.5d0
             lb=qturb*(1.+5.*qty)/an
          endif
        else
          lb=1.d30
        endif
        lscale(i)=l0*ls*lb/(l0*ls+l0*lb+ls*lb)
      end do

      return
      end subroutine getl


      subroutine dflux(lmonin,ustar,vsurf,z0m,z0h,z0q,zgs, 1
     2                 cm,ch,cq,itype)
!@sum   dflux computes (dimensionless) surface fluxes of momemtun,
!@+     heat and moisture (drag coefficient, Stanton number,
!@+     and Dalton number)
!@auth  Ye Cheng/G. Hartke
!@ver   1.0
!@var lmonin = Monin-Obukhov length (m)
!@var ustar  = friction speed (sqrt of surface momentum flux) (m/sec)
!@var vsurf  = total surface wind speed, used to limit ustar -> cm
!@var zgs    = height of the surface layer (m)
!@var itype  = integer identifying surface type
!@var z0m   = momentum roughness length, prescribed (itype=3,4) (m)
!@var z0m   = roughness length for momentum, computed (itype=1,2)
!@var cm    = drag coefficient for momentum
!@var ch    = Stanton number
!@var cq    = Dalton number
!@var z0h   = roughness length for temperature (m)
!@var z0q   = roughness length for water vapor (m)
      implicit none

      real*8,  intent(in) :: lmonin,ustar,vsurf,zgs
      integer,  intent(in) :: itype
      real*8,  intent(inout) :: z0m
      real*8,  intent(out) :: cm,ch,cq,z0h,z0q

      real*8, parameter :: nu=1.5d-5,num=0.135d0*nu,nuh=0.395d0*nu,
     *     nuq=0.624d0*nu

      real*8 r0q,beta,zgsbyl,z0mbyl,z0hbyl,z0qbyl,cmn,chn,cqn,dpsim
     *     ,dpsih,dpsiq,xms,xm0,xhs,xh0,xqs,xq0,dm,dh,dq,lzgsbyz0m
     *     ,lzgsbyz0h,lzgsbyz0q

      if ((itype.eq.1).or.(itype.eq.2)) then
c *********************************************************************
c  Compute roughness lengths using rough surface formulation:
        z0m=num/ustar+0.018d0*ustar*ustar*bygrav
        if (z0m.gt.0.2d0) z0m=0.2d0
c       if (ustar.le.0.02d0) then
          z0h=nuh/ustar + 1.4d-5
          if (z0h.gt.0.5852d0) z0h=0.5852d0
          z0q=nuq/ustar + 1.3d-4
          if (z0q.gt.0.92444d0) z0q=0.92444d0
c         else
c         r0q=(ustar*z0m/nu)**0.25d0
c         z0h=7.4*z0m*exp(-2.4604d0*r0q)
c         z0q=7.4*z0m*exp(-2.2524d0*r0q)
c         if (ustar.lt.0.2d0) then
c           beta=(ustar-0.02d0)/0.18d0
c           z0h=(1.-beta)*nuh/ustar+beta*z0h
c           z0q=(1.-beta)*nuq/ustar+beta*z0q
c         endif
c       endif
c *********************************************************************
      else
c *********************************************************************
c  For land and land ice, z0m is specified. For z0h and z0q,
c    empirical evidence suggests:
        z0h=z0m*.13533528d0    ! = exp(-2.)
        z0q=z0h
c *********************************************************************
      endif

      lzgsbyz0m = log(zgs/z0m)
      lzgsbyz0h = log(zgs/z0h)
      lzgsbyz0q = log(zgs/z0q)
      cmn=kappa*kappa/(lzgsbyz0m**2)
      chn=kappa*kappa/(lzgsbyz0m*lzgsbyz0h)
      cqn=kappa*kappa/(lzgsbyz0m*lzgsbyz0q)

      zgsbyl=zgs/lmonin
      z0mbyl=z0m/lmonin
      z0hbyl=z0h/lmonin
      z0qbyl=z0q/lmonin

c *********************************************************************
c  Now compute DPSI, which is the difference in the psi functions
c    computed at zgs and the relevant roughness height:
c *********************************************************************

      if (lmonin.gt.0.) then
c *********************************************************************
c  Here the atmosphere is stable with respect to the ground:
        dpsim=-gamams*(zgsbyl-z0mbyl)
        dpsih= sigma1*lzgsbyz0h-sigma*gamahs*(zgsbyl-z0hbyl)
        dpsiq= sigma1*lzgsbyz0q-sigma*gamahs*(zgsbyl-z0qbyl)
c *********************************************************************
        else
c *********************************************************************
c  Here the atmosphere is unstable with respect to the ground:
        xms  =    (1.-gamamu*zgsbyl)**0.25d0
        xm0  =    (1.-gamamu*z0mbyl)**0.25d0
        xhs  =sqrt(1.-gamahu*zgsbyl)
        xh0  =sqrt(1.-gamahu*z0hbyl)
        xqs  =sqrt(1.-gamahu*zgsbyl)
        xq0  =sqrt(1.-gamahu*z0qbyl)
        dpsim=log((1.+xms)*(1.+xms)*(1.+xms*xms)/
     2           ((1.+xm0)*(1.+xm0)*(1.+xm0*xm0)))-
     3        2.*(atan(xms)-atan(xm0))
        dpsih=sigma1*lzgsbyz0h+2.*sigma*log((1.+xhs)/(1.+xh0))
        dpsiq=sigma1*lzgsbyz0q+2.*sigma*log((1.+xqs)/(1.+xq0))
c *********************************************************************
      endif

      dm=      1./(1.-min(dpsim/lzgsbyz0m,.9d0))**2
      dh=sqrt(dm)/(1.-min(dpsih/lzgsbyz0h,.9d0))
      dq=sqrt(dm)/(1.-min(dpsiq/lzgsbyz0q,.9d0))

      cm=XCDpbl*dm*cmn
      ch=dh*chn
      cq=dq*cqn
      if (cm.gt.cmax) cm=cmax
      if (ch.gt.cmax) ch=cmax
      if (cq.gt.cmax) cq=cmax

      if (cm.lt.cmin) cm=cmin
      if (ch.lt.cmin) ch=cmin
      if (cq.lt.cmin) cq=cmin

      return
      end subroutine dflux


      subroutine simil(u,t,q,z,ustar,tstar,qstar, 1
     2                 z0m,z0h,z0q,lmonin,tg,qg)
!@sum   simil calculates the similarity solutions for wind speed,
!@+     virtual potential temperature, and moisture mixing ratio
!@+     at height z.
!@auth  Ye Cheng/G. Hartke
!@ver   1.0
!@var     z       height above ground at which solution is computed (m)
!@var     ustar   friction speed (m/sec)
!@var     tstar   temperature scale (K)
!@var     qstar   moisture scale
!@var     z0m     momentum roughness height (m)
!@var     z0h     temperature roughness height (m)
!@var     z0q     moisture roughness height (m)
!@var     lmonin  Monin-Obukhov length scale (m)
!@var     tg      ground temperature (K)
!@var     qg      ground moisture mixing ratio
!@var     u       computed similarity solution for wind speed (m/sec)
!@var     t       computed similarity solution for virtual potential
!@+               temperature (K)
!@var     q       computed similarity solution for moisture mixing ratio
      implicit none

      real*8,  intent(in) :: z,ustar,tstar,qstar,z0m,z0h,z0q
      real*8,  intent(in) :: lmonin,tg,qg
      real*8,  intent(out) :: u,t,q

      real*8 zbyl,z0mbyl,z0hbyl,z0qbyl,dpsim,dpsih,dpsiq,xm,xm0,xh,xh0
     *     ,xq,xq0,lzbyz0m,lzbyz0h,lzbyz0q

      zbyl  =z  /lmonin
      z0mbyl=z0m/lmonin
      z0hbyl=z0h/lmonin
      z0qbyl=z0q/lmonin
      lzbyz0m=log(z/z0m)
      lzbyz0h=log(z/z0h)
      lzbyz0q=log(z/z0q)
c *********************************************************************
c  Now compute DPSI, which is the difference in the psi functions
c    computed at zgs and the relevant roughness height:
c *********************************************************************

      if (lmonin.gt.0.) then
c *********************************************************************
c  Here the atmosphere is stable with respect to the ground:
        dpsim=-gamams*(zbyl-z0mbyl)
        dpsih= sigma1*lzbyz0h-sigma*gamahs*(zbyl-z0hbyl)
        dpsiq= sigma1*lzbyz0q-sigma*gamahs*(zbyl-z0qbyl)
c *********************************************************************
        else
c *********************************************************************
c  Here the atmosphere is unstable with respect to the ground:
        xm   =    (1.-gamamu*  zbyl)**0.25d0
        xm0  =    (1.-gamamu*z0mbyl)**0.25d0
        xh   =sqrt(1.-gamahu*  zbyl)
        xh0  =sqrt(1.-gamahu*z0hbyl)
        xq   =sqrt(1.-gamahu*  zbyl)
        xq0  =sqrt(1.-gamahu*z0qbyl)
        dpsim=log((1.+xm )*(1.+xm )*(1.+xm *xm )/
     2           ((1.+xm0)*(1.+xm0)*(1.+xm0*xm0)))-
     3        2.*(atan(xm )-atan(xm0))
        dpsih=sigma1*lzbyz0h+2.*sigma*log((1.+xh)/(1.+xh0))
        dpsiq=sigma1*lzbyz0q+2.*sigma*log((1.+xq)/(1.+xq0))
c *********************************************************************
      endif

      u=   (ustar/kappa)*(lzbyz0m-dpsim)
      t=tg+(tstar/kappa)*(lzbyz0h-dpsih)
      q=qg+(qstar/kappa)*(lzbyz0q-dpsiq)

      return
      end subroutine simil


      subroutine griddr(z,zhat,xi,xihat,dz,dzh,z1,zn,bgrid,n,ierr) 2,3
!@sum Computes altitudes on the vertical grid. The XI coordinates are
!@+   uniformly spaced and are mapped in a log-linear fashion onto the
!@+   Z grid. (The Z's are the physical coords.) Also computes the
!@+   altitudes on the secondary grid, ZHAT(I), and the derivatives
!@+   dxi/dz evaluated at both all Z(I) and ZHAT(I).
!@auth  Ye Cheng/G. Hartke
!@ver   1.0
c     Grids:
c
c                n   - - - - - - - - - - - - -
c                    -------------------------  n-1
c                n-1 - - - - - - - - - - - - -
c                    -------------------------  j+1
c     (main,z)   j+1 - - - - - - - - - - - - -
c     (z)            -------------------------  j     (secondary, zhat)
c                j   - - - - - - - - - - - - -
c                    -------------------------  j-1
c                j-1 - - - - - - - - - - - - -
c                    -------------------------    2
c                2   - - - - - - - - - - - - -
c                    -------------------------    1
c                1   - - - - - - - - - - - - -
c
c     dz(j)==zhat(j)-zhat(j-1), dzh(j)==z(j+1)-z(j)
!@var bgrid determines how strongly non-linear the
!@+   mapping is. BGRID=0 gives linear mapping. Increasing BGRID
!@+   packs more points into the bottom of the layer.
!@+   Bgrid is calculated in the beginning of this subroutine every
!@+   time this subroutine is called. zs is the first grid height
!@+   and dzs is the grid separation near the surface. The values of
!@+   parameters byzs(=1/zs) and bydzs(=1/dzs) are obtained by
!@+   calling an old version of this subroutine with bgrid=0.2927
!@+   and has been tested as appropriate for pe(2)=934mb.
!@+   Now we impose that for other pe(2)s, zs and dzs be the same
!@+   as when pe(2)=934mb, so as to maintain the balance between
!@+   the accuracy and stability.
!@+   The new value of bgrid is then calculated below which
!@+   also depends on ztop (i.e., depends on pe(2)).
!@var  z       height of main grids (meter)
!@var  zhat    height of secondary grids (meter)
!@var  xi      an uniformly spaced coordinate mapped to z
!@var  xihat   an uniformly spaced coordinate mapped to zhat
!@var  dz   dxi/(dxi/dz)
!@var  dzh  dxi/(dxi/dzh)
!@var  dxi  (ztop - zbottom)/(n-1)
!@var  ierr Error reporting flag
      implicit none

      integer, intent(in) :: n    !@var n  array dimension
      real*8, dimension(n), intent(out) :: z,xi,dz
      real*8, dimension(n-1), intent(out) :: zhat,xihat,dzh
      integer, intent(out) :: ierr
      real*8, intent(in) :: z1,zn
      real*8, intent(out) :: bgrid

      real*8, parameter ::  tolz=1d-3
     &  ,byzs=1.d0/10.d0,bydzs=1.d0/4.7914d0
      real*8 z1pass,znpass,b,xipass,lznbyz1
      common /grids_99/z1pass,znpass,b,xipass,lznbyz1
!$OMP  THREADPRIVATE(/GRIDS_99/)
      real*8, external :: fgrid2
      real*8 rtsafe
      integer i,j,iter  !@var i,j,iter loop variable
      real*8 dxi,zmin,zmax,dxidz,dxidzh

      z1pass=z1
      znpass=zn
      dxi=(zn-z1)/float(n-1)
      bgrid=max((dxi*bydzs-1.)/((zn-z1)*byzs-log(zn/z1)),0.d0)

      b=bgrid
      zmin=z1
      zmax=zn

      do i=1,n-1
        xi(i)=z1+(zn-z1)*float(i-1)/float(n-1)
        xihat(i)=z1+(zn-z1)*(float(i)-0.5)/float(n-1)
      end do
      xi(n)=zn
      z(1)=z1

      lznbyz1 = log(zn/z1)
      dxidz=1.+bgrid*((zn-z1)/z1-lznbyz1)
      dz(1)=dxi/dxidz
      xipass=xihat(1)
      zhat(1)=rtsafe(fgrid2,zmin,zmax,tolz,ierr)
      if (ierr.gt.0) return
      dxidzh=1.+bgrid*((zn-z1)/zhat(1)-lznbyz1)
      dzh(1)=dxi/dxidzh

      do i=2,n-1
        xipass=xi(i)
        z(i)=rtsafe(fgrid2,zmin,zmax,tolz,ierr)
        if (ierr.gt.0) return
        xipass=xihat(i)
        zhat(i)=rtsafe(fgrid2,zmin,zmax,tolz,ierr)
        if (ierr.gt.0) return
        dxidz=1.+bgrid*((zn-z1)/z(i)-lznbyz1)
        dxidzh=1.+bgrid*((zn-z1)/zhat(i)-lznbyz1)
        dz(i)=dxi/dxidz
        dzh(i)=dxi/dxidzh
      end do
      z(n)=zn
      dxidz=1.+bgrid*((zn-z1)/zn-lznbyz1)
      dz(n)=dxi/dxidz

      return
      end subroutine griddr


      subroutine tfix(t,z,ttop,tgrnd,lmonin,tstar,ustar,khs,ts_guess,n) 1
!@sum   tfix
!@auth  Ye Cheng/G. Hartke
!@ver   1.0

      implicit none

      integer, intent(in) :: n    !@var n  array dimension
      real*8, dimension(n),intent(in) :: z
      real*8, dimension(n),intent(inout) :: t
      real*8, intent(in) :: ttop,tgrnd,ustar,khs,ts_guess
      real*8, intent(inout) :: lmonin,tstar

      real*8 dtdz
      integer i  !@var i loop variable

      t(1)=ts_guess
      do i=2,n-1
        t(i)=t(1)+(      end do

      dtdz   = (      tstar  = khs*dtdz/ustar
      if (abs(tstar).gt.smax*abs(      if (abs(tstar).lt.smin*abs(
      lmonin = ustar*ustar*tgrnd/(kappa*grav*tstar)
      if(abs(lmonin).lt.teeny) lmonin=sign(teeny,lmonin)
      
      return
      end subroutine tfix


      subroutine ccoeff0 1
!@sum   ccoeff0 sets/calculates model coefficients for the
!@+     Giss 2000 turbulence model (level2 2/2.5)
!@auth  Ye Cheng
!@ver   1.0
      implicit none

      ! temperary variable
      real*8 :: del

      prt=     0.82d0
      b1=     19.3d0
      b2=     15.8d0
      b123=b1**(2./3.)
      g1=       .1070d0
      g2=       .0032d0
      g3=       .0864d0
      g4=       .1000d0
      g5=     11.04d0
      g6=       .786d0
      g7=       .643d0
      g8=       .547d0
c
      d1=(7.*g4/3+g8)/g5
      d2=(g3**2-g2**2/3.)-1./(4.*g5**2)*(g6**2-g7**2)
      d3=g4/(3.*g5**2)*(4.*g4+3.*g8)
      d4=g4/(3.*g5**2)*(g2*g6-3.*g3*g7-g5*(g2**2-g3**2))
     &   +g8/g5*(g3**2-g2**2/3.)
      d5=-1./(4.*g5**2)*(g3**2-g2**2/3)*(g6**2-g7**2)
      s0=g1/2.
      s1=-g4/(3.*g5**2)*(g6+g7)+2.*g4/(3.*g5)*(g1-g2/3.-g3)
     &   +g1/(2.*g5)*g8
      s2=-g1/(8.*g5**2)*(g6**2-g7**2)
      s4=2./(3.*g5)
      s5=2.*g4/(3.*g5**2)
      s6=2./(3.*g5)*(g3**2-g2**2/3)-g1/(2.*g5)*(g3-g2/3.)
     &   +g1/(4*g5**2)*(g6-g7)

c     find rimax:

      c1=s5+2*d3
      c2=s1-s6-2*d4
      c3=-s2+2.*d5
      c4=s4+2.*d1
      c5=-s0+2.*d2

      rimax=(c2+sqrt(c2**2-4.*c1*c3))/(2.*c1)
      rimax=int(rimax*1000.)/1000.

c     find ghmin,ghmax,gmmax0:

      del=(s4+2.*d1)**2-8.*(s5+2.*d3)
      ghmin=(-s4-2.*d1+sqrt(del))/(2.*(s5+2.*d3))
      ghmin=int(ghmin*10000.)/10000.
      ghmax=(b1*0.53d0)**2
      gmmax0=(b1*0.34d0)**2  ! not in use yet

      ! for level 3 model only:
      g0=2./3.
      d1_3=(7.*g4/3.)/g5
      d2_3=d2
      d3_3=g4/(3.*g5**2)*(4.*g4)
      d4_3=g4/(3.*g5**2)*(g2*g6-3*g3*g7-g5*(g2**2-g3**2))
      d5_3=d5
      s0_3=s0
      s1_3=-g4/(3.*g5**2)*(g6+g7)+2.*g4/(3.*g5)*(g1-g2/3.-g3)
      s2_3=s2
      s3_3=g0*g4/g5*(g3+g2/3.+1./(2.*g5)*(g6+g7))
      s4_3=s4
      s5_3=s5
      s6_3=s6

      return
      end subroutine ccoeff0


      subroutine getk(km,kh,kq,ke,gma,gha,u,v,t,e,lscale,dzh,n) 2,1
!@sum   getk calculates eddy diffusivities Km, Kh and Ke
!@+     Giss 2000 turbulence model at level 2.5
!@auth  Ye Cheng
!@ver   1.0
c     Grids:
c
c                n   - - - - - - - - - - - - -
c                    -------------------------  n-1
c                n-1 - - - - - - - - - - - - -
c                    -------------------------  j+1
c     (main,z)   j+1 - - - - - - - - - - - - -
c     (z)            -------------------------  j     (secondary, zhat)
c                j   - - - - - - - - - - - - -
c                    -------------------------  j-1
c                j-1 - - - - - - - - - - - - -
c                    -------------------------    2
c                2   - - - - - - - - - - - - -
c                    -------------------------    1
c                1   - - - - - - - - - - - - -
c
c     dz(j)==zhat(j)-zhat(j-1), dzh(j)==z(j+1)-z(j)
c     at main: u,v,t,q,ke
c     at edge: e,lscale,km,kh,gm,gh
!@sum getk computes the turbulent viscosity, Km, and turbulent
!@+   conductivity, Kh, and turbulent diffusivity , Ke,
!@+   using the GISS second order closure model (2000)
!@+   at main: u,v,t,q,ke
!@+   at secondary: e,lscale,km,kh,gma,gha
!@auth  Ye Cheng/G. Hartke
!@ver   1.0
!@var u,v,t,e,lscale,t_real z-profiles
!@var dz(j) zhat(j)-zhat(j-1)
!@var dzh(j)  z(j+1)-z(j)
!@var km turbulent viscosity for u and v equations
!@var kh turbulent conductivity for t and q equations
!@var ke turbulent diffusivity for e equation
!@var gma normalized velocity gradient, tau**2*as2
!@var gha normalized temperature gradient, tau**2*an2
!@var tau B1*lscale/sqrt(2*e)
!@var as2 shear squared, (dudz)**2+(dvdz)**2
!@var an2 Brunt-Vaisala frequency, grav/T*dTdz
!@var se stability constant for e, adjustable

      USE CONSTANT, only : teeny

      implicit none

      integer, intent(in) :: n    !@var n  array dimension
      real*8, dimension(n), intent(in) :: u,v,t
      real*8, dimension(n-1), intent(in) :: e,lscale,dzh
      real*8, dimension(n-1), intent(out) :: km,kh,kq,ke,gma,gha

      real*8, parameter :: se=0.1d0,kmax=100.d0
     &  ,kmmin=1.5d-5,khmin=2.5d-5,kqmin=2.5d-5,kemin=1.5d-5
      real*8 :: an2,dudz,dvdz,as2,ell,den,qturb,tau,gh,gm,gmmax,sm,sh
     &  ,sq,sq_by_sh,taue
      integer :: i,j  !@var i,j loop variable

      do i=1,n-1
        an2=2.*grav*(        dudz=(        dvdz=(        as2=dudz*dudz+dvdz*dvdz
        ell=lscale(i)
        qturb=sqrt(2.*e(i))
        tau=B1*ell/max(qturb,teeny)
        gh=tau*tau*an2
        gm=tau*tau*as2
        if(gh.lt.ghmin) gh=ghmin
        if(gh.gt.ghmax) gh=ghmax
        gmmax=(1+d1*gh+d3*gh*gh)/(d2+d4*gh)
        if(gm.gt.gmmax) gm=gmmax
        den=1.+d1*gh+d2*gm+d3*gh*gh+d4*gh*gm+d5*gm*gm
        sm=(s0+s1*gh+s2*gm)/den
        sh=(s4+s5*gh+s6*gm)/den
        sq=sh
        taue=tau*e(i)
        km(i)=min(max(taue*sm,kmmin),kmax)
        kh(i)=min(max(taue*sh,khmin),kmax)
        kq(i)=min(max(taue*sq,kqmin),kmax)
        ke(i)=min(max(taue*se,kemin),kmax)
        gma(i)=gm
        gha(i)=gh
      end do
      return
      end subroutine getk


      subroutine e_eqn(esave,e,u,v,t,km,kh,ke,lscale, 1,1
     &                     dz,dzh,ustar,dtime,n)
!@sum e_eqn integrates differential eqn for e (tridiagonal method)
!@+   between the surface and the first GCM layer.
!@+   The boundary conditions at the bottom are:
!@+   e(1)=(1/2)*B1**(2/3)*ustar**2
!@+   at the top, dedz is continuous
!@auth Ye Cheng/G. Hartke
!@ver  1.0
!@var u z-profle of west-east   velocity component
!@var v z-profle of south-north velocity component
!@var t z-profle of virtual potential temperature
!@var km z-profile of turbulent viscosity
!@var kh z-profile of turbulent conductivity
!@var ke z-profile of turbulent diffusion in eqn for e
!@var lscale z-profile of the turbulent dissipation length scale
!@var z vertical grids (main, meter)
!@var dz(j) zhat(j)-zhat(j-1)
!@var dzh(j)  z(j+1)-z(j)
!@var dtime time step
!@var ustar friction velocity at the surface
!@var n number of vertical subgrid main layers
      implicit none

      integer, intent(in) :: n    !@var n  array dimension
      real*8, dimension(n) :: sub,dia,sup,rhs

      real*8, intent(in) :: ustar, dtime
      real*8, dimension(n), intent(in) :: u,v,t,dz
      real*8, dimension(n-1), intent(in) :: esave,km,kh,ke,lscale,dzh
      real*8, dimension(n-1), intent(inout) :: e

      real*8 :: an2,dudz,dvdz,as2,qturb,ri,aa,bb,cc,gm,tmp
      integer :: j !@var j loop variable
c
c     sub(j)*e_jm1_kp1+dia(j)*e_j_kp1+sup(j)*e_jp1_kp1 = rhs(j)
c     from kirk:/u/acyxc/papers/2ndOrder/maple/phik.1,
c       sub(j)=-dtime*aj*P1a_jm1_k/dxi^2
c       sup(j)=-dtime*aj*P1a_jp1_k/dxi^2
c       dia(j)=1-(sub(j)+dia(j))+dtime*P3_j_k
c       rhs(j)=PHI_j_k + dtime*P4_j_k
c       where P1a == P1*a
c       aj=dxi/dz(j) if j refers to the primary grid
c       aj=dxi/dzh(j) if j refers to the secondary grid
c       now for e_eqn j refers to the secondary grid, therefore:
c       aj = dxi/dzh(j)
c       a(j-1/2) = dxi/dzh(j-1/2)=dxi/dz(j)
c       a(j+1/2) = dxi/dzh(j+1/2)=dxi/dz(j+1)
c       p1(j-1/2)=0.5*(ke(j)+ke(j-1))
c       p1(j+1/2)=0.5*(ke(j)+ke(j+1))
c       ke(j)=sq*lscale(j)*qturb, qturb=sqrt(2.*e(j))
c
      do j=2,n-2
          qturb=sqrt(2.*e(j))
          sub(j)=-dtime*0.5d0*(ke(j)+ke(j-1))/(dzh(j)*dz(j))
          sup(j)=-dtime*0.5d0*(ke(j)+ke(j+1))/(dzh(j)*dz(j+1))
          dia(j)=1.-(sub(j)+sup(j))+dtime*2*qturb/(b1*lscale(j))
          an2=2*grav*(          dudz=(          dvdz=(          as2=dudz*dudz+dvdz*dvdz
          rhs(j)=esave(j)+dtime*(km(j)*as2-kh(j)*an2)
       end do

      dia(1)=1.
      sup(1)=0.
      rhs(1)=0.5d0*b123*ustar*ustar

      j=n-1
      an2=2.*grav*(      dudz=(      dvdz=(      as2=max(dudz*dudz+dvdz*dvdz,teeny)
      ri=an2/as2
      if(ri.gt.rimax) ri=rimax
      aa=c1*ri*ri-c2*ri+c3
      bb=c4*ri+c5
      cc=2.d0
      if(abs(aa).lt.1d-8) then
        gm= -cc/bb
      else
        tmp=bb*bb-4.*aa*cc
        gm=(-bb-sqrt(tmp))/(2.*aa)
      endif
      sub(n-1)=0.
      dia(n-1)=1.
      rhs(n-1)=max(0.5d0*(B1*lscale(j))**2*as2/max(gm,teeny),teeny)

c     sub(n-1)=-1.
c     dia(n-1)=1.
c     rhs(n-1)=0.

      call TRIDIAG(sub,dia,sup,rhs,e,n-1)

      do j=1,n-1
         e(j)=min(max(      end do

      Return
      end subroutine e_eqn


      subroutine t_eqn(u,v,t0,t,q,z,kh,kq,dz,dzh,ch,usurf,tgrnd 1,1
     &                ,ttop,dtime,n)
!@sum t_eqn integrates differential eqn for t (tridiagonal method)
!@+   between the surface and the first GCM layer.
!@+   The boundary conditions at the bottom are:
!@+   kh * dt/dz = ch * usurf * (t - tg) - deltx * t * wq
!@+   at the top, the virtual potential temperature is prescribed.
!@auth Ye Cheng/G. Hartke
!@ver  1.0
!@var u z-profle of west-east   velocity component
!@var v z-profle of south-north velocity component
!@var t z-profle of virtual potential temperature
!@var q z-profle of specific humidity
!@var t0 z-profle of t at previous time step
!@var kh z-profile of heat conductivity
!@var kq z-profile of moisture diffusivity
!@var z vertical grids (main, meter)
!@var dz(j) zhat(j)-zhat(j-1)
!@var dzh(j)  z(j+1)-z(j)
!@var ch  dimensionless heat flux at surface (stanton number)
!@var usurf effective surface velocity
!@var tgrnd virtual potential temperature at the ground
!@var ttop virtual potential temperature at the first GCM layer
!@var dtime time step
!@var n number of vertical subgrid main layers
      implicit none

      integer, intent(in) :: n
      real*8, dimension(n) :: sub,dia,sup,rhs

      real*8, dimension(n), intent(in) :: u,v,t0,q,z,dz
      real*8, dimension(n-1), intent(in) :: dzh,kh,kq
      real*8, dimension(n), intent(inout) :: t
      real*8, intent(in) :: ch,tgrnd
      real*8, intent(in) :: ttop,dtime,usurf

      real*8 :: facth,factx,facty
      integer :: i,j,iter  !@var i,j,iter loop variable

      do i=2,n-1
         sub(i)=-dtime/(dz(i)*dzh(i-1))*kh(i-1)
         sup(i)=-dtime/(dz(i)*dzh(i))*kh(i)
         dia(i)=1.-(sub(i)+sup(i))
      end do

      factx=(dpdxr-dpdxr0)/(      facty=(dpdyr-dpdyr0)/(      do i=2,n-1
        rhs(i)=t0(i)-dtime*t(i)*bygrav*(      end do

      facth  = ch*usurf*dzh(1)/kh(1)

      dia(1) = 1.+facth
     &        +deltx*kq(1)/kh(1)*(      sup(1) = -1.
      rhs(1) = facth*tgrnd

      dia(n)  = 1.
      sub(n)  = 0.
      rhs(n) = ttop

      call TRIDIAG(sub,dia,sup,rhs,t,n)

      return
      end subroutine t_eqn


      subroutine q_eqn(q0,q,kq,dz,dzh,cq,usurf,qgrnd,qtop,dtime,n 1,2
     &     ,flux_max,fr_sat)
!@sum q_eqn integrates differential eqn q (tridiagonal method)
!@+   between the surface and the first GCM layer.
!@+   The boundary conditions at the bottom are:
!@+   kq * dq/dz = min ( cq * usurf * (q - qg) ,
!@+         fr_sat * cq * usurf * (q - qg) + ( 1 - fr_sat ) * flux_max )
!@+   at the top, the moisture is prescribed.
!@auth Ye Cheng/G. Hartke
!@ver  1.0
!@var q z-profle of specific humidity
!@var q0 z-profle of q at previous time step
!@var kq z-profile of moisture diffusivity
!@var dz(j) zhat(j)-zhat(j-1)
!@var dzh(j)  z(j+1)-z(j)
!@var cq  dimensionless moisture flux at surface (dalton number)
!@var usurf effective surface velocity
!@var qgrnd specific humidity at the ground
!@var qtop specific humidity at the first GCM layer
!@var dtime time step
!@var n number of vertical subgrid main layers
!@var flux_max maximal flux from the unsaturated soil
!@var fr_sat fraction of the saturated soil
      implicit none

      integer, intent(in) :: n
      real*8, dimension(n) :: sub,dia,sup,rhs

      real*8, dimension(n), intent(in) :: q0,dz
      real*8, dimension(n-1), intent(in) :: dzh,kq
      real*8, dimension(n), intent(out) :: q
      real*8, intent(in) :: cq,qgrnd,qtop,dtime,usurf
      real*8, intent(in) :: flux_max,fr_sat

      real*8 :: factq
      integer :: i  !@var i loop variable

      do i=2,n-1
         sub(i)=-dtime/(dz(i)*dzh(i-1))*kq(i-1)
         sup(i)=-dtime/(dz(i)*dzh(i))*kq(i)
         dia(i)=1.-(sub(i)+sup(i))
      end do

      do i=2,n-1
        rhs(i)=q0(i)
      end do

      factq  = cq*usurf*dzh(1)/kq(1)

      dia(1) = 1.+factq
      sup(1) = -1.
      rhs(1)= factq*qgrnd

      dia(n)  = 1.
      sub(n)  = 0.
      rhs(n) = qtop

      call TRIDIAG(sub,dia,sup,rhs,q,n)

c**** Now let us check if the computed flux doesn't exceed the maximum
c**** for unsaturated fraction

      if ( fr_sat .ge. 1. ) return   ! all soil is saturated
      if ( cq * usurf * (qgrnd - q(1)) .le. flux_max ) return

c**** Flux is too high, have to recompute with the following boundary
c**** conditions at the bottom:
c**** kq * dq/dz = fr_sat * cq * usurf * (q - qg)
c****              + ( 1 - fr_sat ) * flux_max

      dia(1) = 1. + fr_sat*factq
      sup(1) = -1.
      rhs(1)= fr_sat*factq*qgrnd + (1.-fr_sat)*flux_max*dzh(1)/kq(1)

      call TRIDIAG(sub,dia,sup,rhs,q,n)

      return
      end subroutine q_eqn


      subroutine tr_eqn(tr0,tr,kq,dz,dzh,sfac,constflx,trtop,,1

     *     dtime,n)
!@sum tr_eqn integrates differential eqn for tracers (tridiag. method)
!@+   between the surface and the first GCM layer.
!@+   The boundary conditions at the bottom are:
!@+   kq * dtr/dz = sfac * trs - constflx
!@+   i.e. for moisture, sfac=cq*usurf, constflx=cq*usurf*qg
!@+        to get:  kq * dq/dz = cq * usurf * (qs - qg)
!@+   This should be flexible enough to deal with most situations.
!@+   at the top, the tracer conc. is prescribed.
!@auth Ye Cheng/G. Hartke
!@ver  1.0
!@var tr z-profle of tracer concentration
!@var tr0 z-profle of tr at previous time step
!@var kq z-profile of moisture diffusivity
!@var dz(j) zhat(j)-zhat(j-1)
!@var dzh(j)  z(j+1)-z(j)
!@var sfac factor multiplying surface tracer conc. in b.c.
!@var constflx the constant component of the surface tracer flux
!@var trtop tracer concentration at the first GCM layer
!@var dtime time step
!@var n number of vertical subgrid main layers

      implicit none

      integer, intent(in) :: n
      real*8, dimension(n) :: sub,dia,sup,rhs

      real*8, dimension(n), intent(in) :: tr0,dz
      real*8, dimension(n-1), intent(in) :: dzh,kq
      real*8, dimension(n), intent(out) :: tr
      real*8, intent(in) :: sfac,constflx,trtop,dtime

      real*8 :: facttr
      integer :: i  !@var i loop variable

      do i=2,n-1
        sub(i)=-dtime/(dz(i)*dzh(i-1))*kq(i-1)
        sup(i)=-dtime/(dz(i)*dzh(i))*kq(i)
        dia(i)=1.-(sub(i)+sup(i))
      end do

      do i=2,n-1
        rhs(i)=tr0(i)
      end do

      facttr  = dzh(1)/kq(1)

      dia(1) = 1.+sfac*facttr
      sup(1) = -1.
      rhs(1)= facttr*constflx

      dia(n)  = 1.
      sub(n)  = 0.
      rhs(n) = trtop

      call TRIDIAG(sub,dia,sup,rhs,tr,n)



      return
      end subroutine tr_eqn


      subroutine uv_eqn(u0,v0,u,v,z,km,dz,dzh, 1,2
     2                  ustar,cm,z0m,utop,vtop,dtime,coriol,
     3                  ug,vg,uocean,vocean,n)
!@sum uv_eqn integrates differential eqns for u & v (tridiagonal method)
!@+   between the surface and the first GCM layer.
!@+   The boundary conditions at the bottom are:
!@+   km * du/dz = cm * usurf * u
!@+   km * dv/dz = cm * usurf * v
!@+    at the top, the winds are prescribed.
!@auth Ye Cheng/G. Hartke
!@ver  1.0
!@var u z-profle of west-east   velocity component
!@var v z-profle of south-north velocity component
!@var u0 z-profle of u at previous time step
!@var v0 z-profle of v at previous time step
!@var km z-profile of turbulent viscosity
!@var kh z-profile of turbulent conductivity
!@var dz(j) zhat(j)-zhat(j-1)
!@var dzh(j)  z(j+1)-z(j)
!@var ustar friction velocity at the surface
!@var cm  dimensionless  momentum flux at surface (drag coeff.)
!@var ch  dimensionless heat flux at surface (stanton number)
!@var cq  dimensionless moisture flux at surface (dalton number)
!@var z0m  roughness height for momentum (if itype=1 or 2)
!@var utop due east component of the wind at the first GCM layer
!@var vtop due north component of the wind at the first GCM layer
!@var dtime time step
!@var coriol the Coriolis parameter
!@var ug due east component of the geostrophic wind
!@var vg due north component of the geostrophic wind
!@var n number of vertical subgrid main layers
      implicit none

      integer, intent(in) :: n
      real*8, dimension(n) :: sub,dia,sup,rhs,rhs1

      real*8, dimension(n), intent(in) :: u0,v0,z,dz
      real*8, dimension(n), intent(inout) :: u,v
      real*8, dimension(n-1), intent(in) :: km,dzh
      real*8, intent(in) :: ustar,cm,z0m,utop,vtop,dtime,coriol,ug,vg
     &                     ,uocean,vocean

      real*8 :: factx,facty,dpdx,dpdy,usurf,factor
      integer :: i,j,iter  !@var i,j,iter loop variable

      do i=2,n-1
         sub(i)=-dtime/(dz(i)*dzh(i-1))*km(i-1)
         sup(i)=-dtime/(dz(i)*dzh(i))*km(i)
         dia(i)=1.-(sub(i)+sup(i))
      end do

      factx=(dpdxr-dpdxr0)/(      facty=(dpdyr-dpdyr0)/(      do i=2,n-1
        dpdx=factx*(        dpdy=facty*(        rhs(i)=u0(i)+dtime*(coriol*v(i)-dpdx)
        rhs1(i)=v0(i)-dtime*(coriol*u(i)+dpdy)
c       rhs(i)=u0(i)+dtime*coriol*(v(i)-vg)
c       rhs1(i)=v0(i)-dtime*coriol*(u(i)-ug)
      end do

      usurf  = sqrt((u(1)-uocean)**2+(      factor = cm*usurf*dzh(1)/km(1)

      dia(1) = 1.+factor
      sup(1) = -1.
      rhs(1) = factor*uocean
      rhs1(1) = factor*vocean

      dia(n) = 1.
      sub(n) = 0.
      rhs(n)  = utop

      rhs1(1)  = 0.
      rhs1(n)  = vtop

      call TRIDIAG(sub,dia,sup,rhs,u,n)
      call TRIDIAG(sub,dia,sup,rhs1,v,n)

      return
      end subroutine uv_eqn


      subroutine t_eqn_sta(t,q,kh,kq,dz,dzh,ch,usurf,tgrnd,ttop,n) 1,1
!@sum  t_eqn_sta computes the static solutions of t
!@+    between the surface and the first GCM layer.
!@+    The boundary conditions at the bottom are:
!@+    kh * dt/dz = ch * usurf * (t - tg)
!@+    at the top, virtual potential temperature is prescribed.
!@auth Ye Cheng/G. Hartke
!@ver  1.0
!@var u z-profle of west-east   velocity component
!@var v z-profle of south-north velocity component
!@var t z-profle of virtual potential temperature
!@var q z-profle of specific humidity
!@var kh z-profile of heat conductivity
!@var kq z-profile of specific humidity
!@var z vertical grids (main, meter)
!@var dz(j) zhat(j)-zhat(j-1)
!@var dzh(j)  z(j+1)-z(j)
!@var ch  dimensionless heat flux at surface (stanton number)
!@var usurf effective surface velocity
!@var tgrnd virtual potential temperature at the ground
!@var ttop virtual potential temperature at the first GCM layer
!@var n number of vertical main subgrid layers
      implicit none

      integer, intent(in) :: n
      real*8, dimension(n) :: sub,dia,sup,rhs

      real*8, dimension(n), intent(in) :: q,dz
      real*8, dimension(n), intent(inout) :: t
      real*8, dimension(n-1), intent(in) :: kh,kq,dzh
      real*8, intent(in) :: ch,tgrnd,ttop,usurf

      real*8 :: facth
      integer :: i !@var i loop variable

      do i=2,n-1
         sub(i)=-1./(dz(i)*dzh(i-1))*kh(i-1)
         sup(i)=-1./(dz(i)*dzh(i))*kh(i)
         dia(i)=-(sub(i)+sup(i))
      end do

      do i=2,n-1
        rhs(i)=0.
      end do

      facth  = ch*usurf*dzh(1)/kh(1)

      dia(1) = 1.+facth
     &        +deltx*kq(1)/kh(1)*(      sup(1) = -1.
      rhs(1) = facth*tgrnd

      dia(n)  = 1.
      sub(n)  = 0.
      rhs(n) = ttop

      call TRIDIAG(sub,dia,sup,rhs,t,n)

      return
      end subroutine t_eqn_sta


      subroutine q_eqn_sta(q,kq,dz,dzh,cq,usurf,qgrnd,qtop,n) 1,1
!@sum  q_eqn_sta computes the static solutions of q
!@+    between the surface and the first GCM layer.
!@+    The boundary conditions at the bottom are:
!@+    kq * dq/dz = cq * usurf * (q - qg)
!@+    at the top, moisture is prescribed.
!@auth Ye Cheng/G. Hartke
!@ver  1.0
!@var u z-profle of west-east   velocity component
!@var v z-profle of south-north velocity component
!@var t z-profle of virtual potential temperature
!@var q z-profle of specific humidity
!@var kq z-profile of moisture diffusivity
!@var z vertical grids (main, meter)
!@var dz(j) zhat(j)-zhat(j-1)
!@var dzh(j)  z(j+1)-z(j)
!@var cq  dimensionless moisture flux at surface (dalton number)
!@var usurf effective surface velocity
!@var qgrnd specific humidity at the ground
!@var qtop specific humidity at the first GCM layer
!@var n number of vertical main subgrid layers
      implicit none

      integer, intent(in) :: n
      real*8, dimension(n) :: sub,dia,sup,rhs

      real*8, dimension(n), intent(in) :: dz
      real*8, dimension(n), intent(inout) :: q
      real*8, dimension(n-1), intent(in) :: kq,dzh
      real*8, intent(in) :: cq,qgrnd,qtop,usurf

      real*8 :: factq
      integer :: i  !@var i loop variable

      do i=2,n-1
         sub(i)=-1./(dz(i)*dzh(i-1))*kq(i-1)
         sup(i)=-1./(dz(i)*dzh(i))*kq(i)
         dia(i)=-(sub(i)+sup(i))
      end do

      do i=2,n-1
        rhs(i)=0.
      end do

      factq  = cq*usurf*dzh(1)/kq(1)

      dia(1) = 1.+factq
      sup(1) = -1.
      rhs(1)= factq*qgrnd

      dia(n)  = 1.
      sub(n)  = 0.
      rhs(n) = qtop

      call TRIDIAG(sub,dia,sup,rhs,q,n)

      return
      end subroutine q_eqn_sta


      subroutine uv_eqn_sta(u,v,z,km,dz,dzh, 1,2
     2            ustar,cm,utop,vtop,coriol,uocean,vocean,n)
!@sum  uv_eqn_sta computes the static solutions of the u and v
!@+    between the surface and the first GCM layer.
!@+    The boundary conditions at the bottom are:
!@+    km * du/dz = cm * usurf * u
!@+    km * dv/dz = cm * usurf * v
!@+    at the top, the winds are prescribed.
!@auth Ye Cheng/G. Hartke
!@ver  1.0
!@var u z-profle of west-east   velocity component
!@var v z-profle of south-north velocity component
!@var z vertical height at the main grids (meter)
!@var zhat vertical height at the secondary grids (meter)
!@var km z-profile of turbulent viscosity
!@var kh z-profile of turbulent conductivity
!@var dz(j) zhat(j)-zhat(j-1)
!@var dzh(j)  z(j+1)-z(j)
!@var ustar friction velocity at the surface
!@var cm  dimensionless  momentum flux at surface (drag coeff.)
!@var ch  dimensionless heat flux at surface (stanton number)
!@var cq  dimensionless moisture flux at surface (dalton number)
!@var utop due east component of the wind at the first GCM layer
!@var vtop due north component of the wind at the first GCM layer
!@var coriol the Coriolis parameter
!@var n number of vertical subgrid main layers
      implicit none

      integer, intent(in) :: n
      real*8, dimension(n) :: sub,dia,sup,rhs,rhs1

      real*8, dimension(n), intent(in) :: z,dz
      real*8, dimension(n), intent(inout) :: u,v
      real*8, dimension(n-1), intent(in) :: km,dzh
      real*8, intent(in) :: ustar,cm,utop,vtop,coriol,uocean,vocean

      real*8 :: factx,facty,dpdx,dpdy,usurf,factor
      integer :: i,j,iter  !@var i,j,iter loop variable

      do i=2,n-1
          sub(i)=-1./(dz(i)*dzh(i-1))*km(i-1)
          sup(i)=-1./(dz(i)*dzh(i))*km(i)
          dia(i)=-(sub(i)+sup(i))
       end do

      factx=(dpdxr-dpdxr0)/(      facty=(dpdyr-dpdyr0)/(c      write(99,*) factx,facty
      do i=2,n-1
        dpdx=factx*(        dpdy=facty*(        rhs(i)=(coriol*v(i)-dpdx)
        rhs1(i)=-(coriol*u(i)+dpdy)
c       rhs(i)=coriol*(v(i)-vg)
c       rhs1(i)=-coriol*(u(i)-ug)
      end do

      usurf  = sqrt((u(1)-uocean)**2+(      factor = cm*usurf*dzh(1)/km(1)

      dia(1) = 1.+factor
      sup(1) = -1.
      rhs(1) = factor*uocean
      rhs1(1) = factor*vocean

      dia(n) = 1.
      sub(n) = 0.
      rhs(n)  = utop

      rhs1(1)  = 0.
      rhs1(n)  = vtop

      call TRIDIAG(sub,dia,sup,rhs,u,n)
      call TRIDIAG(sub,dia,sup,rhs1,v,n)

      return
      end subroutine uv_eqn_sta


      subroutine level2(e,u,v,t,lscale,dzh,n) 1,1
!@sum  level2 computes the turbulent kinetic energy e using the
!@+    GISS 2000 turbulence model at level 2
!@auth  Ye Cheng/G. Hartke
!@ver   1.0
!@var e z-profle of turbulent kinetic energy
!@var u z-profle of west-east   velocity component
!@var v z-profle of south-north velocity component
!@var t z-profle of virtual potential temperature
!@var lscale z-profile of the turbulent dissipation length scale
!@var z vertical grids (main, meter)
!@var dzh(j)  z(j+1)-z(j)
!@var n number of vertical subgrid main layers

      USE CONSTANT, only : teeny

      implicit none

      integer, intent(in) :: n     !@var n array dimension
      real*8, dimension(n), intent(in)   :: u,v,t
      real*8, dimension(n-1), intent(in) :: lscale,dzh
      real*8, dimension(n-1), intent(out) :: e
      real*8 :: dudz,dvdz,as2,an2,ri,gm
      real*8 :: aa,bb,cc,tmp
      integer :: i    !@var i loop variable

      do i=1,n-1
        an2=2.*grav*(        dudz=(        dvdz=(        as2=max(dudz*dudz+dvdz*dvdz,teeny)
        ri=an2/as2
        if(ri.gt.rimax) ri=rimax
        aa=c1*ri*ri-c2*ri+c3
        bb=c4*ri+c5
        cc=2.d0
        if(abs(aa).lt.1d-8) then
          gm= -cc/bb
        else
          tmp=bb*bb-4.*aa*cc
          gm=(-bb-sqrt(tmp))/(2.*aa)
        endif
        e(i)=0.5d0*(B1*lscale(i))**2*as2/max(gm,teeny)
        e(i)=min(max(      end do

      return
      end subroutine level2


      subroutine inits(tgrnd,qgrnd,zgrnd,zgs,ztop,utop,vtop, 1,14
     2          ttop,qtop,coriol,cm,ch,cq,ustar,uocean,vocean,
     3          ilong,jlat,itype)
!@sum  inits initializes the winds, virtual potential temperature,
!@+    and humidity using static solutions of the GISS 2000
!@+    turbulence model at level 2
!@var  n number of sub-grid levels for the PBL
!@var  tgrnd virtual potential temperature of ground,at roughness height
!@var  qgrnd  moisture at the ground, at the roughness height
!@var  zgrnd
!@var  zgs  height of the surface layer (nominally 10 m)
!@var  ztop height of the first model layer, approx 200 m if lm=9
!@var  utop  x component of wind at the top of the layer
!@var  vtop  y component of wind at the top of the layer
!@var  ttop  virtual potential temperature at the top of the layer
!@var  qtop  moisture at the top of the layer
!@var  coriol the Coriolis parameter
!@var  cm  dimensionless  momentum flux at surface (drag coeff.)
!@var  ch  dimensionless heat flux at surface (stanton number)
!@var  cq  dimensionless moisture flux at surface (dalton number)
!@var  bgrid log-linear gridding parameter
!@var  ustar  friction speed
!@var  ilong  longitude identifier
!@var  jlat  latitude identifier
!@var  itype  surface type

!@var  iprint longitude for diagnostics
!@var  jprint latitude for diagnostics
      implicit none

      real*8, intent(in) :: tgrnd,qgrnd,zgrnd,zgs,ztop,utop,vtop,ttop
     *     ,qtop,coriol,uocean,vocean
      real*8, intent(out) :: cm,ch,cq,ustar
      integer, intent(in) :: ilong,jlat,itype

      real*8, dimension(n-1) :: km,kh,kq,ke,gm,gh
      real*8, dimension(n) :: z,dz,xi,usave,vsave,tsave,qsave
      real*8, dimension(n-1) :: zhat,xihat,dzh,lscale,esave
      real*8 :: lmonin,bgrid,z0m,z0h,z0q,hemi,psi1,psi0,psi
     *     ,usurf,tstar,qstar,ustar0,dtime,test
     *     ,wstar3fac,wstar3,wstar2h,usurfq,usurfh
      integer, save :: iter_count=0
      integer, parameter ::  itmax=100
      integer, parameter ::  iprint=0,jprint=41 ! set iprint>0 to debug
      real*8, parameter ::  w=0.50,tol=1d-3
      integer :: i,j,iter,ierr  !@var i,j,iter loop variable

c**** special threadprivate common block (compaq compiler stupidity)
      real*8, dimension(n) :: u,v,t,q
      real*8, dimension(n-1) :: e
      common/pbluvtq/u,v,t,q,e

!$OMP  THREADPRIVATE (/pbluvtq/)
C**** end special threadprivate common block

      dbl=1000.d0 !initial guess of dbl
      z0m=zgrnd
      z0h=z0m
      z0q=z0m

      call griddr(z,zhat,xi,xihat,dz,dzh,zgs,ztop,bgrid,n,ierr)
      if (ierr.gt.0) then
        print*,"In inits: i,j,itype =",ilong,jlat,itype,tgrnd,qgrnd
        call stop_model("PBL error in inits",255)
      end if

c Initialization for iteration:
      if (coriol.le.0.) then
        hemi=-1.
        else
        hemi= 1.
      endif
      if (tgrnd.gt.ttop) then
        psi1=hemi*5.*radian
        else
        psi1=hemi*15.*radian
      endif
      psi0=atan2(vtop,utop+teeny)
      if(psi0.lt.0.) psi0=psi0+2.*pi
      psi=psi0+psi1
      usurf=0.4d0*sqrt(utop*utop+vtop*vtop)
      u(1)=usurf*cos(psi)
      v(1)=usurf*sin(psi)
      t(1)=tgrnd-0.5d0*(tgrnd-ttop)
      q(1)=qgrnd-0.5d0*(qgrnd-qtop)
      e(1)=1.d-3
        usave(1)=u(1)
        vsave(1)=v(1)
        tsave(1)=t(1)
        qsave(1)=q(1)
        esave(1)=e(1)
      u(n)=utop
      v(n)=vtop
      t(n)=ttop
      q(n)=qtop
      e(n-1)=1.d-3

      do i=2,n-1
        u(i)=u(1)+(        v(i)=v(1)+(        t(i)=t(1)+(        q(i)=q(1)+(        e(i)=e(1)+(        usave(i)=u(i)
        vsave(i)=v(i)
        tsave(i)=t(i)
        qsave(i)=q(i)
        esave(i)=e(i)
      end do

      ustar0=0.
      do iter=1,itmax

        if(iter.eq.1) then
          call getl1(e,zhat,dzh,lscale,n)
        else
          call getl(e,u,v,t,zhat,dzh,lmonin,ustar,lscale,dbl,n)
        endif

        call getk(km,kh,kq,ke,gm,gh,u,v,t,e,lscale,dzh,n)
        call stars(ustar,tstar,qstar,lmonin,tgrnd,qgrnd,
     2             u,v,t,q,z,z0m,z0h,z0q,cm,ch,cq,
     3             km,kh,kq,dzh,itype,n)
        !! dbl=.375d0*sqrt(ustar*abs(lmonin)/omega)
        !@+ M.J.Miller et al. 1992, J. Climate, 5(5), 418-434, Eqs(6-7),
        !@+ for heat and mositure
        if(          wstar3fac=-dbl*grav*2.*(          wstar2h = (wstar3fac*kh(1))**twoby3
        else
          wstar2h = 0.
        endif
        usurfh  = sqrt((u(1)-uocean)**2+(        usurfq  = usurfh

        call t_eqn_sta(t,q,kh,kq,dz,dzh,ch,usurfh,tgrnd,ttop,n)

        call q_eqn_sta(q,kq,dz,dzh,cq,usurfq,qgrnd,qtop,n)

        call uv_eqn_sta(u,v,z,km,dz,dzh,ustar,cm,utop,vtop,coriol,
     &                  uocean,vocean,n)

        call tcheck(t,tgrnd,n)
        call tcheck(q,qgrnd,n)
        call ucheck(u,v,z,ustar,lmonin,z0m,hemi,psi0,psi1,n)

        test=abs(2.*(ustar-ustar0)/(ustar+ustar0))
        if (test.lt.tol) exit

        call level2(e,u,v,t,lscale,dzh,n)

        do i=1,n-1
          u(i)=w*usave(i)+(1.-w)*u(i)
          v(i)=w*vsave(i)+(1.-w)*v(i)
          t(i)=w*tsave(i)+(1.-w)*t(i)
          q(i)=w*qsave(i)+(1.-w)*q(i)
          e(i)=w*esave(i)+(1.-w)*e(i)
          usave(i)=u(i)
          vsave(i)=v(i)
          tsave(i)=t(i)
          qsave(i)=q(i)
          esave(i)=e(i)
        end do
        ustar0=ustar

      end do
c     iter_count=iter_count+min(iter,itmax)
c     write(96,*) "iter_count in inits =",iter_count, iter,dbl

c     call check1(ustar,1,ilong,jlat,1)

      if ((ilong.eq.iprint).and.(jlat.eq.jprint)) then
        call output(u,v,t,q,e,lscale,z,zhat,dzh,
     2              km,kh,kq,ke,gm,gh,cm,ch,cq,z0m,z0h,z0q,
     3              ustar,tstar,qstar,lmonin,tgrnd,qgrnd,
     4              utop,vtop,ttop,qtop,
     5              dtime,bgrid,ilong,jlat,iter,itype,n)
      endif

      return
      end subroutine inits


      subroutine tcheck(t,tgrnd,n) 2
!@sum   tcheck checks for reasonable temperatures
!@auth  Ye Cheng/G. Hartke
!@ver   1.0
c ----------------------------------------------------------------------
c This routine makes sure that the temperature remains within
c  reasonable bounds during the initialization process. (Sometimes the
c  the computed temperature iterated out in left field someplace,
c  *way* outside any reasonable range.) This routine keeps the temp
c  between the maximum and minimum of the boundary temperatures.
c ----------------------------------------------------------------------
      implicit none

      integer, intent(in) :: n    !@var n array dimension
      real*8, dimension(n),intent(inout) ::  t !@var t temperature array
      real*8, intent(in) :: tgrnd  !@var tgrnd ground temperature

      integer, dimension(20) :: imax,imin
      integer nmin,nmax
      real*8 tmin,tmax,dt
      integer i                 !@var i  loop and dummy variables

      if (tgrnd.lt.t(n)) then
        tmin=tgrnd
        tmax=t(n)
      else
        tmin=t(n)
        tmax=tgrnd
      endif
      nmin=0
      nmax=0

      do i=1,n-1
        if (          nmin=nmin+1
          imin(nmin)=i
        endif
        if (          nmax=nmax+1
          imax(nmax)=i
        endif
      end do

      dt=0.01*(tmax-tmin)
      if (nmin.gt.0) then
        do i=1,nmin
          t(imin(i))=tmin+float(i)*dt
       end do
      endif

      if (nmax.gt.0) then
        do i=1,nmax
          t(imax(i))=tmax-float(i)*dt
       end do
      endif

      return
      end subroutine tcheck


      subroutine ucheck(u,v,z,ustar,lmonin,z0m,hemi,psi0,psi1,n) 1
!@sum  ucheck makes sure that the winds remain within reasonable
!@+    bounds during the initialization process. (Sometimes the computed
!@+    wind speed iterated out in left field someplace, *way* outside
!@+    any reasonable range.) Tests and corrects both direction and
!@+    magnitude of the wind rotation with altitude. Tests the total
!@+    wind speed via comparison to similarity theory. Note that it
!@+    works from the top down so that it can assume that at level (i),
!@+    level (i+1) displays reasonable behavior.
!@auth  Ye Cheng/G. Hartke
!@ver   1.0
      implicit none
      real*8, parameter :: psistb=15.*radian, psiuns=5.*radian

      integer, intent(in) :: n  !@var n array dimension
      real*8, dimension(n),intent(in) :: z
      real*8, dimension(n),intent(inout) :: u,v

      real*8, intent(in) :: lmonin,z0m,hemi,psi0,psi1,ustar

      real*8 psilim,psirot,angle,utotal,x,x0,dpsim,psiu,zbyl,z0byl,utest
      integer i                 !@var i  loop and dummy variables

      if (lmonin.ge.0.) then
        psilim=psistb
        else
        psilim=psiuns
      endif

c First, check the rotation of the wind vector:

      do i=n-1,1,-1
        psiu=atan2(        if(psiu.lt.0.) psiu=psiu+2.*pi
        psirot=psiu-psi0
c --------------------------------------------------------------------
c Next, check the magnitude of the wind. If the gradient is incorrect,
c  set the wind magnitude to that given by similarity theory:

        utotal=sqrt(        utest =sqrt(        if (utotal.gt.utest) then
          zbyl=z(i)/lmonin
          z0byl=z0m/lmonin
          if (lmonin.gt.0.) then
            dpsim=-gamams*(zbyl-z0byl)
            else
            x  = (1.-gamamu* zbyl)**0.25d0
            x0 = (1.-gamamu*z0byl)**0.25d0
            dpsim=log((1.+x )*(1.+x )*(1.+x *x )/
     2               ((1.+x0)*(1.+x0)*(1.+x0*x0)))-
     3            2.*(atan(x)-atan(x0))
          endif
          utotal=(ustar/kappa)*(log(          if (utotal.gt.utest) utotal=0.95d0*utest
          u(i)=utotal*cos(psiu)
          v(i)=utotal*sin(psiu)
        endif

        if (hemi*psirot.lt.0.) then
          angle=psi0+psi1*float(n-i)/float(n-1)
          u(i)=utotal*cos(angle)
          v(i)=utotal*sin(angle)
          go to 100
        endif

        if (hemi*psirot.gt.psilim) then
          angle=psi0+hemi*psilim*float(n-i)/float(n-1)
          u(i)=utotal*cos(angle)
          v(i)=utotal*sin(angle)
        endif

 100  end do

      return
      end subroutine ucheck


      subroutine check1(a,n,ilong,jlat,id),1
!@sum   check1 checks for NaN'S and INF'S in real 1-D arrays.
!@auth  Ye Cheng/G. Hartke
!@ver   1.0
      implicit none

      integer, intent(in) :: n      !@var n  array dimension
      integer, intent(in) :: ilong  !@var ilong  longitude identifier
      integer, intent(in) :: jlat   !@var jlat  latitude identifier
      integer, intent(in) :: id     !@var n  integer id
      real*8, dimension(n), intent(in) ::  a !@var a  real array

      integer i,k !@var i,k loop and dummy variables
      character*16 str  !@var str  output string

      do i=1,n
        write(str,'(e16.8)')a(i)
        k=index(str,'N')+index(str,'n')
        if (k.ne.0) then
          write (99,1000) ilong,jlat,i,id,a(i)
          if (id.lt.100) call stop_model('check1',255)
        endif
      end do

      return
 1000 format (1x,'check1:  ilong = ',6x,i3,/,1x,
     2            '         jlat  = ',6x,i3,/,1x,
     3            '         i     = ',6x,i3,/,1x,
     4            '         id    = ',6x,i3,/,1x,
     5            '         value = ',1pe11.4,/)
      end subroutine check1


      subroutine output(u,v,t,q,e,lscale,z,zhat,dzh, 2,1
     2                  km,kh,kq,ke,gm,gh,cm,ch,cq,z0m,z0h,z0q,
     3                  ustar,tstar,qstar,lmonin,tgrnd,qgrnd,
     4                  utop,vtop,ttop,qtop,
     5                  dtime,bgrid,ilong,jlat,iter,itype,n)
!@sum   output produces output for diagnostic purposes
!@auth  Ye Cheng/G. Hartke
!@ver   1.0
!@calls simil
      implicit none
      real*8, parameter :: degree=1./radian

      integer, intent(in) :: n,itype,iter,jlat,ilong
      real*8,  intent(in) :: lscale(n-1),lmonin

      real*8, dimension(n), intent(in) :: u,v,t,q,z
      real*8, dimension(n-1), intent(in) :: e,zhat,dzh,km,kh,kq,ke,gm,gh

      real*8, intent(in) :: cm,ch,cq,z0m,z0h,z0q,ustar,tstar,qstar,tgrnd
     *     ,qgrnd,utop,vtop,ttop,qtop,dtime,bgrid

      integer i  !@var i loop variable

      real*8 :: psitop,psi,utotal,utot1,utot2,sign,shear,tgrad,qgrad
     *     ,egrad,uflux,hflux,qflux,bvfrq2,shear2,rich,dqdz,dtdz,dudz
     *     ,phim,phih,dudzs,dtdzs,dqdzs,uratio,tratio,qratio,dbydzh
     *     ,tgradl,prod,utest,ttest,qtest

      write (99,5000) ilong,jlat,itype,dtime,iter,ustar,tstar,qstar,
     2                lmonin,tgrnd,qgrnd,cm,ch,cq,z0m,z0h,z0q,
     3                utop,vtop,ttop,qtop,
     4                bgrid
      write (99,1000)
      psitop=atan2(      if (psitop.lt.0.) psitop=psitop+360.
      do i=1,n
        psi=atan2(        if (psi.lt.0.) psi=psi+360.
        psi=psi-psitop
        utotal=sqrt(        call simil(utest,ttest,qtest,z(i),ustar,tstar,qstar,
     2             z0m,z0h,z0q,lmonin,tgrnd,qgrnd)
        write (99,3000) i,z(i),u(i),v(i),psi,utotal,utest,
     2                    t(i),ttest,q(i),qtest
      end do
      write (99,9000)
      write (99,2000)
      do i=1,n-1
        utot1=u(i+1)*u(i+1)+v(i+1)*v(i+1)
        utot2=u(i)  *u(i)  +v(i)  *v(i)
        if (utot1.ge.utot2) then
          sign=1.
          else
          sign=-1.
        endif
        shear =sqrt((u(i+1)-u(i))*(     2             +(        tgrad =(        qgrad =(        egrad =(        uflux=km(i)*shear*sign
        hflux=kh(i)*tgrad
        qflux=kh(i)*qgrad
        bvfrq2=grav*log(        shear2=((u(i+1)-u(i))/dzh(i))**2+
     2         ((v(i+1)-v(i))/dzh(i))**2+1d-8
        rich=bvfrq2/shear2
        write (99,3000) i,zhat(i),e(i),lscale(i),km(i),kh(i),
     2                    gm(i),gh(i),uflux,hflux,qflux,rich
      end do
      write (99,9000)
      write (99,7000)
      do i=1,n-1
        utot1=sqrt(        utot2=sqrt(  u(i)*  u(i)+  v(i)*  v(i))
        dudz=(utot1-utot2)/dzh(i)
        dtdz=(        dqdz=(        if (lmonin.lt.0.) then
          phim = 1./((1.-gamamu*zhat(i)/lmonin)**0.25)
          phih = sigma/sqrt(1.-gamahu*zhat(i)/lmonin)
          else
          phim = 1.+gamams*zhat(i)/lmonin
          phih = sigma*(1.+gamahs*zhat(i)/lmonin)
        endif
        dudzs=ustar*phim/(kappa*zhat(i))
        dtdzs=tstar*phih/(kappa*zhat(i))
        dqdzs=qstar*phih/(kappa*zhat(i))
        uratio=dudzs/(dudz+teeny)
        tratio=dtdzs/(dtdz+teeny)
        qratio=dqdzs/(dqdz+teeny)

        dbydzh=1./dzh(i)
        shear2=dbydzh*dbydzh*((u(i+1)-u(i))**2+
     2                        (        tgradl=dbydzh*log(        prod=km(i)*shear2-grav*kh(i)*tgradl

        write (99,3000) i,zhat(i),dudz,dudzs,dtdz,dtdzs,dqdz,dqdzs,
     2                            uratio,tratio,qratio,prod
      end do
      write (99,9000)
      write (99,9000)
c ----------------------------------------------------------------------
      return
1000  format (1x,'   i','      z     ','      U     ','      V     ',
     2                   '     psi    ','    Utotal  ','    Utest   ',
     2                   '      T     ','    Ttest   ','      Q     ',
     2                   '    Qtest   ',/)
2000  format (1x,'   i','     zhat   ','      E     ','      L     ',
     2                   '      Km    ','      Kh    ',
     3                   '      Gm    ','      Gh    ','    uflux   ',
     4                   '    hflux   ','    qflux   ','    rich #  ',/)
3000  format (1x,i4,11(1x,1pe11.4))
5000  format (1x,'i      = ',9x,i2,/,1x,
     2            'j      = ',9x,i2,/,1x,
     2            'itype  = ',9x,i2,/,1x,
     2            'dtime  = ',1pe11.4,/,1x,
     3            'iter   = ',7x,i4,/,1x,
     4            'ustar  = ',1pe11.4,/,1x,
     5            'tstar  = ',1pe11.4,/,1x,
     5            'qstar  = ',1pe11.4,/,1x,
     5            'lmonin = ',1pe11.4,/,1x,
     5            'tgrnd  = ',1pe11.4,/,1x,
     5            'qgrnd  = ',1pe11.4,/,1x,
     6            'cm     = ',1pe11.4,/,1x,
     6            'ch     = ',1pe11.4,/,1x,
     6            'cq     = ',1pe11.4,/,1x,
     6            'z0m    = ',1pe11.4,/,1x,
     7            'z0h    = ',1pe11.4,/,1x,
     8            'z0q    = ',1pe11.4,/,1x,
     6            'utop   = ',1pe11.4,/,1x,
     6            'vtop   = ',1pe11.4,/,1x,
     7            'ttop   = ',1pe11.4,/,1x,
     8            'qtop   = ',1pe11.4,/,1x,
     8            'bgrid  = ',1pe11.4,/)
7000  format (1x,'   i','     zhat   ','     dudz   ','   dudz sim ',
     2                   '     dtdz   ','   dtdz sim ','     dqdz   ',
     3                   '   dqdz sim ','    uratio  ','    tratio  ',
     4                   '    qratio  ','  production',/)
9000  format (1x)
      end subroutine output

      END MODULE SOCPBL


      subroutine fgrid2(z,f,df)
!@sum  fgrid2 computes functional relationship of z and xi + derivative
!@+    fgrid2 will be called in function rtsafe(fgrid2,x1,x2,xacc)
!@auth Ye Cheng/G. Hartke
!@ver  1.0
      implicit none
      real*8, intent(in) :: z
      real*8, intent(out) :: f,df
      real*8 z1,zn,bgrid,xi,lznbyz1
      common /grids_99/z1,zn,bgrid,xi,lznbyz1
!$OMP  THREADPRIVATE(/GRIDS_99/)

      f=z+bgrid*((zn-z1)*log(z/z1)-(z-z1)*lznbyz1)-xi
      df=1.+bgrid*((zn-z1)/z-lznbyz1)

      return
      end subroutine fgrid2


      function rtsafe(funcd,x1,x2,xacc,ierr) 3
!@sum   rtsafe use Newton-Rapheson + safeguards to solve F(x)=0
!@auth  Numerical Recipes
!@ver   1.0
      integer,parameter :: maxit=100
      real*8, intent(in) :: x1,x2,xacc
      integer, intent(out) :: ierr
      real*8 rtsafe
      external funcd
      integer j
      real*8 df,dx,dxold,f,fh,fl,temp,xh,xl

      ierr = 0
      call funcd(x1,fl,df)
      call funcd(x2,fh,df)
      if(fl*fh.gt.0.) then
        ierr = 1
        print*, 'Error: root must be bracketed in rtsafe'
      end if
      if(fl.eq.0.)then
        rtsafe=x1
        return
      else if(fh.eq.0.)then
        rtsafe=x2
        return
      else if(fl.lt.0.)then
        xl=x1
        xh=x2
      else
        xh=x1
        xl=x2
      endif
      rtsafe=.5*(x1+x2)
      dxold=abs(x2-x1)
      dx=dxold
      call funcd(rtsafe,f,df)
      do j=1,MAXIT
        if(((rtsafe-xh)*df-f)*((rtsafe-xl)*df-f).ge.0..or. abs(2.*
     *       f).gt.abs(dxold*df) ) then
          dxold=dx
          dx=0.5*(xh-xl)
          rtsafe=xl+dx
          if(xl.eq.rtsafe)return
        else
          dxold=dx
          dx=f/df
          temp=rtsafe
          rtsafe=rtsafe-dx
          if(temp.eq.rtsafe)return
        endif
        if(abs(dx).lt.xacc) return
        call funcd(rtsafe,f,df)
        if(f.lt.0.) then
          xl=rtsafe
        else
          xh=rtsafe
        endif
      end do
      ierr = 1
      print*, 'Error: rtsafe exceeding maximum iterations'
      return
      END

C**** Functions and subroutines for sub-gridscale wind distribution calc


      real*8 function sig(tke,mdf,dbl,delt,ch,ws,tsv) 54,1
!@sum calculate sigma for sub-grid scale wind distribution      
!@auth Reha Cakmur/Gavin Schmidt
      use constant, only : by3,grav
      implicit none
!@var tke local turbulent kinetic energy at surface (J/kg)
!@var mdf downdraft mass flux from moist convection (m/s)
!@var dbl boundary layer height (m)
!@var delt difference between surface and ground T (ts-tg) (K)
!@var tsv virtual surface temperature (K)
!@var ch local heat exchange coefficient 
!@var ws grid box mean wind speed (m/s)
      real*8, intent(in):: tke,mdf,dbl,delt,tsv,ch,ws
      real*8 wtke,wm,wd

C**** TKE contribution  ~ sqrt( 2/3 * tke)
      wtke=sqrt(2d0*tke*by3)

C**** dry convection/turbulence contribution ~(Qsens*g*H/rho*cp*T)^(1/3)
      if (delt.lt.0d0) then
        wd=(-delt*ch*ws*grav*dbl/tsv)**by3
      else
        wd=0.d0
      endif
C**** moist convection contribution beta/frac_conv = 200.
      wm=200.d0*mdf
C**** sigma
      sig=wm+wd+wtke
      return
      end


      subroutine integrate_sgswind(sig,wt,wmin,wmax,ws,icase,wint),4
!@sum Integrate sgswind distribution for different cases
!@auth Reha Cakmur/Gavin Schmidt
      use constant, only : by3
      implicit none
!@var sig distribution parameter (m/s)
!@var wt threshold velocity (m/s)
!@var ws resolved wind speed (m/s)
!@var wmin,wmax min and max wind speed for integral
      real*8, intent(in):: sig,wt,ws,wmin,wmax
      real*8, intent(out) :: wint
      integer, intent(in) :: icase
      integer, parameter :: nstep = 100
      integer i
      real*8 x,sgsw,bysig2

C**** integrate distribution from wmin to wmax using Simpsons' rule
C**** depending on icase, integral is done on w, w^2 or w^3
C**** Integration maybe more efficient with a log transform (putting
C**** more points near wmin), but that's up to you...
C**** Use approximate value for small sig and unresolved delta function
      if ((wmax-wmin).lt.sig*dble(nstep)) then
        bysig2=1./(sig*sig)
        wint=0
        do i=1,nstep+1
          x=wmin+(wmax-wmin)*(i-1)/dble(nstep)
          if (i.eq.1.or.i.eq.nstep+1) then
            wint=wint+sgsw(x,ws,wt,bysig2,icase)*by3
          elseif (mod(i,2).eq.0) then
            wint=wint+4d0*sgsw(x,ws,wt,bysig2,icase)*by3
          else
            wint=wint+2d0*sgsw(x,ws,wt,bysig2,icase)*by3
          end if
        end do
        wint=wint*(wmax-wmin)/dble(nstep)
      else                      ! approximate delta function
        if (ws.ge.wmin.and.ws.le.wmax) then
          wint = ws**(icase-1)*(ws-wt)
        else
          wint = 0.
        end if
      end if

      return
      end


      real*8 function sgsw(x,ws,wt,bysig2,icase) 3,2
!@sum sgsw function to be integrated for sgs wind calc
!@auth Reha Cakmur/Gavin Schmidt
      use constant, only : pi
      implicit none
      real*8, intent(in) :: x,ws,wt,bysig2
      integer, intent(in) :: icase
      real*8 :: besf,exx,bessi0

      besf=x*ws*bysig2
      if (besf.lt.200d0 ) then
        exx=exp(-0.5d0*(x*x+ws*ws)*bysig2)
        sgsw=(x)**(icase)*(x-wt)*bysig2*BESSI0(besf)*exx
      else ! use bessel function expansion for large besf
        sgsw=(x)**(icase)*(x-wt)*bysig2*exp(-0.5d0*(x-ws)**2*bysig2)
     *       /sqrt(besf*2*pi)
      end if
      
      return
      end


      real*8 function bessi0(x) 1
!@sum Bessel's function I0
!@auth  Numerical Recipes
      implicit none
      real*8 ax,y
      real*8, parameter:: p1=1.0d0, p2=3.5156229d0, p3=3.0899424d0,
     *     p4=1.2067492d0, p5=0.2659732d0, p6=0.360768d-1,
     *     p7=0.45813d-2
      real*8, parameter:: q1=0.39894228d0,q2=0.1328592d-1,
     *     q3=0.225319d-2, q4=-0.157565d-2, q5=0.916281d-2,
     *     q6=-0.2057706d-1, q7=0.2635537d-1, q8=-0.1647633d-1,
     *     q9=0.392377d-2
      real*8, intent(in)::x

      if (abs(x).lt.3.75d0) then
        y=(x/3.75d0)**2.d0
        bessi0=(p1+y*(p2+y*(p3+y*(p4+y*(p5+y*(p6+y*p7))))))
      else
        ax=abs(x)
        y=3.75d0/ax
        bessi0=(exp(ax)/sqrt(ax))*(q1+y*(q2+y*(q3+y*(q4
     *       +y*(q5+y*(q6+y*(q7+y*(q8+y*q9))))))))
      end if
      end function bessi0