!
!##################################################################
!##################################################################
!######                                                      ######
!######                SUBROUTINE SFCPHYSICS                 ######
!######                                                      ######
!######                     Developed by                     ######
!######     Center for Analysis and Prediction of Storms     ######
!######                University of Oklahoma                ######
!######                                                      ######
!##################################################################
!##################################################################
!


SUBROUTINE sfcphysics(nx,ny,nz,nstyps,                                  & 1,66
           u,v,w,ptprt,pprt,qv,qr,                                      &
           rhostr,ptbar,pbar,qvbar,                                     &
           x,y,z, zp, j1,j2,j3,j3inv, prcrate,                          &
           soiltyp,stypfrct,vegtyp,lai,roufns,veg,                      &
           tsfc,tsoil,wetsfc,wetdp,wetcanp,snowdpth,qvsfc,              &
           radsw,rnflx,                                                 &
           cdm,cdh,cdq,usflx,vsflx,ptsflx,qvsflx,pbldpth)
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  Calculate the surface momentum, heat and surface moisture fluxes.
!  These fluxes will be passed into turbulent mixing subroutines.
!  The drag coefficients are returned in cdm, cdh and cdq.
!
!  The soil-vegetation model is also called to update its state
!  variables.
!
!-----------------------------------------------------------------------
!
!  AUTHOR: Ming Xue
!  4/07/1992.
!
!  MODIFICATION HISTORY:
!
!  5/06/92 (M. Xue)
!  Added full documentation.
!
!  6/4/92 (M. Xue and H. Jin)
!  Further facelift.
!
!  9/14/1992 (M. Xue)
!  Different drag coefficients defined for momentum, temperature and
!  moisture. Definition of ptsfc0 and qvsfc0 changed. qvbar added to
!  the argument list.
!
!  9/28/1992 (M. Xue)
!  Sign correction in all the flux formulations.
!
!  8/28/1993 (M. Xue and Hao jin)
!  Modified the flux formulations considering the terrain effect.
!
!  01/24/94 (Y. Liu and V. Wong)
!  Added the surface energy and moisture budget model to predict the
!  ground temperature and moistures.
!
!  9/04/1994 (K. Brewster, V. Wong and X. Song)
!  Facelift.
!
!  11/01/1994 (Y. Liu)
!  Subroutine name of soil model were changed from SFCEBM to SOILEBM
!  and the arguments passed were changed to 2-D arrays instead of 3-D
!  temporary arrays.
!
!  12/07/1994 (Y. Liu)
!  Cleaned up internal documentation and an empty DO loop of labeled
!  360 in subroutine SFCFLXSD.
!
!  01/28/1995 (V. Wong and X. Song)
!  Add the option of stability and roughness dependent surface layer
!  parameterization for both land and sea surfaces.
!
!  02/07/1995 (Y. Liu)
!  Fixed a bug in the calculation of surface precipitation,
!  tem1(i,j,4).
!
!  02/27/95 (V. Wong, Y. Liu and X. Song)
!  Delete subroutines of calculating C_u and C_pt at the neutral state
!  Make predicting procedure more consistent. Add a limiting parameter
!  "blimit" to prevent the drag coefficients from becoming too small
!  in the stable region.
!
!  03/02/1995 (Y. Liu)
!  Added the minimum surface total wind speed to guarantee the basic
!  heat and moisture fluxes.
!
!  02/07/96 (V.Wong and X.Song)
!  Changed the calculation of roughness zo over the sea according to
!  Clarke(1970).  The new formulation is insensitive to the value of
!  the depth. of the first model layer above the ground.
!  Added kvwtr to denote the Von Karman constant over the sea.
!  Set a lower limiter, zolimit, for zo, and an upper limiter, z1limit,
!  for depth of the surface layer z1.
!
!  12/08/98 (Donghai Wang and Vince Wang)
!  Added the snow cover.
!
!  2/25/1999 (M. Xue)
!  Reorganized this subroutine from original SFCFLX.
!  Calculations of drag coefficients and fluxes are grouped into
!  into new subroutines DRAGCOEF and SFCFLX, respectively.
!  Named local arrays are defined for many 2D variables for better
!  readability.
!
!  2000/02/29 (Gene Bassett, Yang Kun, Vince Wong)
!  Corrected typo: "c8=gammahl*stabp" from equation (15) in Byun (1990).
!
!  2001/12/05 (Yunheng Wang)
!  Merged M. Xue's version to the official release.
!
!  2001/12/10 (Ming Xue)
!  Streamlined the usage of work arrays.
!
!-----------------------------------------------------------------------
!
!  INPUT:
!
!    nx       Number of grid points in the x-direction (east/west)
!    ny       Number of grid points in the y-direction (north/south)
!    nz       Number of grid points in the vertical
!
!    u        x component of velocity at a given time level (m/s)
!    v        y component of velocity at a given time level (m/s)
!    w        z component of velocity at a given time level (m/s)
!
!    ptprt    Perturbation potential temperature at a given time
!             level (K)
!    pprt     Perturbation pressure (Pascal)
!    qv       Water vapor specific humidity at a given time level
!             (kg/kg)
!    qr       Rain water mixing ratio at a given time level (kg/kg)
!    pbldpth  Planetary boundary layer depth (m) at a given time level
!    rhostr   Base state density rhobar times j3 (kg/m**3)
!    ptbar    Base state potential temperature (K)
!    pbar     Base state pressure (Pascal)
!    qvbar    Base state water vapor specific humidity (kg/kg)
!
!    x        The x-coord. of the physical and computational grid.
!             Defined at u-point.
!    y        The y-coord. of the physical and computational grid.
!             Defined at v-point.
!    zp       The physical height coordinate defined at w-point
!
!    j1       Coordinate transformation Jacobian -d(zp)/d(x)
!    j2       Coordinate transformation Jacobian -d(zp)/d(y)
!    j3       Coordinate transformation Jacobian, d(zp)/d(z)
!    j3inv    Inverse of j3
!
!    soiltyp  Soil type
!    vegtyp   Vegetation type
!    lai      Leaf Area Index
!    roufns   Surface roughness
!    veg      Vegetation fraction
!
!    radsw    Solar radiation reaching the surface
!    rnflx    Net radiation flux absorbed by the surface
!
!  OUTPUT:
!
!    cdm      Drag coefficient for momentum defined as scalar point
!    cdh      Drag coefficient for heat
!    cdq      Drag coefficient for moisture
!    usflx    Surface flux of u-momentum
!    vsflx    Surface flux of v-momentum
!    ptsflx   Surface flux of heat (K*kg/(m**2*s))
!    qvsflx   Surface flux of moisture (K*kg/(m**2*s))
!    pbldpth  Planetary boundary layer depth (m)
!
!  INPUT & OUTPUT:
!
!    tsfc     Skin temperature at the ground or ocean surface (K)
!    tsoil    Deep soil temperature (K) (in 1 m deep layer)
!    wetsfc   Surface soil moisture
!    wetdp    Deep soil moisture
!    wetcanp  Canopy moisture
!    snowdpth Snow depth (m)
!    qvsfc    Effective specific humidity at sfc.
!
!  WORK ARRAY:
!
!    ptsfc    Potential temperature at the ground level (K)
!    psfc     Surface pressure (Pascal)
!    wdspsfc  Wind speed at the surface
!    rhosfc   Base-state air density at the surface
!    pt1      Potential temperature at the first model level AGL
!    tair     Temperature of surface air
!    qvair    Specific humidity of surface air
!
!    tem1     3-dimensional array (nz >= 4) to store:
!    tem2     3-dimensional array (nz >= 4) to store:
!    tem3     3-dimensional array (nz >= 4) to store:
!    tem4     3-dimensional array (nz >= 4) to store:
!    tem5     3-dimensional array (nz >= 4) to store:
!
!-----------------------------------------------------------------------
!

!
!-----------------------------------------------------------------------
!
!  Variable Declarations.
!
!-----------------------------------------------------------------------
!
  IMPLICIT NONE

  INTEGER :: nx,ny,nz       ! The number grid points in 3 directions

  INCLUDE 'timelvls.inc'

  REAL :: u      (nx,ny,nz) ! Total u-velocity (m/s)
  REAL :: v      (nx,ny,nz) ! Total v-velocity (m/s)
  REAL :: w      (nx,ny,nz) ! Total w-velocity (m/s)

  REAL :: ptprt  (nx,ny,nz) ! Perturbation potential temperature (K)
  REAL :: pprt   (nx,ny,nz) ! Perturbation pressure (Pascal)
  REAL :: qv     (nx,ny,nz) ! Water vapor specific humidity (kg/kg)
  REAL :: qr     (nx,ny,nz) ! Rain water mixing ratio (kg/kg)

  REAL :: rhostr (nx,ny,nz) ! Base state density rhobar times j3.
  REAL :: ptbar  (nx,ny,nz) ! Base state potential temperature (K)
  REAL :: pbar   (nx,ny,nz) ! Base state pressure (Pascal)
  REAL :: qvbar  (nx,ny,nz) ! Base state water vapor specific
                            ! humidity (kg/kg)

  REAL :: x      (nx)       ! The x-coord. of the physical and
                            ! computational grid.
                            ! Defined at u-point.
  REAL :: y      (ny)       ! The y-coord. of the physical and
                            ! computational grid.
                            ! Defined at v-point.
  REAL :: z      (nz)       ! The z-coord. of the physical and
                            ! computational grid. Defined at w-point.

  REAL :: zp     (nx,ny,nz) ! The height of the terrain.

  REAL :: j1     (nx,ny,nz) ! Coordinate transformation
                            ! Jacobian -d(zp)/d(x)
  REAL :: j2     (nx,ny,nz) ! Coordinate transformation
                            ! Jacobian -d(zp)/d(y)
  REAL :: j3     (nx,ny,nz) ! Coordinate transformation
                            ! Jacobian  d(zp)/d(z)
  REAL :: j3inv  (nx,ny,nz) ! Inverse of j3

  REAL :: prcrate(nx,ny)    ! precipitation rate (kg/(m**2*s))

  INTEGER :: nstyps                ! Number of soil types
  INTEGER :: soiltyp(nx,ny,nstyps) ! Soil type at each point
  REAL :: stypfrct(nx,ny,nstyps)! Fraction of soil types
  INTEGER :: vegtyp (nx,ny)        ! Vegetation type at each point
  REAL :: lai    (nx,ny)           ! Leaf Area Index
  REAL :: roufns (nx,ny)           ! Surface roughness
  REAL :: veg    (nx,ny)           ! Vegetation fraction

  REAL :: tsfc   (nx,ny,0:nstyps)  ! Skin temperature at the ground 
                                   ! or ocean surface (K)
                                   ! (in 1 cm top layer)
  REAL :: tsoil  (nx,ny,0:nstyps)  ! Deep soil temperature (K)
                                   ! (in 1 m deep layer)
  REAL :: wetsfc (nx,ny,0:nstyps)  ! Surface soil moisture
  REAL :: wetdp  (nx,ny,0:nstyps)  ! Deep soil moisture
  REAL :: wetcanp(nx,ny,0:nstyps)  ! Canopy water amount

  REAL :: qvsfc  (nx,ny,0:nstyps)  ! Effective S.H. at sfc.
                                   ! at the ground level (K)
  REAL :: snowdpth(nx,ny)          ! Snow depth

  REAL :: radsw  (nx,ny)    ! Solar radiation flux reaching the surface
  REAL :: rnflx  (nx,ny)    ! Net radiation flux absorbed by surface

  REAL :: cdm    (nx,ny)    ! Drag coefficient for momentum
                            ! defined as scalar point
  REAL :: cdh    (nx,ny)    ! Drag coefficient for heat
  REAL :: cdq    (nx,ny)    ! Drag coefficient for moisture

  REAL :: usflx  (nx,ny)    ! surface flux of u-momentum
                            ! (kg/(m*s**2))
  REAL :: vsflx  (nx,ny)    ! surface flux of v-momentum
                            ! (kg/(m*s**2))
  REAL :: ptsflx (nx,ny)    ! surface flux of heat (K*kg/(m**2*s))
  REAL :: qvsflx (nx,ny)    ! surface flux of moisture (kg/(m**2*s))

  REAL :: pbldpth(nx,ny,nt) ! Planetary boundary layer depth (m)

  REAL :: ptsfc  (nx,ny)    ! Potential temperature at the ground level (K)
  REAL :: psfc   (nx,ny)    ! Surface pressure (Pascal)
  REAL :: wdspsfc(nx,ny)    ! Wind speed at the surface
  REAL :: rhosfc (nx,ny)    ! Base-state air density at the surface
  REAL :: pt1    (nx,ny)    ! Potential temperature at 1st model level AGL
  REAL :: tair   (nx,ny)    ! Temperature of surface air
  REAL :: qvair  (nx,ny)    ! Specific humidity of surface air

