pkg/atm_phys/diffusivity_mod.F90

b2ea1d2979 Jean*0001 module diffusivity_mod
                0002
                0003 !=======================================================================
                0004 !
                0005 !                          DIFFUSIVITY MODULE
                0006 !
                0007 !     Routines for computing atmospheric diffusivities in the
                0008 !       planetary boundary layer and in the free atmosphere
                0009 !
                0010 !=======================================================================
                0011
                0012 use gcm_params_mod, only: gcm_LEN_MBUF, gcm_SQZ_R, gcm_stdMsgUnit
                0013 use constants_mod, only : grav, vonkarm, cp_air, rdgas, rvgas
                0014
                0015 !use       fms_mod, only:  error_mesg, FATAL, file_exist,   &
                0016 !                          check_nml_error, open_namelist_file,      &
                0017 !                          mpp_pe, mpp_root_pe, close_file, stdlog,  &
                0018 !                          write_version_number
                0019
                0020 use monin_obukhov_mod, only : mo_diff
                0021
                0022 implicit none
                0023 private
                0024
                0025 ! public interfaces
                0026 !=======================================================================
                0027
                0028  public diffusivity, pbl_depth, molecular_diff
                0029
                0030 !=======================================================================
                0031
                0032 ! form of iterfaces
                0033
                0034 !=======================================================================
                0035 ! subroutine diffusivity (t, q, u, v, p_full, p_half, z_full, z_half,
                0036 !                         u_star, b_star, h, k_m, k_t)
                0037
                0038 ! input:
                0039
                0040 !        t     : real, dimension(:,:,:) -- (:,:,pressure), third index running
                0041 !                          from top of atmosphere to bottom
                0042 !                 temperature (K)
                0043 !
                0044 !        q     : real, dimension(:,:,:)
                0045 !                 water vapor specific humidity (nondimensional)
                0046 !
                0047 !        u     : real, dimension(:,:)
                0048 !                 zonal wind (m/s)
                0049 !
                0050 !        v     : real, dimension(:,:,:)
                0051 !                 meridional wind (m/s)
                0052 !
                0053 !        z_full  : real, dimension(:,:,:
                0054 !                 height of full levels (m)
                0055 !                 1 = top of atmosphere; size(p_half,3) = surface
                0056 !                 size(z_full,3) = size(t,3)
                0057 !
                0058 !        z_half  : real, dimension(:,:,:)
                0059 !                 height of  half levels (m)
                0060 !                 size(z_half,3) = size(t,3) +1
                0061 !              z_half(:,:,size(z_half,3)) must be height of surface!
                0062 !                                  (if you are not using eta-model)
                0063 !
                0064 !        u_star: real, dimension(:,:)
                0065 !                friction velocity (m/s)
                0066 !
                0067 !        b_star: real, dimension(:,:)
                0068 !                buoyancy scale (m/s**2)
                0069
                0070 !   (u_star and b_star can be obtained by calling
                0071 !     mo_drag in monin_obukhov_mod)
                0072
                0073 ! output:
                0074
                0075 !        h     : real, dimension(:,:,)
                0076 !                 depth of planetary boundary layer (m)
                0077 !
                0078 !        k_m   : real, dimension(:,:,:)
                0079 !                diffusivity for momentum (m**2/s)
                0080 !
                0081 !                defined at half-levels
                0082 !                size(k_m,3) should be at least as large as size(t,3)
                0083 !                only the returned values at
                0084 !                      levels 2 to size(t,3) are meaningful
                0085 !                other values will be returned as zero
                0086 !
                0087 !        k_t   : real, dimension(:,:,:)
                0088 !                diffusivity for temperature and scalars (m**2/s)
                0089 !
                0090 !
                0091 !=======================================================================
                0092
                0093 !--------------------- version number ----------------------------------
                0094
15388b052e Jean*0095 character(len=128) :: version = '$Id: diffusivity_mod.F90,v 1.2 2017/08/11 20:48:50 jmc Exp $'
b2ea1d2979 Jean*0096 character(len=128) :: tag = '$Name:  $'
                0097
                0098 !=======================================================================
                0099
                0100 !  DEFAULT VALUES OF NAMELIST PARAMETERS:
                0101
                0102 logical :: fixed_depth         = .false.
                0103 logical :: do_virtual_non_mcm  = .false.  ! applies only to non-'mcm' pbl scheme
                0104 real    :: depth_0             =  5000.0
                0105 real    :: frac_inner          =  0.1
                0106 real    :: rich_crit_pbl       =  1.0
                0107 real    :: entr_ratio          =  0.2
                0108 real    :: parcel_buoy         =  2.0
                0109 real    :: znom                =  1000.0
                0110 logical :: free_atm_diff       = .false.
                0111 logical :: free_atm_skyhi_diff = .false.
                0112 logical :: pbl_mcm             = .false.
                0113 real    :: rich_crit_diff      =  0.25
                0114 real    :: mix_len             = 30.
                0115 real    :: rich_prandtl        =  1.00
                0116 real    :: background_m        =  0.0
                0117 real    :: background_t        =  0.0
                0118 logical :: ampns               = .false. ! include delta z factor in
                0119                                          ! defining ri ?
                0120 real    :: ampns_max           = 1.0E20  ! limit to reduction factor
                0121                                          ! applied to ri due to delta z
                0122                                          ! factor
                0123
                0124 namelist /diffusivity_nml/ fixed_depth, depth_0, frac_inner,&
                0125                            rich_crit_pbl, entr_ratio, parcel_buoy,&
                0126                            znom, free_atm_diff, free_atm_skyhi_diff,&
                0127                            pbl_mcm, rich_crit_diff, mix_len, rich_prandtl,&
                0128                            background_m, background_t, ampns, ampns_max, &
                0129                            do_virtual_non_mcm
                0130
                0131 !=======================================================================
                0132
                0133 !  OTHER MODULE VARIABLES
                0134
                0135 real    :: small  = 1.e-04
                0136 real    :: gcp    = grav/cp_air
                0137 logical :: init   = .false.
                0138 real    :: beta   = 1.458e-06
                0139 real    :: rbop1  = 110.4
                0140 real    :: rbop2  = 1.405
                0141
                0142 real, parameter :: d608 = (rvgas-rdgas)/rdgas
                0143
                0144 contains
                0145
                0146 !=======================================================================
                0147
                0148 subroutine diffusivity_init(myThid)
                0149
                0150 integer, intent(in) :: myThid
                0151 integer :: unit, ierr, io
                0152
                0153 integer         :: iUnit
                0154 CHARACTER*(gcm_LEN_MBUF) :: msgBuf
                0155
                0156 !------------------- read namelist input -------------------------------
                0157
                0158 ! read namelist and copy to logfile
                0159
                0160 !    _BARRIER
                0161 !    _BEGIN_MASTER(myThid)
                0162      CALL BARRIER(myThid)
                0163      IF ( myThid.EQ.1 ) THEN
                0164
                0165      WRITE(msgBuf,'(A)') 'DIFFUSIVITY_INIT: opening data.atm_gray'
                0166      CALL PRINT_MESSAGE( msgBuf, gcm_stdMsgUnit, gcm_SQZ_R, myThid )
                0167      CALL OPEN_COPY_DATA_FILE(                                      &
                0168                            'data.atm_gray', 'DIFFUSIVITY_INIT',       &
                0169                            iUnit,                                   &
                0170                            myThid )
                0171 !    Read parameters from open data file
                0172      READ(UNIT=iUnit,NML=diffusivity_nml)
                0173      WRITE(msgBuf,'(A)')                                            &
                0174           'DIFFUSIVITY_INIT: finished reading data.atm_gray'
                0175      CALL PRINT_MESSAGE( msgBuf, gcm_stdMsgUnit, gcm_SQZ_R, myThid )
                0176 !    Close the open data file
15388b052e Jean*0177 #ifdef SINGLE_DISK_IO
b2ea1d2979 Jean*0178      CLOSE(iUnit)
15388b052e Jean*0179 #else
                0180      CLOSE(iUnit,STATUS='DELETE')
                0181 #endif /* SINGLE_DISK_IO */
b2ea1d2979 Jean*0182
                0183 !     if (file_exist('input.nml')) then
                0184 !        unit = open_namelist_file ()
                0185 !        ierr=1; do while (ierr /= 0)
                0186 !           read  (unit, nml=diffusivity_nml, iostat=io, end=10)
                0187 !           ierr = check_nml_error(io,'diffusivity_nml')
                0188 !        enddo
                0189 ! 10     call close_file (unit)
                0190
                0191 !------------------- dummy checks --------------------------------------
                0192          if (frac_inner .le. 0. .or. frac_inner .ge. 1.) &
                0193             CALL PRINT_ERROR( 'diffusivity_init'//  &
                0194             'frac_inner must be between 0 and 1', myThid)
                0195 !           call error_mesg ('diffusivity_init',  &
                0196 !           'frac_inner must be between 0 and 1', FATAL)
                0197          if (rich_crit_pbl .lt. 0.) &
                0198             CALL PRINT_ERROR( 'diffusivity_init'//  &
                0199            'rich_crit_pbl must be greater than or equal to zero', myThid )
                0200 !           call error_mesg ('diffusivity_init',  &
                0201 !          'rich_crit_pbl must be greater than or equal to zero', FATAL)
                0202          if (entr_ratio .lt. 0.) &
                0203             CALL PRINT_ERROR( 'diffusivity_init'//  &
                0204             'entr_ratio must be greater than or equal to zero', myThid )
                0205 !           call error_mesg ('diffusivity_init',  &
                0206 !           'entr_ratio must be greater than or equal to zero', FATAL)
                0207          if (znom .le. 0.) &
                0208             CALL PRINT_ERROR( 'diffusivity_init'//  &
                0209             'znom must be greater than zero', myThid )
                0210 !           call error_mesg ('diffusivity_init',  &
                0211 !           'znom must be greater than zero', FATAL)
                0212          if (.not.free_atm_diff .and. free_atm_skyhi_diff)&
                0213             CALL PRINT_ERROR( 'diffusivity_init'//  &
                0214             'free_atm_diff must be set to true if '//&
                0215             'free_atm_skyhi_diff = .true.', myThid )
                0216 !           call error_mesg ('diffusivity_init',  &
                0217 !           'free_atm_diff must be set to true if '//&
                0218 !           'free_atm_skyhi_diff = .true.', FATAL)
                0219          if (rich_crit_diff .le. 0.) &
                0220             CALL PRINT_ERROR( 'diffusivity_init'//  &
                0221             'rich_crit_diff must be greater than zero', myThid )
                0222 !           call error_mesg ('diffusivity_init',  &
                0223 !           'rich_crit_diff must be greater than zero', FATAL)
                0224          if (mix_len .lt. 0.) &
                0225             CALL PRINT_ERROR( 'diffusivity_init'//  &
                0226             'mix_len must be greater than or equal to zero', myThid )
                0227 !           call error_mesg ('diffusivity_init',  &
                0228 !           'mix_len must be greater than or equal to zero', FATAL)
                0229          if (rich_prandtl .lt. 0.) &
                0230             CALL PRINT_ERROR( 'diffusivity_init'//  &
                0231             'rich_prandtl must be greater than or equal to zero', myThid )
                0232 !           call error_mesg ('diffusivity_init',  &
                0233 !           'rich_prandtl must be greater than or equal to zero', FATAL)
                0234          if (background_m .lt. 0.) &
                0235             CALL PRINT_ERROR( 'diffusivity_init'//  &
                0236             'background_m must be greater than or equal to zero', myThid )
                0237 !           call error_mesg ('diffusivity_init',  &
                0238 !           'background_m must be greater than or equal to zero', FATAL)
                0239          if (background_t .lt. 0.) &
                0240             CALL PRINT_ERROR( 'diffusivity_init'//  &
                0241             'background_t must be greater than or equal to zero', myThid )
                0242 !           call error_mesg ('diffusivity_init',  &
                0243 !           'background_t must be greater than or equal to zero', FATAL)
                0244          if (ampns_max .lt. 1.) &
                0245             CALL PRINT_ERROR( 'diffusivity_init'//  &
                0246             'ampns_max must be greater than or equal to one', myThid )
                0247 !           call error_mesg ('diffusivity_init',  &
                0248 !           'ampns_max must be greater than or equal to one', FATAL)
                0249          if (ampns .and. .not. free_atm_skyhi_diff) &
                0250             CALL PRINT_ERROR( 'diffusivity_init'//  &
                0251             'ampns is only valid when free_atm_skyhi_diff is & also true', myThid )
                0252 !           call error_mesg ('diffusivity_init',  &
                0253 !           'ampns is only valid when free_atm_skyhi_diff is &
                0254 !                  & also true', FATAL)
                0255
                0256 !     endif  !end of reading input.nml
                0257
                0258 !---------- output namelist to log-------------------------------------
                0259
                0260 !     call write_version_number(version, tag)
                0261 !     if ( mpp_pe() == mpp_root_pe() ) write (stdlog(), nml=diffusivity_nml)
                0262
                0263       init = .true.
                0264
                0265      ENDIF
                0266      CALL BARRIER(myThid)
                0267
                0268 return
                0269 end subroutine diffusivity_init
                0270
                0271 !=======================================================================
                0272
                0273 subroutine diffusivity(t, q, u, v, p_full, p_half, z_full, z_half,  &
                0274                        u_star, b_star, h, k_m, k_t, myThid, kbot)
                0275
                0276 real,    intent(in),           dimension(:,:,:) :: t, q, u, v
                0277 real,    intent(in),           dimension(:,:,:) :: p_full, p_half
                0278 real,    intent(in),           dimension(:,:,:) :: z_full, z_half
                0279 real,    intent(in),           dimension(:,:)   :: u_star, b_star
                0280 real,    intent(inout),        dimension(:,:,:) :: k_m, k_t
                0281 real,    intent(out),          dimension(:,:)   :: h
                0282 integer, intent (in)                            :: myThid
                0283 integer, intent(in), optional, dimension(:,:)   :: kbot
                0284
                0285 real, dimension(size(t,1),size(t,2),size(t,3))  :: svcp,z_full_ag, &
                0286                                                    k_m_save, k_t_save
                0287 real, dimension(size(t,1),size(t,2),size(t,3)+1):: z_half_ag
                0288 real, dimension(size(t,1),size(t,2))            :: z_surf
                0289 integer                                         :: i,j,k,nlev,nlat,nlon
                0290
                0291 if(.not.init) call diffusivity_init(myThid)
                0292
                0293 nlev = size(t,3)
                0294
                0295 k_m_save = k_m
                0296 k_t_save = k_t
                0297
                0298 !compute height of surface
                0299 if (present(kbot)) then
                0300    nlat = size(t,2)
                0301    nlon = size(t,1)
                0302    do j=1,nlat
                0303    do i=1,nlon
                0304           z_surf(i,j) = z_half(i,j,kbot(i,j)+1)
                0305    enddo
                0306    enddo
                0307 else
                0308    z_surf(:,:) = z_half(:,:,nlev+1)
                0309 end if
                0310
                0311 !compute density profile, and heights relative to surface
                0312 do k = 1, nlev
                0313
                0314   z_full_ag(:,:,k) = z_full(:,:,k) - z_surf(:,:)
                0315   z_half_ag(:,:,k) = z_half(:,:,k) - z_surf(:,:)
                0316
                0317   if (do_virtual_non_mcm) then
                0318    svcp(:,:,k)  =   t(:,:,k)*(1. + d608*q(:,:,k)) + gcp*(z_full_ag(:,:,k))
                0319   else
                0320    svcp(:,:,k)  =   t(:,:,k) + gcp*(z_full_ag(:,:,k))
                0321   endif
                0322
                0323 end do
                0324 z_half_ag(:,:,nlev+1) = z_half(:,:,nlev+1) - z_surf(:,:)
                0325
                0326 if(fixed_depth)  then
                0327    h = depth_0
                0328 else
                0329    call pbl_depth(svcp,u,v,z_full_ag,u_star,b_star,h, myThid, kbot=kbot)
                0330 end if
                0331
                0332 if(pbl_mcm) then
                0333    call diffusivity_pbl_mcm (u,v, t, p_full, p_half, &
                0334                              z_full_ag, z_half_ag, h, k_m, k_t, myThid)
                0335 else
                0336    call diffusivity_pbl  (svcp, u, v, z_half_ag, h, u_star, b_star,&
                0337                        k_m, k_t, myThid, kbot=kbot)
                0338 end if
                0339 if(free_atm_diff) &
                0340    call diffusivity_free (svcp, u, v, z_full_ag, z_half_ag, h, k_m, k_t, myThid)
                0341
                0342 k_m = k_m + k_m_save
                0343 k_t = k_t + k_t_save
                0344
                0345 !NOTE THAT THIS LINE MUST FOLLOW DIFFUSIVITY_FREE SO THAT ENTRAINMENT
                0346 !K's DO NOT GET OVERWRITTEN IN DIFFUSIVITY_FREE SUBROUTINE
                0347 if(entr_ratio .gt. 0. .and. .not. fixed_depth) &
                0348     call diffusivity_entr(svcp,z_full_ag,h,u_star,b_star,k_m,k_t, myThid)
                0349
                0350 !set background diffusivities
                0351 if(background_m.gt.0.0) k_m = max(k_m,background_m)
                0352 if(background_t.gt.0.0) k_t = max(k_t,background_t)
                0353
                0354 return
                0355 end subroutine diffusivity
                0356
                0357 !=======================================================================
                0358
                0359 subroutine pbl_depth(t, u, v, z, u_star, b_star, h, myThid, kbot)
                0360
                0361 real,   intent(in) ,           dimension(:,:,:) :: t, u, v, z
                0362 real,   intent(in) ,           dimension(:,:)   :: u_star,b_star
                0363 real,   intent(out),           dimension(:,:)   :: h
                0364 integer, intent(in)                             :: myThid
                0365 integer,intent(in) , optional, dimension(:,:)   :: kbot
                0366
                0367 real,    dimension(size(t,1),size(t,2),size(t,3))  :: rich
                0368 real,    dimension(size(t,1),size(t,2))            :: ws,k_t_ref,&
                0369                                                       h_inner,tbot
                0370 real                                               :: rich1, rich2,&
                0371                                                       h1,h2,svp,t1,t2
                0372 integer, dimension(size(t,1),size(t,2))            :: ibot
                0373 integer                                            :: i,j,k,nlon,&
                0374                                                       nlat, nlev
                0375
                0376 nlev = size(t,3)
                0377 nlat = size(t,2)
                0378 nlon = size(t,1)
                0379
                0380 !assign ibot, compute tbot (virtual temperature at lowest level)
                0381 if (present(kbot)) then
                0382     ibot(:,:) = kbot
                0383     do j = 1,nlat
                0384     do i = 1,nlon
                0385           tbot(i,j) = t(i,j,ibot(i,j))
                0386     enddo
                0387     enddo
                0388 else
                0389     ibot(:,:) = nlev
                0390     tbot(:,:) = t(:,:,nlev)
                0391 end if
                0392
                0393 !compute richardson number for use in pbl depth of neutral/stable side
                0394 do k = 1,nlev
                0395   rich(:,:,k) =  z(:,:,k)*grav*(t(:,:,k)-tbot(:,:))/tbot(:,:)&
                0396                 /(u(:,:,k)*u(:,:,k) + v(:,:,k)*v(:,:,k) + small )
                0397 end do
                0398
                0399 !compute ws to be used in evaluating parcel buoyancy
                0400 !ws = u_star / phi(h_inner,u_star,b_star)  .  To find phi
                0401 !a call to mo_diff is made.
                0402
                0403 h_inner(:,:)=frac_inner*znom
                0404 call mo_diff(h_inner, u_star, b_star, ws, k_t_ref, myThid )
                0405 ws = max(small,ws/vonkarm/h_inner)
                0406
                0407 do j = 1, nlat
                0408  do i = 1, nlon
                0409
                0410         !do neutral or stable case
                0411         if (b_star(i,j).le.0.) then
                0412
                0413               h1     = z(i,j,ibot(i,j))
                0414               h(i,j) = h1
                0415               rich1  = rich(i,j,ibot(i,j))
                0416               do k = ibot(i,j)-1, 1, -1
                0417                        rich2 = rich(i,j,k)
                0418                        h2    = z(i,j,k)
                0419                        if(rich2.gt.rich_crit_pbl) then
                0420                              h(i,j) = h2 + (h1 - h2)*(rich2 - rich_crit_pbl)&
                0421                                                     /(rich2 - rich1        )
                0422                              go to 10
                0423                        endif
                0424                        rich1 = rich2
                0425                        h1    = h2
                0426               enddo
                0427
                0428         !do unstable case
                0429         else
                0430
                0431               svp    = tbot(i,j)*(1.+ &
                0432                        (parcel_buoy*u_star(i,j)*b_star(i,j)/grav/ws(i,j)) )
                0433               h1     = z(i,j,ibot(i,j))
                0434               h(i,j) = h1
                0435               t1     = tbot(i,j)
                0436               do k = ibot(i,j)-1 , 1, -1
                0437                        h2 = z(i,j,k)
                0438                        t2 = t(i,j,k)
                0439                        if (t2.gt.svp) then
                0440                              h(i,j) = h2 + (h1 - h2)*(t2 - svp)/(t2 - t1 )
                0441                              go to 10
                0442                        end if
                0443                        h1 = h2
                0444                        t1 = t2
                0445               enddo
                0446
                0447         end if
                0448 10 continue
                0449   enddo
                0450 enddo
                0451
                0452 return
                0453 end subroutine pbl_depth
                0454
                0455 !=======================================================================
                0456
                0457 subroutine diffusivity_pbl(t, u, v, z_half, h, u_star, b_star, &
                0458                            k_m, k_t, myThid, kbot)
                0459
                0460 real,    intent(in)  ,           dimension(:,:,:) :: t, u, v, z_half
                0461 real,    intent(in)  ,           dimension(:,:)   :: h, u_star, b_star
                0462 real,    intent(inout) ,           dimension(:,:,:) :: k_m, k_t
                0463 integer, intent (in)                                :: myThid
                0464 integer, intent(in)  , optional, dimension(:,:)   :: kbot
                0465
                0466 real, dimension(size(t,1),size(t,2))              :: h_inner, k_m_ref,&
                0467                                                      k_t_ref, factor
                0468 real, dimension(size(t,1),size(t,2),size(t,3)+1)  :: zm
                0469 real                                              :: h_inner_max
                0470 integer                                           :: i,j, k, kk, nlev
                0471
                0472 nlev = size(t,3)
                0473
                0474 !assign z_half to zm, and set to zero any values of zm < 0.
                0475 !the setting to zero is necessary so that when using eta model
                0476 !below ground half levels will have zero k_m and k_t
                0477 zm = z_half
                0478 if (present(kbot)) then
                0479    where(zm < 0.)
                0480         zm = 0.
                0481    end where
                0482 end if
                0483
                0484 h_inner    = frac_inner*h
                0485 h_inner_max = maxval(h_inner)
                0486
                0487 kk = nlev
                0488 do k = 2, nlev
                0489   if( minval(zm(:,:,k)) < h_inner_max) then
                0490       kk = k
                0491       exit
                0492   end if
                0493 end do
                0494
                0495 k_m = 0.0
                0496 k_t = 0.0
                0497
                0498 call mo_diff(h_inner        , u_star, b_star, k_m_ref         , k_t_ref, myThid )
                0499 call mo_diff(zm(:,:,kk:nlev), u_star, b_star, k_m(:,:,kk:nlev), k_t(:,:,kk:nlev), &
                0500               myThid )
                0501
                0502 do k = 2, nlev
                0503   where(zm(:,:,k) >= h_inner .and. zm(:,:,k) < h)
                0504     factor = (zm(:,:,k)/h_inner)* &
                0505              (1.0 - (zm(:,:,k) - h_inner)/(h - h_inner))**2
                0506     k_m(:,:,k) = k_m_ref*factor
                0507     k_t(:,:,k) = k_t_ref*factor
                0508   end where
                0509 end do
                0510
                0511 return
                0512 end subroutine diffusivity_pbl
                0513
                0514 !=======================================================================
                0515
                0516 subroutine diffusivity_pbl_mcm(u, v, t, p_full, p_half, z_full, z_half, &
                0517                                h, k_m, k_t, myThid)
                0518
                0519 real, intent(in)  , dimension(:,:,:) :: u, v, t, z_full, z_half
                0520 real, intent(in)  , dimension(:,:,:) :: p_full, p_half
                0521 real, intent(in)  , dimension(:,:)   :: h
                0522 real, intent(inout) , dimension(:,:,:) :: k_m, k_t
                0523 integer, intent(in)                    :: myThid
                0524
                0525 integer                                        :: k, nlev
                0526 real, dimension(size(z_full,1),size(z_full,2)) :: elmix, htcrit
                0527 real, dimension(size(z_full,1),size(z_full,2)) :: delta_u, delta_v, delta_z
                0528
                0529 real :: htcrit_ss
                0530 real :: h_ss
                0531 real, dimension(size(z_full,1),size(z_full,2)) :: sig_half, z_half_ss, elmix_ss
                0532
                0533 !  htcrit_ss = height at which mixing length is a maximum (75m)
                0534
                0535 !  h_ss   = height at which mixing length vanishes (4900m)
                0536 !  elmix_ss   = mixing length
                0537
                0538 ! Define some constants:
                0539 !  salaps = standard atmospheric lapse rate (K/m)
                0540 !  tsfc   = idealized global mean surface temperature (15C)
                0541 real :: tsfc = 288.16
                0542 real :: salaps = -6.5e-3
                0543
                0544 nlev = size(z_full,3)
                0545
                0546 k_m = 0.
                0547
                0548 h_ss = depth_0
                0549 htcrit_ss = frac_inner*h_ss
                0550
                0551 do k = 2, nlev
                0552
                0553 ! TK mods 8/13/01:  (code derived from SS)
                0554 ! Compute the height of each half level assuming a constant
                0555 ! standard lapse rate using the above procedure.
                0556 ! WARNING: These should be used with caution.  They will
                0557 !  have large errors above the tropopause.
                0558
                0559 ! In order to determine the height, the layer mean temperature
                0560 ! from the surface to that level is required.  A surface
                0561 ! temperature of 15 deg Celsius and a standard lapse rate of
                0562 ! -6.5 deg/km will be used to estimate an average temperature
                0563 ! profile.
                0564
                0565    sig_half = p_half(:,:,k)/p_half(:,:,nlev+1)
                0566    z_half_ss = -rdgas * .5*(tsfc+tsfc*(sig_half**(-rdgas*salaps/grav))) * alog(sig_half)/grav
                0567
                0568    !compute mixing length as in SS (no geographical variation)
                0569     elmix_ss = 0.
                0570
                0571     where (z_half_ss < htcrit_ss .and. z_half_ss > 0.)
                0572          elmix_ss = vonkarm*z_half_ss
                0573     endwhere
                0574     where (z_half_ss >= htcrit_ss .and. z_half_ss < h_ss)
                0575          elmix_ss = vonkarm*htcrit_ss*(h_ss-z_half_ss)/(h_ss-htcrit_ss)
                0576     endwhere
                0577
                0578    delta_z = rdgas*0.5*(t(:,:,k)+t(:,:,k-1))*(p_full(:,:,k)-p_full(:,:,k-1))/&
                0579              (grav*p_half(:,:,k))
                0580    delta_u =      u(:,:,k-1) -      u(:,:,k)
                0581    delta_v =      v(:,:,k-1) -      v(:,:,k)
                0582
                0583    k_m(:,:,k) =   elmix_ss * elmix_ss *&
                0584                   sqrt(delta_u*delta_u + delta_v*delta_v)/delta_z
                0585
                0586 end do
                0587
                0588 k_t = k_m
                0589
                0590 return
                0591 end subroutine diffusivity_pbl_mcm
                0592
                0593 !=======================================================================
                0594
                0595 subroutine diffusivity_free(t, u, v, z, zz, h, k_m, k_t, myThid)
                0596
                0597 real, intent(in)    , dimension(:,:,:) :: t, u, v, z, zz
                0598 real, intent(in)    , dimension(:,:)   :: h
                0599 real, intent(inout) , dimension(:,:,:) :: k_m, k_t
                0600 integer, intent(in)                    :: myThid
                0601
                0602 real, dimension(size(t,1),size(t,2))   :: dz, b, speed2, rich, fri, &
                0603                                           alpz, fri2
                0604 integer                                :: k
                0605
                0606 do k = 2, size(t,3)
                0607
                0608 !----------------------------------------------------------------------
                0609 !  define the richardson number. set it to zero if it is negative. save
                0610 !  a copy of it for later use (rich2).
                0611 !----------------------------------------------------------------------
                0612   dz     = z(:,:,k-1) - z(:,:,k)
                0613   b      = grav*(t(:,:,k-1)-t(:,:,k))/t(:,:,k)
                0614   speed2 = (u(:,:,k-1) - u(:,:,k))**2 + (v(:,:,k-1) - v(:,:,k))**2
                0615   rich= b*dz/(speed2+small)
                0616   rich = max(rich, 0.0)
                0617
                0618   if (free_atm_skyhi_diff) then
                0619 !---------------------------------------------------------------------
                0620 !   limit the standard richardson number to between 0 and the critical
                0621 !   value (rich2). compute the richardson number factor needed in the
                0622 !   eddy mixing coefficient using this standard richardson number.
                0623 !---------------------------------------------------------------------
                0624     where (rich(:,:) >= rich_crit_diff)
                0625       fri2(:,:) = 0.0
                0626     elsewhere
                0627       fri2(:,:)  = (1.0 - rich/rich_crit_diff)**2
                0628     endwhere
                0629   endif
                0630
                0631 !---------------------------------------------------------------------
                0632 !  if ampns is activated, compute the delta z factor. define rich
                0633 !  including this factor.
                0634 !---------------------------------------------------------------------
                0635   if (ampns) then
                0636     alpz(:,:) = MIN ( (1.  + 1.e-04*(dz(:,:)**1.5)), ampns_max)
                0637     rich(:,:) = rich(:,:) / alpz(:,:)
                0638   endif
                0639
                0640 !---------------------------------------------------------------------
                0641 !   compute the richardson number factor to be used in the eddy
                0642 !   mixing coefficient. if ampns is on, this value includes it; other-
                0643 !   wise it does not.
                0644 !---------------------------------------------------------------------
                0645   fri(:,:)   = (1.0 - rich/rich_crit_diff)**2
                0646
                0647 !---------------------------------------------------------------------
                0648 !   compute the eddy mixing coefficients in the free atmosphere ( zz
                0649 !   > h). in the non-ampns case, values are obtained only when the
                0650 !   standard richardson number is sub-critical; in the ampns case values
                0651 !   are obtained only when the richardson number computed with the
                0652 !   ampns factor is sub critical. when the ampns factor is activated,
                0653 !   it is also included in the mixing coefficient. the value of mixing
                0654 !   for temperature, etc. is reduced dependent on the ri stability
                0655 !   factor calculated without the ampns factor.
                0656 !---------------------------------------------------------------------
                0657   if (free_atm_skyhi_diff) then
                0658
                0659 !---------------------------------------------------------------------
                0660 !   this is the skyhi-like formulation -- possible ampns factor, ratio
                0661 !   of k_m to k_t defined based on computed stability factor.
                0662 !---------------------------------------------------------------------
                0663     if (ampns) then
                0664       where (rich < rich_crit_diff .and. zz(:,:,k) > h)
                0665            k_m(:,:,k) = mix_len*mix_len*sqrt(speed2)*fri(:,:)* &
                0666                         ( 1.  + 1.e-04*(dz(:,:)**1.5))/dz
                0667            k_t(:,:,k) = k_m(:,:,k)* (0.1 + 0.9*fri2(:,:))
                0668       end where
                0669     else
                0670       where (rich < rich_crit_diff .and. zz(:,:,k) > h)
                0671         k_m(:,:,k) = mix_len*mix_len*sqrt(speed2)*fri(:,:)/dz
                0672         k_t(:,:,k) = k_m(:,:,k)* (0.1 + 0.9*fri2(:,:))
                0673       end where
                0674     endif
                0675   else
                0676
                0677 !---------------------------------------------------------------------
                0678 !   this is the non-skyhi-like formulation -- no ampns factor, ratio
                0679 !   of k_m to k_t defined by rich_prandtl.
                0680 !---------------------------------------------------------------------
                0681     where (rich < rich_crit_diff .and. zz(:,:,k) > h)
                0682          k_t(:,:,k) = mix_len*mix_len*sqrt(speed2)*fri(:,:)/dz
                0683          k_m(:,:,k) = k_t(:,:,k)*rich_prandtl
                0684     end where
                0685   end if
                0686 end do
                0687
                0688 end subroutine diffusivity_free
                0689
                0690 !=======================================================================
                0691
                0692 subroutine molecular_diff ( temp, press, k_m, k_t, myThid)
                0693
                0694 real, intent(in),    dimension (:,:,:)  ::  temp, press
                0695 real, intent(inout), dimension (:,:,:)  ::  k_m, k_t
                0696 integer, intent(in)                     :: myThid
                0697
                0698       real, dimension (size(temp,1), size(temp,2)) :: temp_half, &
                0699                                                       rho_half, rbop2d
                0700       integer      :: k
                0701
                0702 !---------------------------------------------------------------------
                0703
                0704       do k=2,size(temp,3)
                0705         temp_half(:,:) = 0.5*(temp(:,:,k) + temp(:,:,k-1))
                0706         rho_half(:,:) = press(:,:,k)/(rdgas*temp_half(:,:) )
                0707         rbop2d(:,:)  = beta*temp_half(:,:)*sqrt(temp_half(:,:))/  &
                0708                        (rho_half(:,:)*(temp_half(:,:)+rbop1))
                0709         k_m(:,:,k) = rbop2d(:,:)
                0710         k_t(:,:,k) = rbop2d(:,:)*rbop2
                0711       end do
                0712
                0713       k_m(:,:,1) = 0.0
                0714       k_t(:,:,1) = 0.0
                0715
                0716 end subroutine molecular_diff
                0717
                0718 !=======================================================================
                0719
                0720 subroutine diffusivity_entr(t, z,  h, u_star, b_star, k_m, k_t, myThid)
                0721
                0722 real, intent(in)    , dimension(:,:,:) :: t, z
                0723 real, intent(in)    , dimension(:,:)   :: h, u_star, b_star
                0724 real, intent(inout) , dimension(:,:,:) :: k_m, k_t
                0725 integer, intent(in)                    :: myThid
                0726
                0727 integer                                :: k, nlev
                0728
                0729 nlev=size(t,3)
                0730
                0731 do k = 2,nlev
                0732     where (b_star .gt. 0. .and. z(:,:,k-1) .gt. h .and. &
                0733                                 z(:,:,k)   .le. h)
                0734         k_t(:,:,k) = (z(:,:,k-1)-z(:,:,k))*entr_ratio*t(:,:,k)* &
                0735                       u_star*b_star/grav/max(small,t(:,:,k-1)-t(:,:,k))
                0736         k_m(:,:,k) = k_t(:,:,k)
                0737     end where
                0738 enddo
                0739 end subroutine diffusivity_entr
                0740
                0741 !=======================================================================
                0742
                0743 end module diffusivity_mod
                0744