pkg/atm_phys/my25_turb_mod.F90

b2ea1d2979 Jean*0001   MODULE MY25_TURB_MOD
                0002
                0003 !=======================================================================
                0004 !   MELLOR-YAMADA LEVEL 2.5 TURBULENCE CLOSURE SCHEME - GFDL VERSION   !
                0005 !=======================================================================
                0006
                0007 ! use Fms_Mod,         ONLY: FILE_EXIST, OPEN_NAMELIST_FILE, ERROR_MESG, FATAL,&
                0008 !                            read_data, write_data, CLOSE_FILE,&
                0009 !                            check_nml_error,mpp_pe, mpp_root_pe, &
                0010 !                            write_version_number, stdlog, open_file
                0011 ! use Tridiagonal_Mod, ONLY: TRI_INVERT, CLOSE_TRIDIAGONAL
                0012  use constants_mod,   only: grav, vonkarm
                0013  use monin_obukhov_mod, only : mo_diff
                0014
                0015 !---------------------------------------------------------------------
                0016  implicit none
                0017  private
                0018 !---------------------------------------------------------------------
                0019
                0020  public :: MY25_TURB, MY25_TURB_INIT, MY25_TURB_END, TKE_SURF
                0021
                0022 !---------------------------------------------------------------------
                0023 ! --- GLOBAL STORAGE
                0024 !---------------------------------------------------------------------
                0025
                0026   real, public, allocatable, dimension(:,:,:) :: TKE
                0027
                0028 !---------------------------------------------------------------------
                0029
                0030  character(len=128) :: version = '$Id: my25_turb_mod.F90,v 1.1 2013/05/08 22:14:14 jmc Exp $'
                0031  character(len=128) :: tag = '$Name:  $'
                0032
                0033  logical :: do_init = .true.
                0034  logical :: init_tke
                0035  integer :: num_total_pts, pts_done
                0036
                0037 !---------------------------------------------------------------------
                0038 ! --- CONSTANTS
                0039 !---------------------------------------------------------------------
                0040
                0041  real :: aa1,   aa2,   bb1,  bb2,  ccc
                0042  real :: ckm1,  ckm2,  ckm3, ckm4, ckm5, ckm6, ckm7, ckm8
                0043  real :: ckh1,  ckh2,  ckh3, ckh4
                0044  real :: cvfq1, cvfq2, bcq
                0045
                0046  real, parameter :: aa1_old =  0.78
                0047  real, parameter :: aa2_old =  0.79
                0048  real, parameter :: bb1_old = 15.0
                0049  real, parameter :: bb2_old =  8.0
                0050  real, parameter :: ccc_old =  0.056
                0051  real, parameter :: aa1_new =  0.92
                0052  real, parameter :: aa2_new =  0.74
                0053  real, parameter :: bb1_new = 16.0
                0054  real, parameter :: bb2_new = 10.0
                0055  real, parameter :: ccc_new =  0.08
                0056  real, parameter :: cc1     =  0.27
                0057  real, parameter :: t00     =  2.7248e2
                0058  real, parameter :: small   =  1.0e-10
                0059
                0060 !---------------------------------------------------------------------
                0061 ! --- NAMELIST
                0062 !---------------------------------------------------------------------
                0063
                0064  real    :: TKEmax       =  5.0
                0065  real    :: TKEmin       =  0.0
                0066  real    :: el0max       =  1.0e6
                0067  real    :: el0min       =  0.0
                0068  real    :: alpha_land   =  0.10
                0069  real    :: alpha_sea    =  0.10
                0070  real    :: akmax        =  1.0e4
                0071  real    :: akmin_land   =  5.0
                0072  real    :: akmin_sea    =  0.0
                0073  integer :: nk_lim       =  2
                0074  integer :: init_iters   =  20
                0075  logical :: do_thv_stab  = .true.
                0076  logical :: use_old_cons = .false.
                0077  real    :: kcrit        =  0.01
                0078
                0079   NAMELIST / my25_turb_nml /                           &
                0080          TKEmax,   TKEmin,   init_iters,               &
                0081          akmax,    akmin_land, akmin_sea, nk_lim,      &
                0082          el0max,   el0min,  alpha_land,  alpha_sea,    &
                0083          do_thv_stab, use_old_cons,                    &
                0084          kcrit
                0085
                0086 !---------------------------------------------------------------------
                0087
                0088  contains
                0089
                0090 !#######################################################################
                0091
                0092  SUBROUTINE MY25_TURB( delt, fracland, phalf, pfull, theta, &
                0093                        um,   vm,       zhalf, zfull, z0,    &
                0094                        TKE,  el0,      el,    akm,   akh,   &
                0095                        mask, kbot,     ustar, bstar, h    )
                0096
                0097 !=======================================================================
                0098 !---------------------------------------------------------------------
                0099 ! Arguments (Intent in)
                0100 !       delt     -  Time step in seconds
                0101 !       fracland -  Fractional amount of land beneath a grid box
                0102 !       phalf    -  Pressure at half levels
                0103 !       pfull    -  Pressure at full levels
                0104 !       theta    -  Potential temperature
                0105 !       um, vm   -  Wind components
                0106 !       zhalf    -  Height at half levels
                0107 !       zfull    -  Height at full levels
                0108 !       z0       -  Roughness length
                0109 !       mask     -  OPTIONAL; floating point mask (0. or 1.) designating
                0110 !                   where data is present
                0111 !       kbot     -  OPTIONAL;lowest model level index (integer);
                0112 !                    at levels > kbot, mask = 0.
                0113 !       ustar    -  OPTIONAL:friction velocity (m/sec)
                0114 !       bstar    -  OPTIONAL:buoyancy scale (m/sec**2)
                0115 !---------------------------------------------------------------------
                0116   real,    intent(in)                   :: delt
                0117   real,    intent(in), dimension(:,:)   :: fracland, z0
                0118   real,    intent(in), dimension(:,:,:) :: phalf, pfull, zhalf, zfull
                0119   real,    intent(in), dimension(:,:,:) :: um, vm, theta
                0120
                0121   integer, intent(in), OPTIONAL, dimension(:,:)   :: kbot
                0122   real,    intent(in), OPTIONAL, dimension(:,:,:) :: mask
                0123   real,    intent(in), OPTIONAL, dimension(:,:)   :: ustar, bstar
                0124
                0125 !---------------------------------------------------------------------
                0126 ! Arguments (Intent in/out)
                0127 !       TKE  -  turbulent kinetic energy
                0128 !---------------------------------------------------------------------
                0129   real, intent(inout), dimension(:,:,:) :: TKE
                0130
                0131 !---------------------------------------------------------------------
                0132 ! Arguments (Intent out)
                0133 !       el0  -  characteristic length scale
                0134 !       el   -  master length scale
                0135 !       akm  -  mixing coefficient for momentum
                0136 !       akh  -  mixing coefficient for heat and moisture
                0137 !         h  -  OPTIONAL, diagnosed depth of planetary boundary
                0138 !                         layer (m)
                0139 !---------------------------------------------------------------------
                0140   real, intent(out), dimension(:,:)   :: el0
                0141   real, intent(out), dimension(:,:,:) :: akm, akh, el
                0142   real, intent(out), OPTIONAL, dimension(:,:) :: h
                0143
                0144 !---------------------------------------------------------------------
                0145 !  (Intent local)
                0146 !---------------------------------------------------------------------
                0147   integer :: ix, jx, kx, i, j, k
                0148   integer :: kxp, kxm, klim, it, itermax
                0149   real    :: cvfqdt, dvfqdt
                0150
                0151   real, dimension(SIZE(um,1),SIZE(um,2)) :: zsfc, x1, x2, akmin
                0152
                0153   real, dimension(SIZE(um,1),SIZE(um,2),SIZE(um,3)-1) ::     &
                0154         dsdzh, shear, buoync, qm2,  qm3, qm4, el2,           &
                0155         aaa,   bbb,   ccc,    ddd,                           &
                0156         xxm1,  xxm2,  xxm3,   xxm4, xxm5
                0157
                0158   real, dimension(SIZE(um,1),SIZE(um,2),SIZE(um,3)) ::       &
                0159         dsdz, qm,  xx1, xx2
                0160
                0161 !====================================================================
                0162
                0163 ! --- Check to see if MY25_TURB has been initialized
                0164 !del  if( do_init ) CALL ERROR_MESG( ' MY25_TURB',     &
                0165 !del                                 ' MY25_TURB_INIT has not been called',&
                0166 !del                                   FATAL )
                0167
                0168 ! --- Set dimensions etc
                0169   ix  = SIZE( um, 1 )
                0170   jx  = SIZE( um, 2 )
                0171   kx  = SIZE( um, 3 )
                0172   kxp = kx + 1
                0173   kxm = kx - 1
                0174
                0175 !====================================================================
                0176 ! --- SURFACE HEIGHT
                0177 !====================================================================
                0178
                0179   if( PRESENT( kbot ) ) then
                0180      do j = 1,jx
                0181      do i = 1,ix
                0182         k = kbot(i,j) + 1
                0183        zsfc(i,j) = zhalf(i,j,k)
                0184      end do
                0185      end do
                0186   else
                0187        zsfc(:,:) = zhalf(:,:,kxp)
                0188   endif
                0189
                0190 !====================================================================
                0191 ! --- D( )/DZ OPERATORS: AT FULL LEVELS & AT HALF LEVELS
                0192 !====================================================================
                0193
                0194    dsdz(:,:,1:kx)  = 1.0 / ( zhalf(:,:,2:kxp) - zhalf(:,:,1:kx) )
                0195   dsdzh(:,:,1:kxm) = 1.0 / ( zfull(:,:,2:kx)  - zfull(:,:,1:kxm) )
                0196
                0197 !====================================================================
                0198 ! --- WIND SHEAR
                0199 !====================================================================
                0200
                0201   xxm1(:,:,1:kxm) = dsdzh(:,:,1:kxm)*( um(:,:,2:kx) - um(:,:,1:kxm) )
                0202   xxm2(:,:,1:kxm) = dsdzh(:,:,1:kxm)*( vm(:,:,2:kx) - vm(:,:,1:kxm) )
                0203
                0204   shear = xxm1 * xxm1 + xxm2 * xxm2
                0205
                0206 !====================================================================
                0207 ! --- BUOYANCY
                0208 !====================================================================
                0209
                0210   xxm1(:,:,1:kxm) = theta(:,:,2:kx) - theta(:,:,1:kxm)
                0211
                0212   if( do_thv_stab ) then
                0213      xxm2(:,:,1:kxm) = 0.5*( theta(:,:,2:kx) + theta(:,:,1:kxm) )
                0214   else
                0215      xxm2(:,:,1:kxm) = t00
                0216   end if
                0217
                0218   buoync = grav * dsdzh * xxm1 / xxm2
                0219
                0220 !====================================================================
                0221 ! --- MASK OUT UNDERGROUND VALUES FOR ETA COORDINATE
                0222 !====================================================================
                0223
                0224   if( PRESENT( mask ) ) then
                0225     where (mask(:,:,2:kx) < 0.1)   ! assume mask(k=1) always eq 1 ?
                0226          TKE(:,:,3:kxp) = 0.0
                0227         dsdz(:,:,2:kx ) = 0.0
                0228        dsdzh(:,:,1:kxm) = 0.0
                0229        shear(:,:,1:kxm) = 0.0
                0230       buoync(:,:,1:kxm) = 0.0
                0231     endwhere
                0232   endif
                0233
                0234 !====================================================================
                0235 ! --- SET ITERATION LOOP IF INITALIZING TKE
                0236 !====================================================================
                0237
                0238             itermax = 1
                0239      if (init_tke) then
                0240             itermax = init_iters
                0241             pts_done = pts_done + ix*jx
                0242             if (pts_done >= num_total_pts) init_tke = .false.
                0243      endif
                0244
                0245 ! $$$$$$$$$$$$$$$$$
                0246   do it = 1,itermax
                0247 ! $$$$$$$$$$$$$$$$$
                0248
                0249 !====================================================================
                0250 ! --- SOME TKE STUFF
                0251 !====================================================================
                0252
                0253   xx1(:,:,1:kx)  = 2.0 * TKE(:,:,2:kxp)
                0254   where (xx1(:,:,1:kx) > 0.0)
                0255     qm(:,:,1:kx)  = SQRT( xx1(:,:,1:kx) )
                0256   elsewhere
                0257     qm(:,:,1:kx)  = 0.0
                0258   endwhere
                0259
                0260   qm2(:,:,1:kxm)  = xx1(:,:,1:kxm)
                0261   qm3(:,:,1:kxm)  =  qm(:,:,1:kxm) * qm2(:,:,1:kxm)
                0262   qm4(:,:,1:kxm)  = qm2(:,:,1:kxm) * qm2(:,:,1:kxm)
                0263
                0264 !====================================================================
                0265 ! --- CHARACTERISTIC LENGTH SCALE
                0266 !====================================================================
                0267
                0268   xx1(:,:,1:kxm) = qm(:,:,1:kxm)*( pfull(:,:,2:kx) - pfull(:,:,1:kxm) )
                0269   do k = 1, kxm
                0270      xx2(:,:,k) = xx1(:,:,k)  * ( zhalf(:,:,k+1) - zsfc(:,:) )
                0271   end do
                0272
                0273   if( PRESENT( kbot ) ) then
                0274        xx1(:,:,kx) = 0.0
                0275        xx2(:,:,kx) = 0.0
                0276      do j = 1,jx
                0277      do i = 1,ix
                0278         k = kbot(i,j)
                0279        xx1(i,j,k)  =  qm(i,j,k) * ( phalf(i,j,k+1) - pfull(i,j,k) )
                0280        xx2(i,j,k)  = xx1(i,j,k) * z0(i,j)
                0281      end do
                0282      end do
                0283   else
                0284        xx1(:,:,kx) =  qm(:,:,kx) * ( phalf(:,:,kxp) - pfull(:,:,kx) )
                0285        xx2(:,:,kx) = xx1(:,:,kx) * z0(:,:)
                0286   endif
                0287
                0288   if (PRESENT(mask)) then
                0289     x1 = SUM( xx1, 3, mask=mask.gt.0.1 )
                0290     x2 = SUM( xx2, 3, mask=mask.gt.0.1 )
                0291   else
                0292     x1 = SUM( xx1, 3 )
                0293     x2 = SUM( xx2, 3 )
                0294   endif
                0295
                0296 !---- should never be equal to zero ----
                0297 !del  if (count(x1 <= 0.0) > 0) CALL ERROR_MESG( ' MY25_TURB',  &
                0298 !del                             'divid by zero, x1 <= 0.0', FATAL)
                0299   el0 = x2 / x1
                0300   el0 = el0 * (alpha_land*fracland + alpha_sea*(1.-fracland))
                0301
                0302   el0 = MIN( el0, el0max )
                0303   el0 = MAX( el0, el0min )
                0304
                0305 !====================================================================
                0306 ! --- MASTER LENGTH SCALE
                0307 !====================================================================
                0308
                0309   do k = 1, kxm
                0310      xx1(:,:,k)  = vonkarm * ( zhalf(:,:,k+1) - zsfc(:,:) )
                0311   end do
                0312
                0313   x1(:,:) = vonkarm * z0(:,:)
                0314
                0315   if( PRESENT( kbot ) ) then
                0316      do j = 1,jx
                0317      do i = 1,ix
                0318         do k = kbot(i,j), kx
                0319           xx1(i,j,k) = x1(i,j)
                0320         end do
                0321      end do
                0322      end do
                0323   else
                0324         xx1(:,:,kx) = x1(:,:)
                0325   endif
                0326
                0327   do k = 1,kx
                0328     el(:,:,k+1) = xx1(:,:,k) / ( 1.0 + xx1(:,:,k) / el0(:,:) )
                0329   end do
                0330     el(:,:,1)   = el0(:,:)
                0331
                0332   el2(:,:,1:kxm) = el(:,:,2:kx) * el(:,:,2:kx)
                0333
                0334 !====================================================================
                0335 ! --- MIXING COEFFICIENTS
                0336 !====================================================================
                0337
                0338   xxm3(:,:,1:kxm) = el2(:,:,1:kxm)*buoync(:,:,1:kxm)
                0339   xxm4(:,:,1:kxm) = el2(:,:,1:kxm)* shear(:,:,1:kxm)
                0340   xxm5(:,:,1:kxm) =  el(:,:,2:kx )*   qm3(:,:,1:kxm)
                0341
                0342 !-------------------------------------------------------------------
                0343 ! --- MOMENTUM
                0344 !-------------------------------------------------------------------
                0345
                0346   xxm1 = xxm5*( ckm1*qm2 + ckm2*xxm3 )
                0347   xxm2 = qm4 + ckm5*qm2*xxm4 + xxm3*( ckm6*xxm4 + ckm7*qm2 + ckm8*xxm3 )
                0348
                0349   xxm2 = MAX( xxm2, 0.2*qm4 )
                0350   xxm2 = MAX( xxm2, small  )
                0351
                0352   akm(:,:,1)    = 0.0
                0353   akm(:,:,2:kx) = xxm1(:,:,1:kxm) / xxm2(:,:,1:kxm)
                0354
                0355   akm = MAX( akm, 0.0 )
                0356
                0357 !-------------------------------------------------------------------
                0358 ! --- HEAT AND MOISTURE
                0359 !-------------------------------------------------------------------
                0360
                0361   xxm1(:,:,1:kxm) = ckh1*xxm5(:,:,1:kxm) - ckh2*xxm4(:,:,1:kxm)*akm(:,:,2:kx)
                0362   xxm2(:,:,1:kxm) = qm2(:,:,1:kxm) + ckh3*xxm3(:,:,1:kxm)
                0363
                0364   xxm1 = MAX( xxm1, ckh4*xxm5 )
                0365   xxm2 = MAX( xxm2, 0.4*qm2   )
                0366   xxm2 = MAX( xxm2, small     )
                0367
                0368   akh(:,:,1)    = 0.0
                0369   akh(:,:,2:kx) = xxm1(:,:,1:kxm) / xxm2(:,:,1:kxm)
                0370
                0371 !-------------------------------------------------------------------
                0372 ! --- BOUNDS
                0373 !-------------------------------------------------------------------
                0374
                0375 ! --- UPPER BOUND
                0376   akm = MIN( akm, akmax )
                0377   akh = MIN( akh, akmax )
                0378
                0379 ! --- LOWER BOUND
                0380 !  where( akm(:,:,1:klim) < small )  akm(:,:,1:klim) = 0.0
                0381 !  where( akh(:,:,1:klim) < small )  akh(:,:,1:klim) = 0.0
                0382
                0383 ! --- LOWER BOUND NEAR SURFACE
                0384
                0385   akmin = akmin_land*fracland + akmin_sea*(1.-fracland)
                0386
                0387   if( PRESENT( kbot ) ) then
                0388      do j = 1,jx
                0389      do i = 1,ix
                0390              klim = kbot(i,j) - nk_lim + 1
                0391      do  k = klim,kbot(i,j)
                0392         akm(i,j,k) = MAX( akm(i,j,k), akmin(i,j) )
                0393         akh(i,j,k) = MAX( akh(i,j,k), akmin(i,j) )
                0394      end do
                0395      end do
                0396      end do
                0397   else
                0398              klim = kx - nk_lim + 1
                0399      do  k = klim,kx
                0400         akm(:,:,k) = MAX( akm(:,:,k), akmin(:,:) )
                0401         akh(:,:,k) = MAX( akh(:,:,k), akmin(:,:) )
                0402      end do
                0403   endif
                0404
                0405 !-------------------------------------------------------------------
                0406 ! --- MASK OUT UNDERGROUND VALUES FOR ETA COORDINATE
                0407 !-------------------------------------------------------------------
                0408
                0409  if( PRESENT( mask ) ) then
                0410      akm(:,:,1:kx) = akm(:,:,1:kx) * mask(:,:,1:kx)
                0411      akh(:,:,1:kx) = akh(:,:,1:kx) * mask(:,:,1:kx)
                0412  endif
                0413
                0414 !====================================================================
                0415 ! --- PROGNOSTICATE TURBULENT KE
                0416 !====================================================================
                0417
                0418   cvfqdt = cvfq1 * delt
                0419   dvfqdt = cvfq2 * delt * 2.0
                0420
                0421 !-------------------------------------------------------------------
                0422 ! --- PART OF LINEARIZED ENERGY DISIIPATION TERM
                0423 !-------------------------------------------------------------------
                0424
                0425   xxm1(:,:,1:kxm) = dvfqdt * qm(:,:,1:kxm) / el(:,:,2:kx)
                0426
                0427 !-------------------------------------------------------------------
                0428 ! --- PART OF LINEARIZED VERTICAL DIFFUSION TERM
                0429 !-------------------------------------------------------------------
                0430
                0431   xx1(:,:,1:kx) = el(:,:,2:kxp) * qm(:,:,1:kx)
                0432
                0433   xx2(:,:,1)    = 0.5*  xx1(:,:,1)
                0434   xx2(:,:,2:kx) = 0.5*( xx1(:,:,2:kx) + xx1(:,:,1:kxm) )
                0435
                0436   xx1 = xx2 * dsdz
                0437
                0438 !-------------------------------------------------------------------
                0439 ! --- IMPLICIT TIME DIFFERENCING FOR VERTICAL DIFFUSION
                0440 ! --- AND ENERGY DISSIPATION TERM
                0441 !-------------------------------------------------------------------
                0442
                0443   aaa(:,:,1:kxm) = -cvfqdt * xx1(:,:,2:kx ) * dsdzh(:,:,1:kxm)
                0444   ccc(:,:,1:kxm) = -cvfqdt * xx1(:,:,1:kxm) * dsdzh(:,:,1:kxm)
                0445   bbb(:,:,1:kxm) =     1.0 - aaa(:,:,1:kxm) -   ccc(:,:,1:kxm)
                0446   bbb(:,:,1:kxm) =           bbb(:,:,1:kxm) +  xxm1(:,:,1:kxm)
                0447   ddd(:,:,1:kxm) =           TKE(:,:,2:kx )
                0448
                0449 ! correction for vertical diffusion of TKE surface boundary condition
                0450
                0451   if (present(kbot)) then
                0452      do j = 1,jx
                0453      do i = 1,ix
                0454           k = kbot(i,j)
                0455           ddd(:,:,k-1) = ddd(:,:,k-1) - aaa(:,:,k-1) * TKE(:,:,k+1)
                0456      enddo
                0457      enddo
                0458   else
                0459           ddd(:,:,kxm) = ddd(:,:,kxm) - aaa(:,:,kxm) * TKE(:,:,kxp)
                0460   endif
                0461
                0462 ! mask out terms below ground
                0463
                0464   if (present(mask)) then
                0465      where (mask(:,:,2:kx) < 0.1) ddd(:,:,1:kxm) = 0.0
                0466   endif
                0467
                0468
                0469 !del  CALL TRI_INVERT( xxm1, ddd, aaa, bbb, ccc )
                0470 !del  CALL CLOSE_TRIDIAGONAL
                0471
                0472 !-------------------------------------------------------------------
                0473 ! --- MASK OUT UNDERGROUND VALUES FOR ETA COORDINATE
                0474 !-------------------------------------------------------------------
                0475
                0476  if( PRESENT( mask ) ) then
                0477     where (mask(:,:,2:kx) < 0.1) xxm1(:,:,1:kxm) = TKE(:,:,2:kx)
                0478  endif
                0479
                0480 !-------------------------------------------------------------------
                0481 ! --- SHEAR AND BUOYANCY TERMS
                0482 !-------------------------------------------------------------------
                0483
                0484   xxm2(:,:,1:kxm) =  delt*( akm(:,:,2:kx)* shear(:,:,1:kxm)    &
                0485                           - akh(:,:,2:kx)*buoync(:,:,1:kxm) )
                0486
                0487 !-------------------------------------------------------------------
                0488 ! --- UPDATE TURBULENT KINETIC ENERGY
                0489 !-------------------------------------------------------------------
                0490
                0491   TKE(:,:,1)    = 0.0
                0492   TKE(:,:,2:kx) = xxm1(:,:,1:kxm) +  xxm2(:,:,1:kxm)
                0493
                0494 !====================================================================
                0495 ! --- BOUND TURBULENT KINETIC ENERGY
                0496 !====================================================================
                0497
                0498   TKE = MIN( TKE, TKEmax )
                0499   TKE = MAX( TKE, TKEmin )
                0500
                0501   if( PRESENT( mask ) ) then
                0502      where (mask(:,:,1:kx) < 0.1) TKE(:,:,2:kxp) = 0.0
                0503   endif
                0504
                0505 !====================================================================
                0506 ! --- COMPUTE PBL DEPTH IF DESIRED
                0507 !====================================================================
                0508
                0509   if (present(h)) then
                0510
                0511       if (.not.present(ustar).or..not.present(bstar)) then
                0512 !del          CALL ERROR_MESG( ' MY25_TURB',     &
                0513 !del              'cannot request pbl depth diagnostic if ustar'// &
                0514 !del              ' and bstar are not also supplied', FATAL )
                0515       end if
                0516       if (present(kbot)) then
                0517           call k_pbl_depth(ustar,bstar,akm,akh,zsfc,zfull,zhalf,&
                0518                            h,kbot=kbot)
                0519       else
                0520           call k_pbl_depth(ustar,bstar,akm,akh,zsfc,zfull,zhalf,h)
                0521       end if
                0522
                0523   end if
                0524 !====================================================================
                0525
                0526 ! $$$$$$$$$$$$$$$$$
                0527   end do
                0528 ! $$$$$$$$$$$$$$$$$
                0529
                0530 !====================================================================
                0531   end SUBROUTINE MY25_TURB
                0532
                0533 !#######################################################################
                0534
                0535   SUBROUTINE MY25_TURB_INIT( ix, jx, kx )
                0536
                0537 !=======================================================================
                0538 ! ***** INITIALIZE MELLOR-YAMADA
                0539 !=======================================================================
                0540 !---------------------------------------------------------------------
                0541 ! Arguments (Intent in)
                0542 !     ix, jx  - Horizontal dimensions for global storage arrays
                0543 !     kx      - Number of vertical levels in model
                0544 !---------------------------------------------------------------------
                0545  integer, intent(in) :: ix, jx, kx
                0546 !---------------------------------------------------------------------
                0547 !  (Intent local)
                0548 !---------------------------------------------------------------------
                0549  integer             :: unit, io, ierr
                0550  real                :: actp, facm
                0551  real, dimension(15) :: au,   tem
                0552
                0553 !=====================================================================
                0554
                0555 !---------------------------------------------------------------------
                0556 ! --- Read namelist
                0557 !---------------------------------------------------------------------
                0558
                0559 !del  if( FILE_EXIST( 'input.nml' ) ) then
                0560 ! -------------------------------------
                0561 !del   unit = OPEN_NAMELIST_FILE ()
                0562 !del   ierr = 1
                0563 !del   do while( ierr .ne. 0 )
                0564 !del   READ ( unit,  nml = my25_turb_nml, iostat = io, end = 10 )
                0565 !del   ierr = check_nml_error (io, 'my25_turb_nml')
                0566 !del   end do
                0567 10 continue
                0568 !del   CALL CLOSE_FILE( unit )
                0569 ! -------------------------------------
                0570 !del  end if
                0571
                0572 !---------------------------------------------------------------------
                0573 ! --- Output version
                0574 !---------------------------------------------------------------------
                0575
                0576
                0577 !del      if ( mpp_pe() == mpp_root_pe() )write (stdlog(), nml=my25_turb_nml)
                0578 !del      call write_version_number(version, tag)
                0579
                0580 !---------------------------------------------------------------------
                0581 ! --- Initialize constants
                0582 !---------------------------------------------------------------------
                0583
                0584   if( use_old_cons ) then
                0585       aa1 = aa1_old
                0586       aa2 = aa2_old
                0587       bb1 = bb1_old
                0588       bb2 = bb2_old
                0589       ccc = ccc_old
                0590   else
                0591       aa1 = aa1_new
                0592       aa2 = aa2_new
                0593       bb1 = bb1_new
                0594       bb2 = bb2_new
                0595       ccc = ccc_new
                0596   end if
                0597
                0598      ckm1 = ( 1.0 - 3.0*ccc )*aa1
                0599      ckm3 =  3.0 * aa1*aa2*    ( bb2 - 3.0*aa2 )
                0600      ckm4 =  9.0 * aa1*aa2*ccc*( bb2 + 4.0*aa1 )
                0601      ckm5 =  6.0 * aa1*aa1
                0602      ckm6 = 18.0 * aa1*aa1*aa2*( bb2 - 3.0*aa2 )
                0603      ckm7 =  3.0 * aa2*        ( bb2 + 7.0*aa1 )
                0604      ckm8 = 27.0 * aa1*aa2*aa2*( bb2 + 4.0*aa1 )
                0605      ckm2 =  ckm3 - ckm4
                0606      ckh1 =  aa2
                0607      ckh2 =  6.0 * aa1*aa2
                0608      ckh3 =  3.0 * aa2*( bb2 + 4.0*aa1 )
                0609      ckh4 =  2.0e-6 * aa2
                0610     cvfq1 = 5.0 * cc1 / 3.0
                0611     cvfq2 = 1.0 / bb1
                0612       bcq = 0.5 * ( bb1**(2.0/3.0) )
                0613
                0614 !---------------------------------------------------------------------
                0615 ! --- Allocate storage for TKE
                0616 !---------------------------------------------------------------------
                0617
                0618   if( ALLOCATED( TKE ) ) DEALLOCATE( TKE )
                0619                            ALLOCATE( TKE(ix,jx,kx+1) )
                0620
                0621 !-------------------------------------------------------------------
                0622   do_init = .false.
                0623 !---------------------------------------------------------------------
                0624 ! --- Input TKE
                0625 !---------------------------------------------------------------------
                0626
                0627 !del  if( FILE_EXIST( 'INPUT/my25_turb.res' ) ) then
                0628
                0629 !del      unit = OPEN_FILE ( file = 'INPUT/my25_turb.res', &
                0630 !del                         form = 'native', action = 'read' )
                0631 !del      call read_data ( unit, TKE )
                0632 !del      CALL CLOSE_FILE( unit )
                0633
                0634       init_tke = .false.
                0635
                0636 !del  else
                0637
                0638       TKE  = TKEmin
                0639
                0640       init_tke      = .true.
                0641       num_total_pts = ix*jx
                0642       pts_done      = 0
                0643
                0644 !del  endif
                0645
                0646
                0647 !=====================================================================
                0648   end SUBROUTINE MY25_TURB_INIT
                0649
                0650 !#######################################################################
                0651
                0652   SUBROUTINE MY25_TURB_END
                0653 !=======================================================================
                0654  integer :: unit
                0655 !=======================================================================
                0656
                0657 !del      unit = OPEN_FILE ( file = 'RESTART/my25_turb.res', &
                0658 !del                         form = 'native', action = 'write' )
                0659 !del      call write_data ( unit, TKE )
                0660 !del      CALL CLOSE_FILE ( unit )
                0661
                0662 !=====================================================================
                0663
                0664   end SUBROUTINE MY25_TURB_END
                0665
                0666 !#######################################################################
                0667
                0668   SUBROUTINE K_PBL_DEPTH(ustar,bstar,akm,akh,zsfc,zfull,zhalf,h,kbot)
                0669 !=======================================================================
                0670
                0671   real,    intent(in), dimension(:,:)   :: ustar, bstar, zsfc
                0672   real,    intent(in), dimension(:,:,:) :: zhalf, zfull
                0673   real,    intent(in), dimension(:,:,:) :: akm, akh
                0674   real,    intent(out),dimension(:,:)   :: h
                0675   integer, intent(in), optional, dimension(:,:)   :: kbot
                0676
                0677 !=======================================================================
                0678   real,    dimension(size(zfull,1),size(zfull,2)) :: km_surf, kh_surf
                0679   real,    dimension(size(zfull,1),size(zfull,2)) :: zhalfhalf
                0680   real, dimension(size(zfull,1),size(zfull,2),size(zfull,3))  :: zfull_ag
                0681   real, dimension(size(zfull,1),size(zfull,2),size(zfull,3)+1):: zhalf_ag
                0682   real, dimension(size(zfull,1),size(zfull,2),size(zfull,3)+1):: diff_tm
                0683   integer, dimension(size(zfull,1),size(zfull,2))            :: ibot
                0684   integer                                         :: i,j,k,nlev,nlat,nlon
                0685
                0686
                0687 nlev = size(zfull,3)
                0688 nlat = size(zfull,2)
                0689 nlon = size(zfull,1)
                0690
                0691 !compute height of surface
                0692 if (present(kbot)) then
                0693    ibot=kbot
                0694 else
                0695    ibot(:,:) = nlev
                0696 end if
                0697
                0698 !compute density profile, and heights relative to surface
                0699 do k = 1, nlev
                0700   zfull_ag(:,:,k) = zfull(:,:,k) - zsfc(:,:)
                0701   zhalf_ag(:,:,k) = zhalf(:,:,k) - zsfc(:,:)
                0702   end do
                0703 zhalf_ag(:,:,nlev+1) = zhalf(:,:,nlev+1) - zsfc(:,:)
                0704
                0705
                0706 !compute height half way between surface and lowest model level
                0707 zhalfhalf=0.5*MINVAL(MAX(zfull_ag,0.0))
                0708
                0709 !compute k's there by a call to mo_diff
                0710 call mo_diff(zhalfhalf,ustar,bstar,km_surf,kh_surf,1)
                0711
                0712 !create combined surface k's and diffusivity matrix
                0713 diff_tm(:,:,nlev+1) = 0.
                0714 diff_tm(:,:,1:nlev) = 0.5*(akm+akh)
                0715 if (present(kbot)) then
                0716 do j=1,nlat
                0717 do i=1,nlon
                0718    diff_tm(i,j,ibot(i,j)+1) = 0.5*(km_surf(i,j)+kh_surf(i,j))
                0719    zhalf_ag(i,j,ibot(i,j)+1) = zhalfhalf(i,j)
                0720 enddo
                0721 enddo
                0722 else
                0723    diff_tm(:,:,nlev+1) = 0.5*(km_surf(:,:)+kh_surf(:,:))
                0724    zhalf_ag(:,:,nlev+1) = zhalfhalf(:,:)
                0725 end if
                0726
                0727
                0728 !determine pbl depth as the height above ground where diff_tm
                0729 !first falls beneath a critical value kcrit.  If the value between
                0730 !ground and level 1 does not exceed kcrit set pbl depth equal to zero.
                0731 !kcrit is a namelist parameter.
                0732
                0733 do j = 1,nlat
                0734 do i = 1,nlon
                0735              if (diff_tm(i,j,ibot(i,j)+1) .gt. kcrit) then
                0736                        k=ibot(i,j)+1
                0737                        do while (k.gt. 2 .and. &
                0738                                  diff_tm(i,j,k-1).gt.kcrit)
                0739                              k=k-1
                0740                        enddo
                0741                        h(i,j) = zhalf_ag(i,j,k) + &
                0742                          (zhalf_ag(i,j,k-1)-zhalf_ag(i,j,k))* &
                0743                          (diff_tm(i,j,k)-kcrit) / &
                0744                          (diff_tm(i,j,k)-diff_tm(i,j,k-1))
                0745              else
                0746                        h(i,j) = 0.
                0747              end if
                0748 enddo
                0749 enddo
                0750
                0751 !-----------------------------------------------------------------------
                0752
                0753
                0754 !=====================================================================
                0755   end SUBROUTINE K_PBL_DEPTH
                0756
                0757 !#######################################################################
                0758
                0759  SUBROUTINE TKE_SURF ( u_star, TKE, kbot )
                0760
                0761 !=======================================================================
                0762 !---------------------------------------------------------------------
                0763 ! Arguments (Intent in)
                0764 !       u_star -  surface friction velocity (m/s)
                0765 !       kbot   -  OPTIONAL;lowest model level index (integer);
                0766 !                 at levels > Kbot, Mask = 0.
                0767 !---------------------------------------------------------------------
                0768   real, intent(in), dimension(:,:)   :: u_star
                0769
                0770   integer, intent(in), OPTIONAL, dimension(:,:) :: kbot
                0771
                0772 !---------------------------------------------------------------------
                0773 ! Arguments (Intent inout)
                0774 !       TKE  -  turbulent kinetic energy
                0775 !---------------------------------------------------------------------
                0776   real, intent(inout), dimension(:,:,:) :: TKE
                0777
                0778 !---------------------------------------------------------------------
                0779 !  (Intent local)
                0780 !---------------------------------------------------------------------
                0781   real, dimension(SIZE(u_star,1),SIZE(u_star,2)) :: x1
                0782   integer  :: ix, jx, kxp, i, j, k
                0783
                0784 !=======================================================================
                0785
                0786   ix  = SIZE( u_star, 1 )
                0787   jx  = SIZE( u_star, 2 )
                0788   kxp = SIZE( TKE,    3 )
                0789
                0790 !---------------------------------------------------------------------
                0791
                0792   x1 = bcq * u_star * u_star
                0793
                0794   if( PRESENT( kbot ) ) then
                0795     do j = 1,jx
                0796     do i = 1,ix
                0797        k = kbot(i,j) + 1
                0798       TKE(i,j,k) = x1(i,j)
                0799     end do
                0800     end do
                0801   else
                0802       TKE(:,:,kxp) = x1(:,:)
                0803   endif
                0804
                0805 !=======================================================================
                0806  end SUBROUTINE TKE_SURF
                0807
                0808 !#######################################################################
                0809   end MODULE MY25_TURB_MOD