!-----------------------------------------------------------------------
! Local work arrays:
!-----------------------------------------------------------------------

  REAL :: shflx  (nx,ny)  ! Sensible heat flux
  REAL :: lhflx  (nx,ny)  ! Latent heat flux
  REAL :: gflx   (nx,ny)  ! Diffusive heat flux from ground surface to deep soil
  REAL :: ct     (nx,ny)  ! Thermal capacity
  REAL :: evaprg (nx,ny)  ! Evaporation from groud surface
  REAL :: evaprtr(nx,ny)  ! Transpiration of the remaining part (1-delta) of leaves
  REAL :: evaprr (nx,ny)  ! Direct evaporation from the fraction delta
  REAL :: qvsat  (nx,ny)  ! Surface specific humidity at saturation
  REAL :: qvsata (nx,ny)  ! qvsat(tair) (kg/kg)
  REAL :: f34    (nx,ny)  ! f34 and surface resistance

  REAL :: temxy1 (nx,ny)  ! Temporary array
  REAL :: temxy2 (nx,ny)  ! Temporary array
  REAL :: temxy3 (nx,ny)  ! Temporary array
  REAL :: temxy4 (nx,ny)  ! Temporary array

  REAL :: temxpy (nx+ny)  ! Temporary array
!
!-----------------------------------------------------------------------
!
!  Misc. local variables:
!
!-----------------------------------------------------------------------
!
  INTEGER :: i, j,is, nsmth

  REAL :: tema,temb
  REAL :: gumove,gvmove
  REAL :: nsfcstsave
  INTEGER :: lfn,tmstrln
  CHARACTER (LEN=80) :: soiloutfl,temchar,timsnd

  INTEGER :: mptag1,mptag2,mptag3,mptag4

  REAL :: tem
!
!-----------------------------------------------------------------------
!
!  Include files:
!
!-----------------------------------------------------------------------
!
  INCLUDE 'globcst.inc'
  INCLUDE 'phycst.inc'
  INCLUDE 'sfcphycst.inc'

  INCLUDE 'mp.inc'            ! Message passing parameters.
  INCLUDE 'bndry.inc'
  INCLUDE 'grid.inc'          ! Grid & map parameters.
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!

  IF ( sfcphy == 0 ) THEN
!
!-----------------------------------------------------------------------
!
!  If sfcphy = 0, usflx, vsflx, ptsflx and qvsflx are set to zero and
!  will not be used in the mixing subroutines.
!  For safety, we set them to zero.
!
!-----------------------------------------------------------------------
!
    DO j=1,ny
      DO i=1,nx
        usflx (i,j) = 0.0
        vsflx (i,j) = 0.0
        ptsflx(i,j) = 0.0
        qvsflx(i,j) = 0.0
      END DO
    END DO

    RETURN

  END IF

!
!-----------------------------------------------------------------------
!
!  Prepare surface arrays to be used later.
!
!-----------------------------------------------------------------------
!
  IF ( grdtrns == 0 ) THEN
    gumove = 0.0
    gvmove = 0.0
  ELSE
    gumove = umove
    gvmove = vmove
  END IF

  CALL set_acct(sfcphy_acct)

  DO j=1,ny-1
    DO i=1,nx-1
      tema = 0.25*(u(i+1,j,1)+u(i,j,1)+u(i+1,j,2)+u(i,j,2))+gumove
      temb = 0.25*(v(i,j+1,1)+v(i,j,1)+v(i,j+1,2)+v(i,j,2))+gvmove
      wdspsfc(i,j)=MAX(vsfcmin,                                         &
                   SQRT(tema*tema+temb*temb+w(i,j,2)*w(i,j,2)))

      psfc  (i,j) = 0.5 * ( pprt(i,j,1) + pbar(i,j,1)                   &
                          + pprt(i,j,2) + pbar(i,j,2) )

      rhosfc(i,j) = 0.5 * ( rhostr(i,j,1)*j3inv(i,j,1)                  &
                          + rhostr(i,j,2)*j3inv(i,j,2) )

      pt1(i,j)=ptbar(i,j,2)+ptprt(i,j,2)

      ptsfc (i,j) = tsfc(i,j,0)*(p0/psfc(i,j))**rddcp
    END DO
  END DO

  IF ( sfcphy == 1 ) THEN
!
!-----------------------------------------------------------------------
!
!  Compute the surface fluxes by calling SFCFLX given drag coefficients
!  ground surface potential temperature and equivalent specific humidity.
!
!-----------------------------------------------------------------------
!

    DO j=1,ny-1
      DO i=1,nx-1
        IF ( soiltyp(i,j,1) == 13 .AND.landwtr == 1) THEN
          cdm(i,j) = cdmwtr
          cdh(i,j) = cdhwtr
          cdq(i,j) = cdqwtr
        ELSE
          cdm(i,j) = cdmlnd
          cdh(i,j) = cdhlnd
          cdq(i,j) = cdqlnd
        END IF
      END DO
    END DO

    CALL sfcflx(nx,ny,nz,u,v,w,ptprt,pprt,qv,ptbar,pbar,qvbar,          &
         wdspsfc,rhosfc,ptsfc,qvsfc(1,1,0), cdm,cdh,cdq,                &
         usflx,vsflx,ptsflx,qvsflx)

  ELSE IF ( sfcphy == 2 ) THEN
!
!-----------------------------------------------------------------------
!
!  Calculate stability-dependent drag coefficients by calling DRAGCOEF.
!  Then calculate the fluxes.
!
!-----------------------------------------------------------------------
!

    CALL dragcoef(nx,ny,nz,                                             &
         zp, ptsfc,qvsfc,soiltyp,roufns,wdspsfc,pt1,                    &
         cdm,cdh,cdq)

    CALL sfcflx(nx,ny,nz,u,v,w,ptprt,pprt,qv,ptbar,pbar,qvbar,          &
         wdspsfc,rhosfc,ptsfc,qvsfc(1,1,0), cdm,cdh,cdq,                &
         usflx,vsflx,ptsflx,qvsflx)

  ELSE IF ( sfcphy == 3 .OR. sfcphy == 4 ) THEN
!
!-----------------------------------------------------------------------
!
!  Before calling SFCEBM, we need to prepare some data for the
!  surface model first.
!
!-----------------------------------------------------------------------
!
    DO j=1,ny-1
      DO i=1,nx-1
        tair(i,j) = 0.5 * ( (ptbar(i,j,1)+ptprt(i,j,1))                 &
                          + (ptbar(i,j,2)+ptprt(i,j,2)) )               &
                        * (psfc(i,j)/p0)**rddcp   ! Temperature
        qvair(i,j) = 0.5 * ( qv(i,j,1) + qv(i,j,2) )
      END DO
    END DO

    DO i=1,nx
      DO j=1,ny
        tsfc   (i,j,0) = 0.0
        tsoil  (i,j,0) = 0.0
        wetsfc (i,j,0) = 0.0
        wetdp  (i,j,0) = 0.0
        wetcanp(i,j,0) = 0.0
        qvsfc  (i,j,0) = 0.0
      END DO
    END DO

    DO j=1,ny
      DO i=1,nx
        usflx (i,j) = 0.0
        vsflx (i,j) = 0.0
        ptsflx(i,j) = 0.0
        qvsflx(i,j) = 0.0
      END DO
    END DO

    IF ( soilinitopt == 1 .AND. nstep == 1 ) THEN
      nsfcstsave = nsfcst
      nsfcst = nint(soiltintv/dtsfc)
    END IF

    DO is=nstyp,1,-1
!
!-----------------------------------------------------------------------
!
!  Specified drag coefficients cdm, cdh and cdq.
!
!-----------------------------------------------------------------------
!

      DO j=1,ny-1
        DO i=1,nx-1
          ptsfc (i,j) = tsfc(i,j,is)*(p0/psfc(i,j))**rddcp
        END DO
      END DO

      IF ( sfcphy == 3 ) THEN

        DO j=1,ny-1
          DO i=1,nx-1
            IF ( soiltyp(i,j,is) == 13 .AND.landwtr == 1) THEN
              cdm(i,j) = cdmwtr
              cdh(i,j) = cdhwtr
              cdq(i,j) = cdqwtr
            ELSE
              cdm(i,j) = cdmlnd
              cdh(i,j) = cdhlnd
              cdq(i,j) = cdqlnd
            END IF
          END DO
        END DO

      END IF

      IF ( sfcphy == 4 ) THEN
!
!-----------------------------------------------------------------------
!
!  Calculate drag coefficients and surface fluxes for each subtypes
!  within a grid cell. The final fluxes are weighted averaged of all
!  types. Drag coefficients are type dependent and are therefore
!  recalculated for each type.
!
!-----------------------------------------------------------------------
!
        CALL dragcoef(nx,ny,nz,zp,                                      &
             ptsfc,qvsfc(1,1,is),soiltyp(1,1,is),roufns,wdspsfc,pt1,    &
             cdm,cdh,cdq)
      END IF
!
!-----------------------------------------------------------------------
!
!  usflx, vsflx, ptsflx and qvsflx are stored in temxy1, temxy2, temxy3 
!  and temxy4, respectively. 
!
!-----------------------------------------------------------------------

      CALL sfcflx(nx,ny,nz,u,v,w,ptprt,pprt,qv,ptbar,pbar,qvbar,        &
           wdspsfc,rhosfc,ptsfc,qvsfc(1,1,is), cdm,cdh,cdq,             &
           temxy1,temxy2,temxy3,temxy4)

      DO i=2,nx-1
        DO j=1,ny-1
          usflx(i,j)  = usflx(i,j)  + temxy1(i,j)*stypfrct(i,j,is)
        END DO
      END DO

      DO i=1,nx-1
        DO j=2,ny-1
          vsflx(i,j)  = vsflx(i,j)  + temxy2(i,j)*stypfrct(i,j,is)
        END DO
      END DO

      DO i=1,nx-1
        DO j=1,ny-1
          ptsflx(i,j) = ptsflx(i,j) + temxy3(i,j)*stypfrct(i,j,is)
          qvsflx(i,j) = qvsflx(i,j) + temxy4(i,j)*stypfrct(i,j,is)
        END DO
      END DO

!-----------------------------------------------------------------------
!  Integrate soil model (soilebm) to update tsfc, qvsfc, tsoil,
!  wetsfc, wetdp, and wetcanp. They will be used by the next step.
!-----------------------------------------------------------------------

      CALL set_acct(soil_acct)

      CALL soilebm(nx,ny,nz,                                            &
           soiltyp(1,1,is),vegtyp,lai,veg,                              &
           tsfc(1,1,is),tsoil(1,1,is),wetsfc(1,1,is),                   &
           wetdp(1,1,is),wetcanp(1,1,is),snowdpth,qvsfc(1,1,is),        &
           wdspsfc,psfc,rhosfc,prcrate,                                 &
           tair,qvair,cdh,cdq,radsw,rnflx,                              &
           shflx,lhflx,gflx,ct,evaprg,evaprtr,evaprr,qvsat,qvsata,f34)

!-----------------------------------------------------------------------
!  If we want to use the sensible and latent heat fluxes calculated 
!  inside the soil model (soilebm) in the atmospheric model, doing the
!  folloing do loop. Otherwise, the flux calculated in sfcflx will be 
!  used, which is based on qvsfc one time step earlier. Otherwise,
!  they should be the same.
!-----------------------------------------------------------------------
!
!     DO i=1,nx-1
!       DO j=1,ny-1
!         ptsflx(i,j) = ptsflx(i,j) + shflx(i,j)*stypfrct(i,j,is)
!         qvsflx(i,j) = qvsflx(i,j) + lhflx(i,j)*stypfrct(i,j,is)
!       END DO
!     END DO

      DO i=1,nx-1
        DO j=1,ny-1
          tsfc   (i,j,0) = tsfc   (i,j,0)                               &
                         + tsfc   (i,j,is)*stypfrct(i,j,is)
          tsoil  (i,j,0) = tsoil  (i,j,0)                               &
                         + tsoil  (i,j,is)*stypfrct(i,j,is)
          wetsfc (i,j,0) = wetsfc (i,j,0)                               &
                         + wetsfc (i,j,is)*stypfrct(i,j,is)
          wetdp  (i,j,0) = wetdp  (i,j,0)                               &
                         + wetdp  (i,j,is)*stypfrct(i,j,is)
          wetcanp(i,j,0) = wetcanp(i,j,0)                               &
                         + wetcanp(i,j,is)*stypfrct(i,j,is)
          qvsfc  (i,j,0) = qvsfc  (i,j,0)                               &
                         + qvsfc  (i,j,is)*stypfrct(i,j,is)
        END DO
      END DO

    END DO  ! loop over is

!-----------------------------------------------------------------------
!  If we want to use the sensible and latent heat fluxes calculated 
!  inside the soil model (soilebm) in the atmospheric model, doing the
!  folloing do loop. Otherwise, the flux calculated in sfcflx will be 
!  used, which is based on qvsfc one time step earlier. Otherwise,
!  they should be the same.
!
!  Convert sensible and latent heat fluxes stored currently in 
!  ptsflx and qvsflx into pt and qv fluxes required by atmospheric model. 
!-----------------------------------------------------------------------
!
!   tem = 1.0/lathv
!   DO i=1,nx-1
!     DO j=1,ny-1
!       ptsflx(i,j)=ptsflx(i,j)*(p0/(pbar(i,j,1)+ptbar(i,j,2)          &
!                   +pprt(i,j,1)+pprt(i,j,2))*0.5)**rddcp
!       qvsflx(i,j)=qvsflx(i,j)*tem
!     ENDDO
!   ENDDO

    IF ( soilinitopt == 1 .AND. nstep == 1 ) THEN
      nsfcst = nsfcstsave

      CALL cvttsnd( curtim, timsnd, tmstrln )

      soiloutfl = runname(1:lfnkey)//'.new_soilvar.'                    &
                  //timsnd(1:tmstrln)
      lfn = lfnkey + 13 + tmstrln

      IF( dirname /= ' ' ) THEN

        temchar = soiloutfl
        soiloutfl = dirname(1:ldirnam)//'/'//temchar
        lfn  = lfn + ldirnam + 1

      END IF

      CALL fnversn(soiloutfl, lfn)

      PRINT *, 'Write soil initial data to ',soiloutfl(1:lfn)

      CALL wrtsoil( nx,ny,nstyps, soiloutfl(1:lfn),dx,dy,               &
               mapproj,trulat1,trulat2,trulon,sclfct,ctrlat,ctrlon,     &
               1,1,1,1,1,1,                                             &
               tsfc,tsoil,wetsfc,wetdp,wetcanp,snowdpth,soiltyp)

      CALL soilcntl( nx,ny, soiloutfl, 1,1,1,1,1,1, x,y)

    END IF

  END IF
!
!-----------------------------------------------------------------------
!
!  Update the boundary zone fluxes.
!
!-----------------------------------------------------------------------
!
  IF (mp_opt > 0) THEN
    CALL acct_interrupt(mp_acct)
    CALL mpsend1dew(usflx,nx,ny,ebc,wbc,1,mptag1,temxpy)
    CALL mpsend1dew(vsflx,nx,ny,ebc,wbc,2,mptag2,temxpy)
    CALL mpsend1dew(ptsflx,nx,ny,ebc,wbc,0,mptag3,temxpy)
    CALL mpsend1dew(qvsflx,nx,ny,ebc,wbc,0,mptag4,temxpy)
    CALL mprecv1dew(usflx,nx,ny,ebc,wbc,1,mptag1,temxpy)
    CALL mprecv1dew(vsflx,nx,ny,ebc,wbc,2,mptag2,temxpy)
    CALL mprecv1dew(ptsflx,nx,ny,ebc,wbc,0,mptag3,temxpy)
    CALL mprecv1dew(qvsflx,nx,ny,ebc,wbc,0,mptag4,temxpy)
    CALL mpsend1dns(usflx,nx,ny,nbc,sbc,1,mptag1,temxpy)
    CALL mpsend1dns(vsflx,nx,ny,nbc,sbc,2,mptag2,temxpy)
    CALL mpsend1dns(ptsflx,nx,ny,nbc,sbc,0,mptag3,temxpy)
    CALL mpsend1dns(qvsflx,nx,ny,nbc,sbc,0,mptag4,temxpy)
    CALL mprecv1dns(usflx,nx,ny,nbc,sbc,1,mptag1,temxpy)
    CALL mprecv1dns(vsflx,nx,ny,nbc,sbc,2,mptag2,temxpy)
    CALL mprecv1dns(ptsflx,nx,ny,nbc,sbc,0,mptag3,temxpy)
    CALL mprecv1dns(qvsflx,nx,ny,nbc,sbc,0,mptag4,temxpy)
  END IF

  CALL acct_interrupt(bc_acct)
  CALL bcu2d(nx,ny, usflx,  ebc,wbc,nbc,sbc)
  CALL bcv2d(nx,ny, vsflx,  ebc,wbc,nbc,sbc)
  CALL bcs2d(nx,ny, ptsflx, ebc,wbc,nbc,sbc)
  CALL bcs2d(nx,ny, qvsflx, ebc,wbc,nbc,sbc)
  CALL acct_stop_inter
!
!-----------------------------------------------------------------------
!
!  Smooth fluxes to remove any 2-dx waves caused by discontinuities
!  in the soil conditions.
!
!-----------------------------------------------------------------------
!
  IF ( smthflx == 1 ) THEN
    DO nsmth=1,numsmth
      CALL smooth9p_nobc( usflx, nx,ny,1,nx,  1,ny-1,temxy1 )
      CALL smooth9p_nobc( vsflx, nx,ny,1,nx-1,1,ny,  temxy1 )
      CALL smooth9p_nobc( ptsflx,nx,ny,1,nx-1,1,ny-1,temxy1 )
      CALL smooth9p_nobc( qvsflx,nx,ny,1,nx-1,1,ny-1,temxy1 )

      IF (mp_opt > 0) THEN
        CALL acct_interrupt(mp_acct)
        CALL mpsend1dew(usflx,nx,ny,ebc,wbc,1,mptag1,temxpy)
        CALL mpsend1dew(vsflx,nx,ny,ebc,wbc,2,mptag2,temxpy)
        CALL mpsend1dew(ptsflx,nx,ny,ebc,wbc,0,mptag3,temxpy)
        CALL mpsend1dew(qvsflx,nx,ny,ebc,wbc,0,mptag4,temxpy)
        CALL mprecv1dew(usflx,nx,ny,ebc,wbc,1,mptag1,temxpy)
        CALL mprecv1dew(vsflx,nx,ny,ebc,wbc,2,mptag2,temxpy)
        CALL mprecv1dew(ptsflx,nx,ny,ebc,wbc,0,mptag3,temxpy)
        CALL mprecv1dew(qvsflx,nx,ny,ebc,wbc,0,mptag4,temxpy)
        CALL mpsend1dns(usflx,nx,ny,nbc,sbc,1,mptag1,temxpy)
        CALL mpsend1dns(vsflx,nx,ny,nbc,sbc,2,mptag2,temxpy)
        CALL mpsend1dns(ptsflx,nx,ny,nbc,sbc,0,mptag3,temxpy)
        CALL mpsend1dns(qvsflx,nx,ny,nbc,sbc,0,mptag4,temxpy)
        CALL mprecv1dns(usflx,nx,ny,nbc,sbc,1,mptag1,temxpy)
        CALL mprecv1dns(vsflx,nx,ny,nbc,sbc,2,mptag2,temxpy)
        CALL mprecv1dns(ptsflx,nx,ny,nbc,sbc,0,mptag3,temxpy)
        CALL mprecv1dns(qvsflx,nx,ny,nbc,sbc,0,mptag4,temxpy)
      END IF
      CALL acct_interrupt(bc_acct)
      CALL bcu2d(nx,ny, usflx,  ebc,wbc,nbc,sbc)
      CALL bcv2d(nx,ny, vsflx,  ebc,wbc,nbc,sbc)
      CALL bcs2d(nx,ny, ptsflx, ebc,wbc,nbc,sbc)
      CALL bcs2d(nx,ny, qvsflx, ebc,wbc,nbc,sbc)
      CALL acct_stop_inter
    END DO
  END IF
!
!-----------------------------------------------------------------------
!
!  Write all surface related fields for diagnostic purpose.
!
!-----------------------------------------------------------------------
!
  IF ( sfcdiag == 1 ) THEN

    CALL set_acct(output_acct)

    CALL soildiag(nx,ny,nz,x,y,z,                                       &
         soiltyp(1,1,1),vegtyp,lai,roufns,veg,zp(1,1,2),                &
         tsfc(1,1,0),tsoil(1,1,0),wetsfc(1,1,0),                        &
         wetdp(1,1,0),wetcanp(1,1,0),qvsfc(1,1,0),                      &
         usflx,vsflx,ptsflx,qvsflx,                                     &
         wdspsfc,psfc,rhosfc,prcrate,                                   &
         tair,qvair,cdh,cdq,cdm,                                        &
         radsw,rnflx,                                                   &
         shflx,lhflx,gflx,ct,evaprg,evaprtr,evaprr,qvsat,qvsata,f34)

  END IF
!
!-----------------------------------------------------------------------
!  Call the subroutine to determine height of PBL top
!-----------------------------------------------------------------------
!
  IF (pbldopt > 0) THEN

    CALL set_acct(sfcphy_acct)

    CALL pbldepth(nx,ny,nz, u,v,w,ptprt,qv,ptbar, zp,j3, ptsflx,        &
                  pbldpth)

  END IF

  RETURN
END SUBROUTINE sfcphysics
!
!##################################################################
!##################################################################
!######                                                      ######
!######                SUBROUTINE SFCFLX                     ######
!######                                                      ######
!######                     Developed by                     ######
!######     Center for Analysis and Prediction of Storms     ######
!######                University of Oklahoma                ######
!######                                                      ######
!##################################################################
!##################################################################
!


SUBROUTINE sfcflx(nx,ny,nz,                                             & 3
           u,v,w,ptprt,pprt,qv,ptbar,pbar,qvbar,                        &
           wdspsfc,rhosfc,ptsfc,qvsfc, cdm,cdh,cdq,                     &
           usflx,vsflx,ptsflx,qvsflx)
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  Calculate the surface momentum, heat and surface moisture fluxes.
!  These fluxes will be passed into turbulent mixing subroutines.
!
!-----------------------------------------------------------------------
!
!  AUTHOR: Ming Xue
!  4/07/1992.
!
!  MODIFICATION HISTORY:
!
!  2/24/1999 (M. Xue)
!  Simplified codes that calculate the sfc fluxes.
!
!-----------------------------------------------------------------------
!
!  INPUT:
!
!    nx       Number of grid points in the x-direction (east/west)
!    ny       Number of grid points in the y-direction (north/south)
!    nz       Number of grid points in the vertical
!
!    u        x component of velocity at a given time level (m/s)
!    v        y component of velocity at a given time level (m/s)
!    w        z component of velocity at a given time level (m/s)
!
!    ptprt    Perturbation potential temperature at a given time
!             level (K)
!    pprt     Perturbation pressure (Pascal)
!    qv       Water vapor specific humidity at a given time level
!             (kg/kg)
!    ptbar    Base state potential temperature (K)
!    pbar     Base state pressure (Pascal)
!    qvbar    Base state specific humidity (kg/kg)
!
!    wdspsfc  Wind speed at the surface
!    rhosfc   Base-state air density at the surface
!    ptsfc    Potential temperature at the ground level (K)
!    qvsfc    Effective specific humidity at sfc.
!    cdm      Drag coefficient for momentum defined as scalar point
!    cdh      Drag coefficient for heat
!    cdq      Drag coefficient for moisture
!
!  OUTPUT:
!
!    usflx    Surface flux of u-momentum
!    vsflx    Surface flux of v-momentum
!    ptsflx   Surface flux of heat (K*kg/(m**2*s))
!    qvsflx   Surface flux of moisture (K*kg/(m**2*s))
!
!-----------------------------------------------------------------------
!

!
!-----------------------------------------------------------------------
!
!  Variable Declarations.
!
!-----------------------------------------------------------------------
!
  IMPLICIT NONE

  INTEGER :: nx,ny,nz          ! The number grid points in 3 directions
  REAL :: u      (nx,ny,nz) ! Total u-velocity (m/s)
  REAL :: v      (nx,ny,nz) ! Total v-velocity (m/s)
  REAL :: w      (nx,ny,nz) ! Total w-velocity (m/s)

  REAL :: ptprt  (nx,ny,nz) ! Perturbation potential temperature (K)
  REAL :: pprt   (nx,ny,nz) ! Perturbation pressure (Pascal)
  REAL :: qv     (nx,ny,nz) ! Water vapor specific humidity (kg/kg)

  REAL :: ptbar  (nx,ny,nz) ! Base state potential temperature (K)
  REAL :: pbar   (nx,ny,nz) ! Base state pressure (Pascal)
  REAL :: qvbar  (nx,ny,nz) ! Base state specific humidity (kg/kg)

  REAL :: ptsfc  (nx,ny)    ! Potential temperature at the ground level (K)
  REAL :: rhosfc (nx,ny)    ! Surface air density (kg/m**3)
  REAL :: wdspsfc(nx,ny)    ! Wind speed at the surface
  REAL :: qvsfc  (nx,ny)    ! Effective specific humidity at surface

  REAL :: cdm    (nx,ny)    ! Drag coefficient for surface momentum flux
  REAL :: cdh    (nx,ny)    ! Drag coefficient for surface heat flux
  REAL :: cdq    (nx,ny)    ! Drag coefficient for surface moisture flux

  REAL :: usflx  (nx,ny)    ! surface flux of u-momentum (kg/(m*s**2))
  REAL :: vsflx  (nx,ny)    ! surface flux of v-momentum (kg/(m*s**2))
  REAL :: ptsflx (nx,ny)    ! surface flux of heat (K*kg/(m**2*s))
  REAL :: qvsflx (nx,ny)    ! surface flux of moisture (kg/(m**2*s))
!
!-----------------------------------------------------------------------
!
!  Misc. local variables:
!
!-----------------------------------------------------------------------
!
  INTEGER :: i,j
  REAL :: usuf, vsuf,qvsfc0
  REAL :: tema
  REAL :: gumove,gvmove
!
!-----------------------------------------------------------------------
!
!  Include files:
!
!  Include file globcst.inc passes cdmlnd,cdmwtr, cdhlnd,cdmwtr,
!  cdqlnd,cdmwtr and UMOVE and VMOVE into this subroutine.
!
!-----------------------------------------------------------------------
!
  INCLUDE 'globcst.inc'
  INCLUDE 'sfcphycst.inc'
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
  IF ( grdtrns == 0 ) THEN
    gumove = 0.0
    gvmove = 0.0
  ELSE
    gumove = umove
    gvmove = vmove
  END IF
!
!-----------------------------------------------------------------------
!
!  Compute the surface flux of u-momentum, i.e. the surface
!  drag in x-direction, according to a simple drag law:
!
!  usflx = rhobar *bar( (w'u') ) = rhobar *Cdm*V*u
!
!  where Cdm is the drag coefficient, V the wind speed at the first
!  scalar level above ground, and u the x-component of wind at the
!  same level.
!
!  usflx is defined one-half dz below the u-point, i.e. at ground level.
!
!-----------------------------------------------------------------------
!
  DO j=1,ny-1
    DO i=2,nx-1
      usuf = 0.5*(u(i,j,2)+u(i,j,1))+gumove
      usflx(i,j)=0.125*(cdm(i,j)+cdm(i-1,j))*                           &
                 (rhosfc(i,j)+rhosfc(i-1,j))*                           &
                 (wdspsfc(i,j)+wdspsfc(i-1,j))*                         &
                 sign(1.0, usuf)*MAX(abs(usuf),vsfcmin)
    END DO
  END DO
!
!-----------------------------------------------------------------------
!
!  Compute the surface flux of v-momentum, i.e. the surface
!  drag in y-direction, according to a simple drag law:
!
!  vsflx = - rhobar *bar( (w'v') ) = rhobar *Cdm*V*v
!
!  where Cdm is the drag coefficient, V the wind speed at the first
!  scalar level above ground, and v the y-component of wind at the
!  same level.
!
!  vsflx is defined one-half dz below the v-point, i.e. at ground level.
!
!-----------------------------------------------------------------------
!
  DO j=2,ny-1
    DO i=1,nx-1
      vsuf = gvmove + 0.5*(v(i,j,2)+v(i,j,1))
      vsflx(i,j)=0.125*(cdm(i,j)+cdm(i,j-1))*                           &
                 (rhosfc(i,j)+rhosfc(i,j-1))*                           &
                 (wdspsfc(i,j)+wdspsfc(i,j-1))*                         &
                 sign(1.0,vsuf)*MAX(abs(vsuf),vsfcmin)
    END DO
  END DO
!
!-----------------------------------------------------------------------
!
!  Compute the surface heat flux according to a simple drag law:
!
!  ptsflx = - rhobar *bar( (w' pt') ) = rhobar *Cdh*V*(pt-ptsfc0)
!
!  where Cdh is the drag coefficient, V the wind speed and pt the
!  potential temperatureat at the first scalar level above ground.
!
!  ptsfc0 is defined as the potential temperature at ground leval.
!  ptsfc0 = ptbar at surface
!
!  ptsflx is defined one-half dz below the scalar point, i.e. at
!  ground level.
!
!  Similarly, the vertical moisture flux at the surface is defined as
!
!  qvsflx = - rhobar *bar( (w' qv') ) = rhobar *Cdq*V*(qv-qvsfc0)
!
!  where qv is the specific humidity at the first model scalar level.
!  qvsfc0 is defined as the mixing ratio at the ground and taken as
!  the base state moisture, qvsfc0 = qvbar, at ground level.
!
!-----------------------------------------------------------------------
!
  DO j=1,ny-1
    DO i=1,nx-1
      IF ( sfcphy == 1 .OR. sfcphy == 3 ) THEN
        qvsfc0=0.5*(qvbar(i,j,1)+qvbar(i,j,2))
      ELSE
        qvsfc0=qvsfc(i,j)
      END IF
      tema = rhosfc(i,j)*MAX(wdspsfc(i,j),vsfcmin)
      ptsflx(i,j)=cdh(i,j)*tema*(ptprt(i,j,2)+ptbar(i,j,2)-ptsfc(i,j))
      qvsflx(i,j)=cdq(i,j)*tema*(qv(i,j,2)-qvsfc0)
    END DO
  END DO

  RETURN
END SUBROUTINE sfcflx
!
!##################################################################
!##################################################################
!######                                                      ######
!######               SUBROUTINE DRAGCOEF                    ######
!######                                                      ######
!######                     Developed by                     ######
!######     Center for Analysis and Prediction of Storms     ######
!######                University of Oklahoma                ######
!######                                                      ######
!##################################################################
!##################################################################
!


SUBROUTINE dragcoef(nx,ny,nz,                                           & 2,6
           zp, ptsfc,qvsfc,soiltyp,roufns,wdspsfc,pt1,                  &
           cdm,cdh,cdq)
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  Calculate the surface momentum, heat and surface moisture fluxes
!  using a stability dependent formulation.
!
!-----------------------------------------------------------------------
!
!  AUTHOR: Ming Xue
!  4/07/1992.
!
!  MODIFICATION HISTORY:
!
!  5/06/92 (M. Xue)
!  Added full documentation.
!
!  6/4/92 (M. Xue and H. Jin)
!  Further facelift.
!
!  9/14/1992 (M. Xue)
!  Different drag coefficients defined for momentum, temperature and
!  moisture. Definition of ptsfc and qvsfc changed. qvbar added to
!  the argument list.
!
!  9/28/1992 (M. Xue)
!  Sign correction in all the flux formulations.
!
!  8/28/1993 (M. Xue and Hao jin)
!  Modified the flux formulations considering the terrain effect.
!
!  9/10/1993 (V. Wong, X. Song and N. Lin)
!  Stability dependent formulation used for usflx, vsflx, and ptsflx.
!  The different flux formulations come from subroutines USTARC and
!  PTSTARC.
!
!  9/04/1994 (K. Brewster, V. Wong and X. Song)
!  Facelift.
!
!  12/7/1994 (Y. Liu)
!  Cleaned up an empty DO loop with a label 360.
!
!  01/28/1995 (V. Wong and X. Song)
!  Add the option of stability and roughness dependent surface layer
!  parameterization for both land and sea surfaces.
!
!  03/02/1995 (Y. Liu)
!  Added the minimum surface total wind speed to guarantee the basic
!  heat and moisture fluxes.
!
!  02/07/96 (V.Wong and X.Song)
!  Changed the calculation of roughness zo over the sea according to
!  Clarke(1970).  The new formulation is insensitive to the value of
!  the depth. of the first model layer above the ground.
!  Added kvwtr to denote the Von Karman constant over the sea.
!  Set a lower limiter, zolimit, for zo, and an upper limiter, z1limit,
!  for depth of the surface layer z1.
!
!  03/17/1997 (Yuhe Liu, Vince Wong and Jing Tian)
!  Added an option for constant cdh (cdq) over water.
!
!  05/29/97 (V. Wong and X. Tan)
!  Modified the formulation considering the height of the surface
!  layer z1 may equal zero.
!
!  02/25/1999 (M. Xue)
!  Isolated from original SFCFLXSD routine. The above listing is for
!  SFCFLXSD.
!
!  01/09/2002 (Y. Wang)
!  Changed the DO loop lower bound from i=2, j=2 to i=1, j=1 in 
!  the calculation of drag coefficient for momentum.
!
!-----------------------------------------------------------------------
!
!  INPUT:
!
!    nx       Number of grid points in the x-direction (east/west)
!    ny       Number of grid points in the y-direction (north/south)
!    nz       Number of grid points in the vertical
!
!    zp       The physical height coordinate defined at w-point of
!             staggered grid.
!
!    ptsfc    Potential temperature at the ground level (K)
!    qvsfc    Effective specific humidity at sfc
!    soiltyp  Soil type at each point
!    roufns   Surface roughness
!
!    wdspsfc  Surface wind speed (at first level about ground) (m/s)
!    pt1      Potential temperature at the first model level AGL
!
!  OUTPUT:
!
!    cdm      Drag coefficient for surface momentum flux (nondimensional)
!    cdh      Drag coefficient for surface heat flux (nondimensional)
!    cdq      Drag coefficient for surface moisture flux (nondimensional)
!
!  Local work arrays:
!
!    c_u      Cu
!    c_pt     Cpt
!    c_uneu   Cu for neutral stratification
!    c_ptneu  Cpt for neutral stratification
!
!-----------------------------------------------------------------------
!
!  Variable Declarations.
!
!-----------------------------------------------------------------------
!
  IMPLICIT NONE

  INTEGER :: nx,ny,nz          ! The number grid points in 3 directions

  REAL :: zp    (nx,ny,nz)     ! The physical height coordinate defined at
                               ! w-point of staggered grid.
  REAL :: ptsfc (nx,ny)        ! Potential temperature at the ground level (K)
  REAL :: qvsfc (nx,ny)        ! Effective surface specific humidity

  INTEGER :: soiltyp(nx,ny)    ! Soil type at each point
  REAL :: roufns(nx,ny)        ! Surface roughness length

  REAL :: wdspsfc(nx,ny)       ! Wind speed at the surface
  REAL :: pt1   (nx,ny)        ! Potential temperature at the first model level AGL

  REAL :: cdm   (nx,ny)        ! Drag coefficient for surface momentum flux
  REAL :: cdh   (nx,ny)        ! Drag coefficient for surface heat flux
  REAL :: cdq   (nx,ny)        ! Drag coefficient for surface moisture flux

  REAL :: c_u   (nx,ny)        ! Cu
  REAL :: c_pt  (nx,ny)        ! Cpt
  REAL :: c_uneu(nx,ny)        ! Cu for neutral stratification
  REAL :: c_ptneu(nx,ny)       ! Cpt for neutral stratification
!
!-----------------------------------------------------------------------
!
!  Misc. local variables:
!
!-----------------------------------------------------------------------
!
  INTEGER :: i, j
  REAL :: check
  REAL :: z1,z1drou,z1droup
!
!-----------------------------------------------------------------------
!
!  Include files:
!
!  Include file globcst.inc passes cdq (moisture transfer coeffcient)
!  and UMOVE and VMOVE into this subroutine.
!
!-----------------------------------------------------------------------
!
  INCLUDE 'globcst.inc'
  INCLUDE 'phycst.inc'
  INCLUDE 'sfcphycst.inc'
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
  DO j=1,ny-1
    DO i=1,nx-1
!
! Calculate C_u and C_pt at the neutral state
!
      z1  = 0.5*(zp(i,j,3)-zp(i,j,2))
      z1  = MIN(z1,z1limit)
      z1drou = z1/roufns(i,j)

      z1droup = (z1+roufns(i,j))/roufns(i,j)

      c_uneu (i,j) = kv/LOG(z1droup)
      c_ptneu(i,j) = c_uneu(i,j) /prantl0l


    END DO
  END DO
!
!-----------------------------------------------------------------------
!
!  A general procedure for computing the drag coefficients:
!
!  1.  Obtain the following quantities:
!      a. The first scalar level above ground level (z1)
!      b. Potential temperature at the height of surface roughness
!         (roufns) - given.
!      c. Potential temperature (pt1) at the first scalar level above
!         ground.
!      d. C_u and C_pt at neutral state in the land case from
!             c_uneu =kv/log10(z1/roughness)
!             c_ptneu=c_uneu/prandtl_number
!         In the water case, call subroutines CUNEUWTR and CPTNEUWTR.
!
!  2.  Given the information from (1), compute the Bulk Richardson
!      number, defined as
!
!         BulkRi = (g/ptsfc) * ( (pt1-ptsfc)*(z1-roufns) / V**2 )
!      where V is the wind speed at the surface.
!
!  3.  After calculation of the BulkRi, then:
!      a. Given the BulkRi and the information in (1), compute C_u at
!         both land and water cases (see subroutines CUC and CUCWTR).
!
!      b. From C_u and C_pt, compute draf coefficients cdh and cdq.
!
!-----------------------------------------------------------------------
!
!  Calculate drag coefficient (defined at scalar point) for momentum
!  fluxes.
!
!-----------------------------------------------------------------------
!
  CALL cuc(nx,ny,nz,1,nx-1,1,ny-1,zp,                                   &
           roufns,wdspsfc,ptsfc,pt1,c_uneu, c_u)

  IF (wtrexist == 1) THEN

    CALL cuneuwtr(nx,ny,nz,1,nx-1,1,ny-1,soiltyp,wdspsfc,               &
                  c_uneu)

    CALL cucwtr  (nx,ny,nz,1,nx-1,1,ny-1,zp,soiltyp,wdspsfc,            &
                  ptsfc,pt1,c_uneu, c_u)

  END IF

  DO j=1,ny-1
    DO i=1,nx-1
      check=ptsfc(i,j)-pt1(i,j)
!
!    Put lower and upper limits on C_u
!
      IF(check >= 0.0) THEN  ! Unstable case
        c_u(i,j)=MIN( c_u(i,j), 2.0*c_uneu(i,j))
      ELSE                   ! Stable case
        c_u(i,j)=MAX( c_u(i,j), blimit*c_uneu(i,j))
      END IF
!
!-----------------------------------------------------------------------
!
!  c_dm=c_U**2
!
!  Change Von Karman constant from 0.4 to 0.35 and since the drag
!  coefficient Cdm is proportional to the square of the constant,
!  the final usflx should be multiplied by (0.35/0.40)**2 = 0.765625
!
!-----------------------------------------------------------------------
!
      cdm(i,j) = 0.765625 * c_u(i,j) * c_u(i,j)

    END DO
  END DO
!
!-----------------------------------------------------------------------
!
!  From C_u and C_pt, compute C_dh = C_u * C_pt and pt flux,
!
!-----------------------------------------------------------------------
!
  CALL cptc(nx,ny,nz,1,nx-1,1,ny-1,zp,roufns,wdspsfc,                   &
            ptsfc,pt1,c_ptneu,c_pt)

  IF ( wtrexist == 1 .AND. cdhwtropt == 0 ) THEN

    CALL cptneuwtr(nx,ny,nz,1,nx-1,1,ny-1,soiltyp,wdspsfc,c_uneu,       &
                   c_ptneu)

    CALL cptcwtr(nx,ny,nz,1,nx-1,1,ny-1,zp,soiltyp,wdspsfc,             &
                 ptsfc,pt1,c_ptneu, c_pt)

  END IF

  DO j=1,ny-1
    DO i=1,nx-1
      check=ptsfc(i,j)-pt1(i,j)
!
!    Impose lower and upper limits on C_u and C_pt.
!
      IF ( soiltyp(i,j) == 13 .AND. cdhwtropt == 1 ) THEN
        cdh(i,j) = cdhwtr
      ELSE
        IF(check >= 0) THEN
          c_pt(i,j) = MIN(c_pt(i,j),3.3333*c_ptneu(i,j) )
        ELSE
          c_pt(i,j) = MAX(c_pt(i,j),blimit*c_ptneu(i,j) )
        END IF

        cdh(i,j) = c_u(i,j)*c_pt(i,j)
      END IF
!
!-----------------------------------------------------------------------
!
!  Change Von Karman constant from 0.4 to 0.35 and since the drag
!  coefficient Cdh is proportional to the square of the constant,
!  the final Cdh should be multiplied by (0.35/0.40)**2 = 0.765625
!
!-----------------------------------------------------------------------
!
      cdh(i,j) = 0.765625*cdh(i,j)
      cdq(i,j) = 0.7*cdh(i,j)

    END DO
  END DO

  RETURN
END SUBROUTINE dragcoef
!
!##################################################################
!##################################################################
!######                                                      ######
!######                    SUBROUTINE CUC                    ######
!######                                                      ######
!######                     Developed by                     ######
!######     Center for Analysis and Prediction of Storms     ######
!######                University of Oklahoma                ######
!######                                                      ######
!##################################################################
!##################################################################
!


SUBROUTINE cuc(nx,ny,nz,ibgn,iend,jbgn,jend,zp,roufns,wspd,ptsfc,       & 1
           pt1, c_uneu,c_u)
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  Compute C_u (friction velocity / wind speed). The quantity C_U is used by
!  the subroutine SFCFLX to obtain surface fluxes for both unstable
!  and stable cases.
!
!-----------------------------------------------------------------------
!
!  AUTHORS: V. Wong, X. Song and N. Lin
!  9/10/1993
!
!  For the unstable case (Bulk Richardson number bulkri < 0 ), the
!  formulation is based on the paper "On the Analytical Solutions of
!  Flux-Profile relationships for the Atmospheric Surface Layer" by
!  D.W. Byun in J. of Applied Meteorlogy, July 1990, pp. 652-657.
!
!  For the stable case, the formulation is based on the
!  paper "A Short History of the Operational PBL - Parameterization
!  at ECMWF" by J.F.Louis, M. Tiedtke and J.F. Geleyn in " Workshop
!  on Planetary Boundary Layer Parameterization", a publication
!  by the European Centre for Medium Range Weather Forecasts,
!  25-27 Nov. 1981.
!
!  MODIFICATION HISTORY:
!
!  9/04/1994 (K. Brewster, V. Wong and X. Song)
!  Facelift.
!
!  2/27/95 (V. Wong and X. Song)
!
!  02/07/96 (V.Wong and X.Song)
!  Set an upper limiter, z1limit, for depth of the surface layer z1.
!
!  05/01/97 (V. Wong and X. Tan)
!  Changed the computation of stabp, use the formulation similar to
!  one used by Bynn(1990), the momentum and thermal roughness
!  lengths are different.
!
!  05/29/97 (V. Wong and X. Tan)
!  Modified the formulation considering the height of the surface
!  layer z1 may equal zero.
!
!-----------------------------------------------------------------------
!
!  INPUT:
!
!    nx       Number of grid points in the x-direction (east/west)
!    ny       Number of grid points in the y-direction (north/south)
!    nz       Number of grid points in the vertical
!
!    ibgn     i-index where evaluation begins.
!    iend     i-index where evaluation ends.
!    jbgn     j-index where evaluation begins.
!    jend     j-index where evaluation ends.
!
!    zp       The physical height coordinate defined at w-point of
!             staggered grid.
!    roufns   Surface roughness
!
!    wspd     Surface wind speed (m/s), defined as
!             sqrt(usuf*usuf + vsuf*vsuf + wsuf*wsuf)
!    ptsfc    Potential temperature at the ground level (K)
!    pt1      Potential temperature at the 1st scalar point above
!             ground level, (K)
!
!    c_uneu   Friction velocity at neutral state
!
!  OUTPUT:
!
!    ustar    Friction velocity
!
!-----------------------------------------------------------------------
!

!
!-----------------------------------------------------------------------
!
!  Variable Declarations.
!
!-----------------------------------------------------------------------
!
  IMPLICIT NONE

  INTEGER :: nx,ny,nz          ! The number grid points in 3 directions

  INTEGER :: ibgn              ! i-index where evaluation begins.
  INTEGER :: iend              ! i-index where evaluation ends.
  INTEGER :: jbgn              ! j-index where evaluation begins.
  INTEGER :: jend              ! j-index where evaluation ends.

  REAL :: zp    (nx,ny,nz)     ! The physical height coordinate
                               ! defined at w-point of staggered grid.
  REAL :: roufns(nx,ny)        ! Surface roughness length

  REAL :: wspd  (nx,ny)        ! Surface wind speed (m/s)

  REAL :: ptsfc(nx,ny)         ! Potential temperature at the ground
                               ! level (K)
  REAL :: pt1   (nx,ny)        ! Potential temperature at the
                               ! 1st scalar
                               ! point above ground level, (K)
  REAL :: c_uneu(nx,ny)      ! Frictional velocity (m/s) at neutral

  REAL :: c_u(nx,ny)        ! Frictional velocity (m/s)
!
!-----------------------------------------------------------------------
!
!  Misc. local variables:
!
!-----------------------------------------------------------------------
!
  INTEGER :: i, j

  REAL :: z1                   ! The height of 1st scalar point above
                               ! ground level (m)
  REAL :: dz0                  ! z1-roufns, where roufns is defined as
                               ! surface roughness length
  REAL :: bulkri               ! Bulk Richardson number

  REAL :: stabp                ! Monin-Obukhov STABility Parameter
                               ! (zeta)
  REAL :: x1,x0,psim           ! Intermediate parameters needed
  REAL :: z1drou,qb3pb2
  REAL :: c7,c8
  REAL :: sb,qb,pb,thetab,tb   ! During  computations
  REAL :: a,b,c,d
  REAL :: tempan,sqrtqb

  REAL :: zt                   ! Thermal roughness length
  REAL :: dzt
  REAL :: ztdrou

  REAL :: z1droup, ztdroup
!
!-----------------------------------------------------------------------
!
!  Include files:
!
!-----------------------------------------------------------------------
!
  INCLUDE 'phycst.inc'
  INCLUDE 'globcst.inc'
  INCLUDE 'sfcphycst.inc'
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!
!-----------------------------------------------------------------------
!
!  Following Byun (1990).
!
!-----------------------------------------------------------------------
!
  DO j=jbgn,jend
    DO i=ibgn,iend

      z1  = 0.5*(zp(i,j,3)-zp(i,j,2))
      z1  = MIN(z1,z1limit)

      zt = ztz0*roufns(i,j)
      dzt = z1-zt
      ztdrou = z1/zt
      ztdroup = (z1+zt)/zt

      dz0 = z1-roufns(i,j)
      z1drou = z1/roufns(i,j)
      z1droup = (z1+roufns(i,j))/roufns(i,j)

      bulkri = (g/ptsfc(i,j))*(pt1(i,j)-ptsfc(i,j))*dz0/                &
               (wspd(i,j)*wspd(i,j))

      IF (bulkri <= 0.0) THEN
!
!-----------------------------------------------------------------------
!
!  Unstable case: See equations (28)-(34) in Byun (1990).
!
!-----------------------------------------------------------------------
!
        bulkri = MAX (bulkri,-10.0)

        sb =bulkri/prantl0l

        qb=oned9*(c1l+c2l*sb*sb)
        pb=oned54*(c3l+c4l*sb*sb)

        qb3pb2=qb**3-pb*pb
        c7 = (z1*dzt*LOG(z1droup)*LOG(z1droup))/(dz0*dz0*LOG(ztdroup))

        IF( qb3pb2 >= 0.0 ) THEN

          sqrtqb = SQRT(qb)
          tempan = MAX( -1.0, MIN( 1.0, pb/(sqrtqb**3) ) )

          thetab=ACOS(tempan)
          stabp =c7*(-2.0*sqrtqb*COS(thetab/3.0)+c5l)

        ELSE

          tb    =(SQRT(-qb3pb2)+ABS(pb))**oned3
          stabp =c7*(-(tb+qb/tb)+c5l)

        END IF
!
!-----------------------------------------------------------------------
!
!  According to equation (14) in Byun (1990).
!
!-----------------------------------------------------------------------
!
        c8=gammaml*stabp
        x1=(1. - c8)**0.25
        x0=(1. - c8/z1drou)**0.25

        psim=2.0*LOG((1.0+x1)/(1.0+x0))+LOG((1.+x1*x1)/(1.+x0*x0))-     &
             2.0*ATAN(x1)+2.0*ATAN(x0)

!
!-----------------------------------------------------------------------
!
!  Compute C_u via equation (10) in Byun (1981).
!
!-----------------------------------------------------------------------
!
        c_u(i,j) =kv/(LOG(z1droup)-psim)

      ELSE
!
!-----------------------------------------------------------------------
!
!  Stable case:
!
!-----------------------------------------------------------------------
!
        a=kv/LOG(z1droup)
        b=5.0
        d=5.0
        c=SQRT(1.0+d*bulkri)

        c_u(i,j) = a/SQRT(1.0+2.0*b*bulkri/c)

      END IF

    END DO
  END DO

  RETURN
END SUBROUTINE cuc
!
!##################################################################
!##################################################################
!######                                                      ######
!######                   SUBROUTINE CPTC                    ######
!######                                                      ######
!######                     Developed by                     ######
!######     Center for Analysis and Prediction of Storms     ######
!######                University of Oklahoma                ######
!######                                                      ######
!##################################################################
!##################################################################
!


SUBROUTINE cptc(nx,ny,nz,ibgn,iend,jbgn,jend,                           & 1
           zp,roufns,wspd,ptsfc,pt1, c_ptneu,c_pt)
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  Compute c_pt (a nondimensional temperature scale)
!
!  The quantity c_pt is used by the subroutine SFCFLX to obtain
!  surface fluxes for both the unstable and stable cases.
!
!-----------------------------------------------------------------------
!
!  AUTHORS: V. Wong, X. Song and N. Lin
!  9/10/1993
!
!  For the unstable case (Bulk Richardson number bulkri < 0 ), the
!  formulation is based on the paper "On the Analytical Solutions of
!  Flux-Profile relationships for the Atmospheric Surface Layer" by
!  D.W. Byun in J. of Applied Meteorlogy, July 1990, pp. 652-657.
!
!  For stable case, the formulation is based on the paper "A Short
!  History of the Operational PBL - Parameterization at ECMWF" by
!  J.F.Louis, M. Tiedtke and J.F. Geleyn in "Workshop on Planetary
!  Boundary Layer Parameterization", a publication by European Centre
!  for Medium Range Weather Forecasts, 25-27 Nov. 1981.
!
!  MODIFICATION HISTORY:
!
!  9/04/1994 (K. Brewster, V. Wong and X. Song)
!  Facelift.
!
!  2/27/95 (V. Wong and X. Song)
!
!  02/07/96 (V.Wong and X.Song)
!  Set an upper limiter, z1limit, for depth of the surface layer z1.
!
!  05/01/97 (V. Wong and X. Tan)
!  Changed the computation of temperature scale
!
!  05/29/97 (V. Wong and X. Tan)
!  Modified the formulation considering the height of the surface
!  layer z1 may equal zero.
!
!-----------------------------------------------------------------------
!
!  INPUT:
!
!    nx       Number of grid points in the x-direction (east/west)
!    ny       Number of grid points in the y-direction (north/south)
!    nz       Number of grid points in the vertical
!
!    ibgn     i-index where evaluation begins.
!    iend     i-index where evaluation ends.
!    jbgn     j-index where evaluation begins.
!    jend     j-index where evaluation ends.
!
!    zp       The physical height coordinate defined at w-point of
!             staggered grid.
!    wspd     Surface wind speed (m/s), defined as
!             sqrt(usuf*usuf + vsuf*vsuf + wsuf*wsuf)
!    roufns   Surface roughness
!
!    ptsfc    Potential temperature at the ground level (K)
!    pt1      Potential temperature at the 1st scalar point above
!             ground level, (K)
!
!    c_ptneu  Temperature scale (K) at neutral state, defined by
!             surface heat flux / friction velocity
!
!  OUTPUT:
!
!    c_pt     Temperature scale (K), defined by
!             surface heat flux / friction velocity
!
!-----------------------------------------------------------------------
!

!
!-----------------------------------------------------------------------
!
!  Variable Declarations.
!
!-----------------------------------------------------------------------
!
  IMPLICIT NONE

  INTEGER :: nx,ny,nz          ! The number grid points in 3 directions

  INTEGER :: ibgn              ! i-index where evaluation begins
  INTEGER :: iend              ! i-index where evaluation ends
  INTEGER :: jbgn              ! j-index where evaluation begins
  INTEGER :: jend              ! j-index where evaluation ends

  REAL :: zp    (nx,ny,nz)     ! The physical height coordinate
                               ! defined at w-point of staggered grid.
  REAL :: roufns(nx,ny)        ! Surface roughness length

  REAL :: wspd  (nx,ny)        ! Surface wind speed (m/s)

  REAL :: ptsfc(nx,ny)         ! Potential temperature at the
                               ! ground level (K)
  REAL :: pt1   (nx,ny)        ! Potential temperature at the
                               ! 1st scalar
                               ! point above ground level, (K)
  REAL :: c_ptneu(nx,ny)       ! Normalized temperature scale (K)
                               ! at neutral state

  REAL :: c_pt  (nx,ny)        ! Temperature scale (K)
!
!-----------------------------------------------------------------------
!
!  Misc. local variables:
!
!-----------------------------------------------------------------------
!
  INTEGER :: i, j

  REAL :: z1                  ! The height of 1st scalar point above
                              ! ground level (m)
  REAL :: dz0                 ! z1-roufns, where roufns is defined as
                              ! surface roughness length
  REAL :: bulkri              ! Bulk Richardson number

  REAL :: stabp               ! Monin-Obukhov STABility Parameter
                              ! (zeta)

  REAL :: y1,y0,psih          ! Intermediate parameters needed
  REAL :: z1drou,qb3pb2
  REAL :: c7,c8
  REAL :: sb,qb,pb,thetab,tb  ! During  computations
  REAL :: a,b,c,d
  REAL :: tempan,sqrtqb

  REAL :: zt                  ! Thermal roughness length
  REAL :: dzt,ztdrou

  REAL :: z1droup,ztdroup
!
!-----------------------------------------------------------------------
!
!  Include files:
!
!-----------------------------------------------------------------------
!
  INCLUDE 'phycst.inc'
  INCLUDE 'globcst.inc'
  INCLUDE 'sfcphycst.inc'
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!
!-----------------------------------------------------------------------
!
!  Following Byun (1990).
!
!-----------------------------------------------------------------------
!
  DO j=jbgn,jend
    DO i=ibgn,iend

      z1  = 0.5*(zp(i,j,3)-zp(i,j,2))
      z1  = MIN(z1,z1limit)

      zt = ztz0*roufns(i,j)
      dzt = z1-zt
      ztdrou = z1/zt
      ztdroup = (z1+zt)/zt

      dz0 = z1-roufns(i,j)
      z1drou = z1/roufns(i,j)
      z1droup = (z1+roufns(i,j))/roufns(i,j)

      bulkri = (g/ptsfc(i,j))*(pt1(i,j)-ptsfc(i,j))*dz0/                &
               (wspd(i,j)*wspd(i,j))

      IF (bulkri <= 0.0) THEN
!
!-----------------------------------------------------------------------
!
!  Unstable case: See equations (28)-(34) in Byun (1990).
!
!-----------------------------------------------------------------------
!
        bulkri = MAX (bulkri,-10.0)

        sb =bulkri/prantl0l

        qb=oned9*(c1l+c2l*sb*sb)
        pb=oned54*(c3l+c4l*sb*sb)

        qb3pb2=qb**3-pb*pb
        c7 = (z1*dzt*LOG(z1droup)*LOG(z1droup))/(dz0*dz0*LOG(ztdroup))

        IF( qb3pb2 >= 0.0 ) THEN

          sqrtqb = SQRT(qb)
          tempan = MAX( -1.0, MIN( 1.0, pb/(sqrtqb**3) ) )

          thetab=ACOS(tempan)
          stabp =c7*(-2.0*SQRT(qb)*COS(thetab/3.0)+c5l)

        ELSE

          tb    =(SQRT(-qb3pb2)+ABS(pb))**oned3
          stabp =c7*(-(tb+qb/tb)+c5l)

        END IF
!
!-----------------------------------------------------------------------
!
!  According to equation (15) in Byun (1990).
!
!-----------------------------------------------------------------------
!
        c8=gammahl*stabp
        y1=SQRT(1.0 - c8)
        y0=SQRT(1.0 - c8/ztdrou)

        psih=2.0*LOG((y1+1.0)/(y0+1.0))
!
!-----------------------------------------------------------------------
!
!  Compute c_pt via equation (11) in Byun (1981).
!
!-----------------------------------------------------------------------
!
        c_pt(i,j)=kv / (prantl0l*(LOG(ztdroup)-psih))

      ELSE
!
!-----------------------------------------------------------------------
!
!  Stable case: See Louis et al (1981).
!
!-----------------------------------------------------------------------
!
        a=kv/LOG(ztdroup)
        b=5.0
        d=5.0
        c=SQRT(1.0+d*bulkri)

        c_pt(i,j) = SQRT(1.0+2.0*b*bulkri/c)
        c_pt(i,j) = a*c_pt(i,j)/(prantl0l*(1.0+3.0*b*bulkri*c))

      END IF

    END DO
  END DO

  RETURN
END SUBROUTINE cptc
!
!##################################################################
!##################################################################
!######                                                      ######
!######                 SUBROUTINE CUCWTR                    ######
!######                                                      ######
!######                     Developed by                     ######
!######     Center for Analysis and Prediction of Storms     ######
!######                University of Oklahoma                ######
!######                                                      ######
!##################################################################
!##################################################################
!


SUBROUTINE cucwtr(nx,ny,nz,ibgn,iend,jbgn,jend,zp,soiltyp,              & 1
           wspd,ptsfc,pt1,c_uwtrneu,c_uwtr)
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  Compute C_u (friction velocity / wind speed) at sea case. The quantity
!  C_uwtr is used by the subroutine SFCFLX to obtain surface fluxes for both
!  unstable and stable cases.
!
!-----------------------------------------------------------------------
!
!  AUTHORS: V. Wong and X. Song
!  8/04/1994
!
!  MODIFICATION HISTORY:
!
!  2/27/95 (V.W. and X.S.)
!
!  1/12/96 (V.W. and X.S.)
!  Changed the calculation related to zo.
!  Added kvwtr to denote the Von Karman constant over the sea;
!  Set a lower limiter for zo, zolimit, and an upper limiter for z1, z1limit.
!
!-----------------------------------------------------------------------
!
!  INPUT:
!
!    nx       Number of grid points in the x-direction (east/west)
!    ny       Number of grid points in the y-direction (north/south)
!    nz       Number of grid points in the vertical
!
!    ibgn     i-index where evaluation begins.
!    iend     i-index where evaluation ends.
!    jbgn     j-index where evaluation begins.
!    jend     j-index where evaluation ends.
!
!    zp       The physical height coordinate defined at w-point of
!             staggered grid.
!    wspd     Surface wind speed (m/s), defined as
!             sqrt(usuf*usuf + vsuf*vsuf + wsuf*wsuf)
!    ptsfc    Potential temperature at the ground level (K)
!    pt1      Potential temperature at the 1st scalar point above
!             ground level, (K)
!
!  OUTPUT:
!
!    c_uwtr   Friction velocity
!
!-----------------------------------------------------------------------
!

!
!-----------------------------------------------------------------------
!
!  Variable Declarations.
!
!-----------------------------------------------------------------------
!
  IMPLICIT NONE

  INTEGER :: nx,ny,nz          ! The number grid points in 3 directions

  INTEGER :: ibgn              ! i-index where evaluation begins.
  INTEGER :: iend              ! i-index where evaluation ends.
  INTEGER :: jbgn              ! j-index where evaluation begins.
  INTEGER :: jend              ! j-index where evaluation ends.

  REAL :: zp    (nx,ny,nz)     ! The physical height coordinate defined at
                               ! w-point of staggered grid.
  INTEGER :: soiltyp(nx,ny)    ! Soil type at each point.

  REAL :: wspd  (nx,ny)        ! Surface wind speed (m/s)

  REAL :: ptsfc(nx,ny)         ! Potential temperature at the ground level
                               ! (K)
  REAL :: pt1   (nx,ny)        ! Potential temperature at the 1st scalar
                               ! point above ground level, (K)
  REAL :: c_uwtrneu(nx,ny)        ! Frictional velocity (m/s)
  REAL :: c_uwtr(nx,ny)        ! Frictional velocity (m/s)
!
!-----------------------------------------------------------------------
!
!  Misc. local variables:
!
!-----------------------------------------------------------------------
!
  INTEGER :: i, j
  REAL :: z1                   ! The height of 1st scalar point above
                               ! ground level (m)
  REAL :: zo                   ! Defined as surface momentum roughness
                               ! length
  REAL :: zt                   ! Defined as surface heat transfer
                               ! roughness length
  REAL :: dzo                  ! z1-zo
  REAL :: dzt                  ! z1-zt
  REAL :: z1dzo                ! z1/zo
  REAL :: z1dzt                ! z1/zt
  REAL :: xcdn                 ! Cdn (sea)
  REAL :: xcdh                 ! Hot-wired value of Cdh (sea)
  REAL :: bulkri               ! Bulk Richardson number

  REAL :: stabp                ! Monin-Obukhov STABility Parameter
                               ! (zeta)

  REAL :: x1,psim           ! Intermediate parameters needed
  REAL :: qb3pb2
  REAL :: c7,c8
  REAL :: sb,qb,pb,thetab,tb   ! during  computations
  REAL :: a,b,c,d
  REAL :: tempan,sqrtqb
!
!-----------------------------------------------------------------------
!
!  Include files:
!
!-----------------------------------------------------------------------
!
  INCLUDE 'phycst.inc'
  INCLUDE 'globcst.inc'
  INCLUDE 'sfcphycst.inc'
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
  xcdh = 1.1E-3

  DO j=jbgn,jend
    DO i=ibgn,iend

      IF (soiltyp (i,j) == 13) THEN

        xcdn = (0.4+0.079*wspd(i,j)) * 1.e-3
        z1  = 0.5*(zp(i,j,3)-zp(i,j,2))
        z1  = MIN(z1,z1limit)
!
! calculate zo and zt
!
        zo =0.032*xcdn*wspd(i,j)*wspd(i,j)/9.8
        zo = MAX(zo,zolimit)
        zt = z1 * EXP( -kvwtr*SQRT(xcdn)/(prantl0w*xcdh) )

        dzo = z1-zo
        z1dzo = z1/zo
        dzt = z1-zt
        z1dzt = z1/zt

        bulkri = (g/ptsfc(i,j))*(pt1(i,j)-ptsfc(i,j))*dzo*dzo/          &
                 (dzt*wspd(i,j)*wspd(i,j))

        IF (bulkri <= 0.0) THEN
!
!-----------------------------------------------------------------------
!
!  Unstable case: A modified formulation, which is similar to
!  equations (28)-(34) in Byun (1990).
!
!-----------------------------------------------------------------------
!
          bulkri = MAX (bulkri,-10.0)

          sb =bulkri/prantl0w

          qb=oned9*(c1w+c2w*sb*sb)
          pb=oned54*(c3w+c4w*sb*sb)

          qb3pb2=qb**3-pb*pb
          c7 = (z1*dzt/(dzo*dzo))*(ALOG(z1dzo)*ALOG(z1dzo)/ALOG(z1dzt))

          IF( qb3pb2 >= 0.0 ) THEN

            sqrtqb = SQRT(qb)
            tempan = MAX( -1.0, MIN( 1.0, pb/(sqrtqb**3) ) )

            thetab=ACOS(tempan)
            stabp =c7*(-2.0*SQRT(qb)*COS(thetab/3.0)+c5w)

          ELSE

            tb    =(SQRT(-qb3pb2)+ABS(pb))**oned3
            stabp =c7*(-(tb+qb/tb)+c5w)

          END IF
!
!-----------------------------------------------------------------------
!
!  According to a modified equation, which is similar to equation (14)
!  in Byun (1990).
!
!-----------------------------------------------------------------------
!
          c8=gammamw*stabp
          x1=(1. - c8)**0.25

          psim=2.0*ALOG(0.5*(1.0+x1))+ALOG(0.5*(1.+x1*x1))-             &
               2.0*ATAN(x1)+ASIN(1.)

!
!-----------------------------------------------------------------------
!
!  Compute c_uwtr via equation (10) in Byun (1981).
!
!-----------------------------------------------------------------------
!
          c_uwtr(i,j) =kv/(LOG(z1dzo)-psim)

        ELSE
!
!-----------------------------------------------------------------------
!
!  Stable case: See Louis et al (1981).
!
!-----------------------------------------------------------------------
!
          a=kv/LOG(z1dzo)
          b=5.0
          d=5.0
          c=SQRT(1.0+d*bulkri)

          c_uwtr(i,j) = a/SQRT(1.0+2.0*b*bulkri/c)

        END IF

      END IF

    END DO
  END DO

  RETURN
END SUBROUTINE cucwtr

!
!##################################################################
!##################################################################
!######                                                      ######
!######                SUBROUTINE CPTCWTR                    ######
!######                                                      ######
!######                     Developed by                     ######
!######     Center for Analysis and Prediction of Storms     ######
!######                University of Oklahoma                ######
!######                                                      ######
!##################################################################
!##################################################################
!


SUBROUTINE cptcwtr(nx,ny,nz,ibgn,iend,jbgn,jend,zp,soiltyp,wspd,        & 1
           ptsfc,pt1,c_ptwtrneu,c_ptwtr)
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  Compute c_ptwtr (the product of ustar and ptstar) at sea case.
!  Note:  a temperature scale defined as surface heat flux divided
!  by the friction velocity.
!
!  The quantity c_ptwtr is used by the subroutine SFCFLX to obtain
!  surface fluxes for both the unstable and stable cases.
!
!-----------------------------------------------------------------------
!
!  AUTHORS: V. Wong and X. Song
!  8/04/1994
!
!  MODIFICATION HISTORY:
!
!  2/27/95 (V.W. and X.S.)
!
!  1/12/96 (V.W. and X.S.)
!  Changed the calculation related to zo over the sea.
!  Added kvwtr to denote the Von Karman constant over the sea;
!  Set a lower limiter for zo, zolimit, and an upper limiter for z1, z1limit.
!
!-----------------------------------------------------------------------
!
!  INPUT:
!
!    nx       Number of grid points in the x-direction (east/west)
!    ny       Number of grid points in the y-direction (north/south)
!    nz       Number of grid points in the vertical
!
!    iend     i-index where evaluation ends.
!    jend     j-index where evaluation ends.
!
!    zp       The physical height coordinate defined at w-point of
!             staggered grid.
!    wspd     Surface wind speed (m/s), defined as
!             sqrt(usuf*usuf + vsuf*vsuf + wsuf*wsuf)
!    ptsfc    Potential temperature at the ground level (K)
!    pt1      Potential temperature at the 1st scalar point above
!             ground level, (K)
!
!
!  OUTPUT:
!
!    c_ptwtr  The product of ustar and ptstar. ptstar is temperature
!             scale (K), defined by surface heat flux / friction velocity
!
!-----------------------------------------------------------------------
!
!
!-----------------------------------------------------------------------
!
!  Variable Declarations.
!
!-----------------------------------------------------------------------
!
  IMPLICIT NONE

  INTEGER :: nx,ny,nz          ! The number grid points in 3 directions

  INTEGER :: ibgn              ! i-index where evaluation begins.
  INTEGER :: iend              ! i-index where evaluation ends.
  INTEGER :: jbgn              ! j-index where evaluation begins.
  INTEGER :: jend              ! j-index where evaluation ends.

  REAL :: zp    (nx,ny,nz)     ! The physical height coordinate defined at
                               ! w-point of staggered grid.
  INTEGER :: soiltyp(nx,ny)    ! Soil type at each point.

  REAL :: wspd  (nx,ny)        ! Surface wind speed (m/s)

  REAL :: ptsfc(nx,ny)         ! Potential temperature at the ground level
                               ! (K)
  REAL :: pt1   (nx,ny)        ! Potential temperature at the 1st scalar
                               ! point above ground level, (K)

  REAL :: c_ptwtrneu(nx,ny)       ! Product of ustar and ptstar
  REAL :: c_ptwtr(nx,ny)       ! Product of ustar and ptstar
!
!-----------------------------------------------------------------------
!
!  Misc. local variables:
!
!-----------------------------------------------------------------------
!
  INTEGER :: i, j

  REAL :: z1                   ! The height of 1st scalar point above
                               ! ground level (m)
  REAL :: zo                   ! Defined as surface momentum roughness
                               ! length
  REAL :: zt                   ! Defined as surface heat transfer
                               ! roughness length
  REAL :: dzo                  ! z1-zo
  REAL :: dzt                  ! z1-zt
  REAL :: z1dzo                ! z1/zo
  REAL :: z1dzt                ! z1/zt
  REAL :: xcdn                 ! Cdn (sea)
  REAL :: xcdh                 ! Hot-wired value of Cdh (sea)

  REAL :: bulkri               ! Bulk Richardson number
  REAL :: stabp                ! Monin-Obukhov STABility Parameter
                               ! (zeta)

  REAL :: y1,y0,psih           ! Intermediate parameters needed
  REAL :: qb3pb2
  REAL :: c7,c8
  REAL :: sb,qb,pb,thetab,tb   ! during  computations
  REAL :: a,b,c,d
  REAL :: tempan,sqrtqb
!
!-----------------------------------------------------------------------
!
!  Include files:
!
!-----------------------------------------------------------------------
!
  INCLUDE 'phycst.inc'
  INCLUDE 'globcst.inc'
  INCLUDE 'sfcphycst.inc'
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
  xcdh = 1.1E-3

  DO j=jbgn,jend
    DO i=ibgn,iend
      IF ( soiltyp(i,j) == 13) THEN

        xcdn = (0.4+0.07*wspd(i,j)) * 1.e-3
        z1  = 0.5*(zp(i,j,3)-zp(i,j,2))
        z1  = MIN(z1,z1limit)
!
! calculate zo and zt
!
        zo =0.032*xcdn*wspd(i,j)*wspd(i,j)/9.8
        zo = MAX(zo,zolimit)
        zt = z1 * EXP( -kv*SQRT(xcdn)/(prantl0w*xcdh) )

        dzo = z1-zo
        z1dzo = z1/zo
        dzt = z1-zt
        z1dzt = z1/zt

        bulkri = (g/ptsfc(i,j))*(pt1(i,j)-ptsfc(i,j))*dzo*dzo/          &
                 (dzt*wspd(i,j)*wspd(i,j))

        IF (bulkri <= 0.0) THEN
!
!-----------------------------------------------------------------------
!
!  Unstable case: A modified formulation, which is similar to
!  equations (28)-(34) in Byun (1990).
!
!-----------------------------------------------------------------------
!
          bulkri = MAX (bulkri,-10.0)

          sb =bulkri/prantl0w

          qb=oned9*(c1w+c2w*sb*sb)
          pb=oned54*(c3w+c4w*sb*sb)

          qb3pb2=qb**3-pb*pb
          c7 = (z1*dzt/(dzo*dzo))*(ALOG(z1dzo)*ALOG(z1dzo)/ALOG(z1dzt))

          IF( qb3pb2 >= 0.0 ) THEN

            sqrtqb = SQRT(qb)
            tempan = MAX( -1.0, MIN( 1.0, pb/(sqrtqb**3) ) )

            thetab=ACOS(tempan)
            stabp =c7*(-2.0*SQRT(qb)*COS(thetab/3.0)+c5w)

          ELSE

            tb    =(SQRT(-qb3pb2)+ABS(pb))**oned3
            stabp =c7*(-(tb+qb/tb)+c5w)

          END IF
!
!-----------------------------------------------------------------------
!
!  According to a modified equation, which is similar to equation (14)
!  in Byun (1990).
!
!-----------------------------------------------------------------------
!
          c8=gammamw*stabp
          y1=SQRT(1.0 - c8)
          y0=SQRT(1.0 - c8/z1dzt)

          psih=2.0*ALOG((y1+1.0)/(y0+1.0))
!
!-----------------------------------------------------------------------
!
!  Compute ptstar via equation (11) in Byun (1981).
!
!-----------------------------------------------------------------------
!
          c_ptwtr(i,j)=kv / (prantl0w*(ALOG(z1dzt)-psih))

        ELSE
!
!-----------------------------------------------------------------------
!
!  Stable case: With the modified formulation in Louis et al (1981).
!
!-----------------------------------------------------------------------
!
          a=kv*kv/(prantl0w*LOG(z1dzo)*LOG(z1dzt))
          b=5.0
          d=5.0
          c=SQRT(1.0+d*bulkri)

          c_ptwtr(i,j) = SQRT(1.0+2.0*b*bulkri/c)
          c_ptwtr(i,j) = a*c_ptwtr(i,j)/(prantl0l*(1.0+3.0*b*bulkri*c))

        END IF
      END IF

    END DO
  END DO

  RETURN
END SUBROUTINE cptcwtr

!
!
!##################################################################
!##################################################################
!######                                                      ######
!######                SUBROUTINE CUNEUWTR                   ######
!######                                                      ######
!######                     Developed by                     ######
!######     Center for Analysis and Prediction of Storms     ######
!######                University of Oklahoma                ######
!######                                                      ######
!##################################################################
!##################################################################
!


SUBROUTINE cuneuwtr(nx,ny,nz,ibgn,iend,jbgn,jend,soiltyp,               & 1
           wspd, c_uneu)
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  Compute c_uneu (friction velocity) at the neutral state. The quantity
!  c_uneu is used by the subroutine SFCFLX to obtain surface fluxes.
!
!-----------------------------------------------------------------------
!
!  AUTHORS: V. Wong and X. Song
!  8/04/1994
!
!-----------------------------------------------------------------------
!
!  INPUT:
!
!    nx       Number of grid points in the x-direction (east/west)
!    ny       Number of grid points in the y-direction (north/south)
!    nz       Number of grid points in the vertical
!
!    ibgn     i-index where evaluation begins.
!    iend     i-index where evaluation ends.
!    jbgn     j-index where evaluation begins.
!    jend     j-index where evaluation ends.
!
!    soiltyp  Soil type at each point
!    wspd     Surface wind speed (m/s), defined as
!             sqrt(usuf*usuf + vsuf*vsuf + wsuf*wsuf)
!
!
!  OUTPUT:
!
!    c_uneu   Friction velocity at sea case.
!
!-----------------------------------------------------------------------
!

!
!-----------------------------------------------------------------------
!
!  Variable Declarations.
!
!-----------------------------------------------------------------------
!
  IMPLICIT NONE

  INTEGER :: nx,ny,nz          ! The number grid points in 3 directions

  INTEGER :: ibgn              ! i-index where evaluation begins.
  INTEGER :: iend              ! i-index where evaluation ends.
  INTEGER :: jbgn              ! j-index where evaluation begins.
  INTEGER :: jend              ! j-index where evaluation ends.

  INTEGER :: soiltyp(nx,ny)    ! Soil type at each point.
  REAL :: wspd  (nx,ny)        ! Surface wind speed (m/s)


  REAL :: c_uneu   (nx,ny)     ! Frictional velocity (m/s)
!
!-----------------------------------------------------------------------
!
!  Misc. local variables:
!
!-----------------------------------------------------------------------
!
  INTEGER :: i, j
!
!-----------------------------------------------------------------------
!
!  Include files:
!
!-----------------------------------------------------------------------
!
  INCLUDE 'phycst.inc'
  INCLUDE 'globcst.inc'
  INCLUDE 'sfcphycst.inc'
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!
  DO j=jbgn,jend
    DO i=ibgn,iend

      IF (soiltyp(i,j) == 13) THEN
        c_uneu(i,j) = SQRT ((0.4+0.079*wspd(i,j)) * 1.e-3)
      END IF

    END DO
  END DO

  RETURN
END SUBROUTINE cuneuwtr
!
!
!##################################################################
!##################################################################
!######                                                      ######
!######               SUBROUTINE CPTNEUWTR                   ######
!######                                                      ######
!######                     Developed by                     ######
!######     Center for Analysis and Prediction of Storms     ######
!######                University of Oklahoma                ######
!######                                                      ######
!##################################################################
!##################################################################
!


SUBROUTINE cptneuwtr(nx,ny,nz,ibgn,iend,jbgn,jend,                      & 1
           soiltyp,wspd,c_uneu,c_ptwtrneu)
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!-----------------------------------------------------------------------
!
!  AUTHORS: V. Wong and X. Song
!  8/04/1994
!
!  MODIFICATION:
!
!  03/17/1997 (Yuhe Liu and Vince Wong)
!  Fixed a bug. Variable c_ptwtrneu was calculated in its inverse.
!
!-----------------------------------------------------------------------
!
!  INPUT:
!
!    nx       Number of grid points in the x-direction (east/west)
!    ny       Number of grid points in the y-direction (north/south)
!    nz       Number of grid points in the vertical
!
!    ibgn     i-index where evaluation begins
!    iend     i-index where evaluation ends
!    jbgn     j-index where evaluation begins
!    jend     j-index where evaluation ends
!
!    wspd     Surface wind speed (m/s), defined as
!             sqrt(usuf*usuf + vsuf*vsuf + wsuf*wsuf)
!
!  OUTPUT:
!
!    c_ptwtrneu   Temperature scale (K), defined by
!             surface heat flux / friction velocity
!
!-----------------------------------------------------------------------
!

!
!-----------------------------------------------------------------------
!
!  Variable Declarations.
!
!-----------------------------------------------------------------------
!
  IMPLICIT NONE

  INTEGER :: nx,ny,nz          ! The number grid points in 3 directions

  INTEGER :: ibgn              ! i-index where evaluation begins.
  INTEGER :: iend              ! i-index where evaluation ends.
  INTEGER :: jbgn              ! j-index where evaluation begins.
  INTEGER :: jend              ! j-index where evaluation ends.

  INTEGER :: soiltyp (nx,ny)
  REAL :: wspd  (nx,ny)        ! Surface wind speed (m/s)

  REAL :: c_uneu (nx,ny)       !
  REAL :: c_ptwtrneu(nx,ny)    ! Temperature scale (K)
!
!-----------------------------------------------------------------------
!
!  Misc. local variables:
!
!-----------------------------------------------------------------------
!
  INTEGER :: i, j
  REAL :: xcdh
!
!-----------------------------------------------------------------------
!
!  Include files:
!
!-----------------------------------------------------------------------
!
  INCLUDE 'phycst.inc'
  INCLUDE 'globcst.inc'
  INCLUDE 'sfcphycst.inc'
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
  xcdh = 1.14E-3

  DO i=ibgn,iend
    DO j=jbgn,jend

      IF (soiltyp(i,j) == 13) THEN
        c_ptwtrneu(i,j) = xcdh/c_uneu(i,j)
      END IF

    END DO
  END DO

  RETURN
END SUBROUTINE cptneuwtr
!
!##################################################################
!##################################################################
!######                                                      ######
!######                 SUBROUTINE PBLDEPTH                  ######
!######                                                      ######
!######                     Developed by                     ######
!######     Center for Analysis and Prediction of Storms     ######
!######                University of Oklahoma                ######
!######                                                      ######
!##################################################################
!##################################################################
!


SUBROUTINE pbldepth(nx,ny,nz, u,v,w,ptprt,qv,ptbar,zp,j3, ptsflx,       & 1,2
           pbldpth)
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  Calcualte time-dependent PBL depth.
!
!-----------------------------------------------------------------------
!
!  AUTHORS: V. Wong and X. Song
!  8/27/1994
!
!  Reference for option pbldopt=3 (not included yet):
!
!  1)" Simple Model of the Daytime Boundary Layer Height" by S-E.
!      Gryning and E. Batchvarova in 9th Symposium on Turbulence and
!      Diffusion, 379-382 (1990).
!  2) "A Rate Equation for the Nocturnal Boundary-Layer Height" by
!      F.T.M. Nieuwatadt and H. Tennekes in J. of Atmospheric Sciences
!      July 1981, pp. 1418-1428.
!
!  MODIFICATIONS:
!
!  3/6/1996 (M. Xue, V. Wong and Y. Liu)
!  Added option 2 for diagnostic calculatin of PBL depth.
!
!  06/06/96 (Ming Xue and Jinxing Zong)
!  PBL depth is now determined to be at the level where environmental
!  virtual potential temperate equals that at the first grid level.
!  Subroutine PBLDEPTH in sfcphy3d.f was modified.
!
!
!-----------------------------------------------------------------------
!
!  INPUT:
!
!    nx       Number of grid points in the x-direction (east/west)
!    ny       Number of grid points in the y-direction (north/south)
!    nz       Number of grid points in the vertical
!
!
!    u        Total u-velocity (m/s)
!    v        Total v-velocity (m/s)
!    w        Total w-velocity (m/s)
!    ptprt    Perturbation potential temperature (K)
!    qv       Specific humidity (kg/kg)
!
!    ptbar    Base state potential temperature (K)
!
!    zp       The physical height coordinate defined at w-point of
!             staggered grid.
!    j3       Coordinate transformation Jacobian d(zp)/d(z)
!
!  OUTPUT:
!
!    pbldpth  Planetary boundary layer depth (m)
!
!  TEMPORARY:
!
!    ptv0     Work array for virtual pot. temp. at k=2 for scalar
!    ptv      Work array for virtual pot. temp.
!    parea    Work array for positive area
!    narea    Work array for negative area
!
!-----------------------------------------------------------------------
!
!  Variable Declarations.
!
!-----------------------------------------------------------------------
!
  IMPLICIT NONE

  INTEGER :: nx,ny,nz          ! The number grid points in 3 directions
  INTEGER :: nt                ! Number of time levels
  PARAMETER( nt=3 )

  REAL :: u      (nx,ny,nz) ! Total u-velocity (m/s)
  REAL :: v      (nx,ny,nz) ! Total v-velocity (m/s)
  REAL :: w      (nx,ny,nz) ! Total w-velocity (m/s)
  REAL :: ptprt  (nx,ny,nz) ! Perturbation potential temperature (K)
  REAL :: qv     (nx,ny,nz) ! Specific humidity (kg/kg)

  REAL :: ptbar  (nx,ny,nz) ! Base state potential temperature (K)

  REAL :: zp     (nx,ny,nz) ! The physical height coordinate defined at
                            ! w-point of staggered grid.
  REAL :: j3     (nx,ny,nz) ! Coordinate transformation Jacobian  d(zp)/d(z)

  REAL :: ptsflx (nx,ny)    ! Surface heat flux (K*kg/(s*m**2))

  REAL :: pbldpth(nx,ny,nt) ! Planetary boundary layer depth (m)

  REAL :: ptv0   (nx,ny)    ! Temporary work array
  REAL :: ptv    (nx,ny)    ! Temporary work array
  REAL :: parea  (nx,ny)    ! Temporary work array
  REAL :: narea  (nx,ny)    ! Temporary work array
!
!-----------------------------------------------------------------------
!
!  Misc. local variables:
!
!-----------------------------------------------------------------------
!
  INTEGER :: i, j, k
  REAL :: zptk0,ptvk0,zptk1,ptvk1,eps, amax,amin
  INTEGER :: pbldepth_set
!
!-----------------------------------------------------------------------
!
!  Include files:
!
!-----------------------------------------------------------------------
!
  INCLUDE 'phycst.inc'
  INCLUDE 'globcst.inc'
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!
  IF ( pbldopt == 1 ) THEN

    DO i=1,nx-1
      DO j=1,ny-1
        pbldpth(i,j,2)=pbldpth0
      END DO
    END DO

  ELSE IF ( pbldopt == 2 ) THEN
!
!-----------------------------------------------------------------------
!
!  Diagnose the PBL depth based on the virtual potential temperature
!  profile. The PBL top is assumed to be at the level where the
!  environmental virtual potential temperature exceeds that at the
!  surface (first level above ground). If the atmosphere is stable
!  right above ground, the PBL depth is set to the thickness of the
!  layer below the first scalar point above ground.
!
!-----------------------------------------------------------------------
!
    DO j=1,ny-1
      DO i=1,nx-1
        ptv0(i,j)=(ptbar(i,j,2)+ptprt(i,j,2))*(1.+0.608*qv(i,j,2))
        parea(i,j) = 0.0
        narea(i,j) = 0.0
        pbldpth(i,j,2)=0.0  ! Set the initial depth to zero.
                            ! pbldpth also acts as a flag.
      END DO
    END DO

    eps = 1.0E-6
!
    DO k=3,nz-1

      DO j=1,ny-1
        DO i=1,nx-1

          IF( ABS(pbldpth(i,j,2)) < eps ) THEN
                              ! Check if pbldpth has been set.

            ptvk1=(ptbar(i,j,k  ) + ptprt(i,j,k  ))                     &
                       * (1. + 0.608*qv(i,j,k  ))
            ptvk0=(ptbar(i,j,k-1) + ptprt(i,j,k-1))                     &
                       * (1. + 0.608*qv(i,j,k-1))

            zptk0 = 0.5*(zp(i,j,k)+zp(i,j,k-1))
            zptk1 = 0.5*(zp(i,j,k+1)+zp(i,j,k))

            IF( ptvk1 > ptv0(i,j) ) THEN

              pbldpth(i,j,2) = zptk0 + (zptk1-zptk0)*                   &
                  (ptv0(i,j)-ptvk0)/(ptvk1-ptvk0) - zp(i,j,2)

            END IF

          END IF

        END DO
      END DO

    END DO

!-----------------------------------------------------------------------
!
!  When no stable layer is found above, the PBL top is at the model top.
!  The minimum PBL top height is at the first scalar point above ground.
!
!-----------------------------------------------------------------------

    pbldepth_set = 1

    DO j=1,ny-1
      DO i=1,nx-1
        IF(pbldpth(i,j,2) == 0.0) THEN
          pbldpth(i,j,2)=zp(i,j,nz-1)-zp(i,j,2)

          WRITE(6,'(/1x,a,a,2i4, 2(/1x,a,a))')                          &
              'Warning: The PBL is found to extend ',                   &
              'the entire model domain at i,j=',i,j,                    &
              'The atmosphere is either neutral or ',                   &
              'unstable throughout the domain depth.',                  &
              'In this case, PBL parameterization, if used, will operate', &
              ' on the entire domain depth.'
          pbldepth_set = 0
        END IF
        pbldpth(i,j,2) =                                                &
               MAX(pbldpth(i,j,2), 0.5*(zp(i,j,3)-zp(i,j,2)))
      END DO
    END DO

    IF(pbldepth_set == 0)THEN
      CALL a3dmax0(ptv0,1,nx,1,nx-1,1,ny,1,ny-1,1,1,1,1,amax,amin)
      WRITE(6,'(1x,2(a,e13.6))')'Max sfc virtual temperature =',amax
    END IF

  ELSE IF ( pbldopt == 3 ) THEN

    WRITE(6,'(/5x,a,i3,a,2(/5x,a))')                                    &
        'PBL depth calculation option=',pbldopt,' not avaiable.',       &
        'please reset parameter pbldopt in input file. ',               &
        'Program stopped in PBLDEPTH.'

    CALL arpsstop('arpsstop called from PBLDEPTH pbldopt incorrect ',1)

  END IF


  RETURN
END SUBROUTINE pbldepth