pkg/atm_phys/monin_obukhov_mod.F90

b2ea1d2979 Jean*0001 module monin_obukhov_mod
                0002
                0003 !=======================================================================
                0004 !
                0005 !                         MONIN-OBUKHOV MODULE
                0006 !
                0007 !          Routines for computing surface drag coefficients
                0008 !                 from data at the lowest model level
                0009 !              and for computing the profile of fields
                0010 !           between the lowest model level and the ground
                0011 !                  using Monin-Obukhov scaling
                0012 !
                0013 !=======================================================================
                0014
                0015 use gcm_params_mod, only: gcm_LEN_MBUF, gcm_SQZ_R, gcm_stdMsgUnit
                0016 use constants_mod, only : grav, vonkarm
                0017 !use       fms_mod, only:  error_mesg, FATAL, file_exist,   &
                0018 !                          check_nml_error, open_namelist_file,      &
                0019 !                          mpp_pe, mpp_root_pe, close_file, stdlog, &
                0020 !                          write_version_number
                0021
                0022 implicit none
                0023 private
                0024
                0025 !=======================================================================
                0026  public mo_drag
                0027  public mo_profile
                0028  public mo_diff
                0029  public stable_mix
                0030 !=======================================================================
                0031
                0032 interface mo_drag
                0033     module procedure  mo_drag_0d, mo_drag_1d, mo_drag_2d
                0034 end interface
                0035
                0036 interface mo_profile
                0037     module procedure  mo_profile_0d,   mo_profile_1d,   mo_profile_2d, &
                0038                       mo_profile_0d_n, mo_profile_1d_n, mo_profile_2d_n
                0039 end interface
                0040
                0041 interface mo_diff
                0042     module procedure  mo_diff_0d_n, mo_diff_0d_1, &
                0043                       mo_diff_1d_n, mo_diff_1d_1, &
                0044                       mo_diff_2d_n, mo_diff_2d_1
                0045 end interface
                0046
                0047 interface stable_mix
                0048     module procedure  stable_mix_0d, stable_mix_1d, &
                0049                       stable_mix_2d, stable_mix_3d
                0050 end interface
                0051
                0052 !--------------------- version number ---------------------------------
                0053
15388b052e Jean*0054 character(len=128) :: version = '$Id: monin_obukhov_mod.F90,v 1.2 2017/08/11 20:48:51 jmc Exp $'
b2ea1d2979 Jean*0055 character(len=128) :: tagname = '$Name:  $'
                0056
                0057 !=======================================================================
                0058
                0059 !  DEFAULT VALUES OF NAMELIST PARAMETERS:
                0060
                0061 real    :: rich_crit      = 2.0
                0062 real    :: drag_min       = 1.e-05
                0063 logical :: neutral        = .false.
                0064 integer :: stable_option  = 1
                0065 real    :: zeta_trans     = 0.5
                0066
                0067 namelist /monin_obukhov_nml/ rich_crit, neutral, drag_min, &
                0068                              stable_option, zeta_trans
                0069
                0070 !=======================================================================
                0071
                0072 !  MODULE VARIABLES
                0073
                0074 real, parameter    :: small  = 1.e-04
                0075 real               :: b_stab, r_crit, sqrt_drag_min, lambda, rich_trans
                0076 logical            :: init = .false.
                0077
                0078 contains
                0079
                0080 !=======================================================================
                0081
                0082 subroutine monin_obukhov_init(myThid)
                0083
                0084 integer, intent(in) :: myThid
                0085 integer :: unit, ierr, io
                0086
                0087 integer         :: iUnit
                0088 CHARACTER*(gcm_LEN_MBUF) :: msgBuf
                0089 !------------------- read namelist input -------------------------------
                0090
                0091 ! read namelist and copy to logfile
                0092
                0093 !    _BARRIER
                0094 !    _BEGIN_MASTER(myThid)
                0095      CALL BARRIER(myThid)
                0096      IF ( myThid.EQ.1 ) THEN
                0097
                0098      WRITE(msgBuf,'(A)') 'MONIN_OBUKHOV_INIT: opening data.atm_gray'
                0099      CALL PRINT_MESSAGE( msgBuf, gcm_stdMsgUnit, gcm_SQZ_R, myThid )
                0100      CALL OPEN_COPY_DATA_FILE(                                      &
                0101                            'data.atm_gray', 'MONIN_OBUKHOV_INIT',       &
                0102                            iUnit,                                   &
                0103                            myThid )
                0104 !    Read parameters from open data file
                0105      READ(UNIT=iUnit,NML=monin_obukhov_nml)
                0106      WRITE(msgBuf,'(A)')                                            &
                0107           'MONIN_OBUKHOV_INIT: finished reading data.atm_gray'
                0108      CALL PRINT_MESSAGE( msgBuf, gcm_stdMsgUnit, gcm_SQZ_R, myThid )
                0109 !    Close the open data file
15388b052e Jean*0110 #ifdef SINGLE_DISK_IO
b2ea1d2979 Jean*0111      CLOSE(iUnit)
15388b052e Jean*0112 #else
                0113      CLOSE(iUnit,STATUS='DELETE')
                0114 #endif /* SINGLE_DISK_IO */
b2ea1d2979 Jean*0115
                0116 !      if (file_exist('input.nml')) then
                0117 !         unit = open_namelist_file ()
                0118 !         ierr=1; do while (ierr /= 0)
                0119 !            read  (unit, nml=monin_obukhov_nml, iostat=io, end=10)
                0120 !            ierr = check_nml_error(io,'monin_obukhov_nml')
                0121 !         enddo
                0122 !  10     call close_file (unit)
                0123 !      endif
                0124
                0125 !---------- output namelist to log-------------------------------------
                0126
                0127 !      call write_version_number(version, tagname)
                0128 !      if ( mpp_pe() == mpp_root_pe() ) write (stdlog(), nml=monin_obukhov_nml)
                0129
                0130 !----------------------------------------------------------------------
                0131
                0132 if(rich_crit.le.0.25)  call PRINT_ERROR( &
                0133         'MONIN_OBUKHOV_INIT in MONIN_OBUKHOV_MOD'// &
                0134          'rich_crit in monin_obukhov_mod must be > 0.25', myThid )
                0135 !if(rich_crit.le.0.25)  call error_mesg( &
                0136 !        'MONIN_OBUKHOV_INIT in MONIN_OBUKHOV_MOD', &
                0137 !        'rich_crit in monin_obukhov_mod must be > 0.25', FATAL)
                0138
                0139 if(drag_min.lt.0.0)  call PRINT_ERROR( &
                0140         'MONIN_OBUKHOV_INIT in MONIN_OBUKHOV_MOD'// &
                0141         'drag_min in monin_obukhov_mod must be >= 0.0', myThid )
                0142 !if(drag_min.lt.0.0)  call error_mesg( &
                0143 !        'MONIN_OBUKHOV_INIT in MONIN_OBUKHOV_MOD', &
                0144 !        'drag_min in monin_obukhov_mod must be >= 0.0', FATAL)
                0145
                0146 if(stable_option < 1 .or. stable_option > 2) call PRINT_ERROR( &
                0147         'MONIN_OBUKHOV_INIT in MONIN_OBUKHOV_MOD'// &
                0148         'the only allowable values of stable_option are 1 and 2', myThid )
                0149 !if(stable_option < 1 .or. stable_option > 2) call error_mesg( &
                0150 !        'MONIN_OBUKHOV_INIT in MONIN_OBUKHOV_MOD', &
                0151 !        'the only allowable values of stable_option are 1 and 2', FATAL)
                0152
                0153 if(stable_option == 2 .and. zeta_trans < 0) call PRINT_ERROR( &
                0154         'MONIN_OBUKHOV_INIT in MONIN_OBUKHOV_MOD'// &
                0155         'zeta_trans must be positive', myThid )
                0156 !if(stable_option == 2 .and. zeta_trans < 0) call error_mesg( &
                0157 !        'MONIN_OBUKHOV_INIT in MONIN_OBUKHOV_MOD', &
                0158 !        'zeta_trans must be positive', FATAL)
                0159
                0160 b_stab = 1.0/rich_crit
                0161 r_crit = 0.95*rich_crit  ! convergence can get slow if one is
                0162                          ! close to rich_crit
                0163
                0164 sqrt_drag_min = 0.0
                0165 if(drag_min.ne.0.0) sqrt_drag_min = sqrt(drag_min)
                0166
                0167 lambda     = 1.0 + (5.0 - b_stab)*zeta_trans   ! used only if stable_option = 2
                0168 rich_trans = zeta_trans/(1.0 + 5.0*zeta_trans) ! used only if stable_option = 2
                0169
                0170 init = .true.
                0171
                0172 !    _END_MASTER(myThid)
                0173 !    _BARRIER
                0174      ENDIF
                0175      CALL BARRIER(myThid)
                0176
                0177 return
                0178 end subroutine monin_obukhov_init
                0179
                0180 !=======================================================================
                0181
                0182 subroutine mo_drag_1d &
                0183          (pt, pt0, z, z0, zt, zq, speed, drag_m, drag_t, drag_q, &
                0184           u_star, b_star, myThid, avail)
                0185
                0186 real, intent(in)   , dimension(:) :: pt, pt0, z, z0, zt, zq, speed
                0187 real, intent(inout), dimension(:) :: drag_m, drag_t, drag_q, u_star, b_star
                0188 integer, intent(in)               :: myThid
                0189 logical, intent(in), optional, dimension(:) :: avail
                0190
                0191 real   , dimension(size(pt)) :: rich, fm, ft, fq, zz
                0192 logical, dimension(size(pt)) :: mask, mask_1, mask_2
                0193
                0194 real   , dimension(size(pt)) :: delta_b, us, bs, qs
                0195
                0196 if(.not.init) call monin_obukhov_init(myThid)
                0197
                0198 mask = .true.
                0199 if(present(avail)) mask = avail
                0200
                0201 where(mask)
                0202    delta_b = grav*(pt0 - pt)/pt0
                0203    rich    = - z*delta_b/(speed*speed + small)
                0204    zz      = max(z,z0,zt,zq)
                0205 else where
                0206    rich = 0.0
                0207 end where
                0208
                0209 if(neutral) then
                0210
                0211   where(mask)
                0212     fm   = log(zz/z0)
                0213     ft   = log(zz/zt)
                0214     fq   = log(zz/zq)
                0215     us   = vonkarm/fm
                0216     bs   = vonkarm/ft
                0217     qs   = vonkarm/fq
                0218     drag_m    = us*us
                0219     drag_t    = us*bs
                0220     drag_q    = us*qs
                0221     u_star = us*speed
                0222     b_star = bs*delta_b
                0223   end where
                0224
                0225 else
                0226
                0227   mask_1 = mask .and. rich <  r_crit
                0228   mask_2 = mask .and. rich >= r_crit
                0229
                0230   where(mask_2)
                0231     drag_m   = drag_min
                0232     drag_t   = drag_min
                0233     drag_q   = drag_min
                0234     us       = sqrt_drag_min
                0235     bs       = sqrt_drag_min
                0236     qs       = sqrt_drag_min
                0237     u_star   = us*speed
                0238     b_star   = bs*delta_b
                0239   end where
                0240   call solve_zeta (rich, zz, z0, zt, zq, fm, ft, fq, mask_1, myThid)
                0241
                0242   where (mask_1)
                0243     us   = max(vonkarm/fm, sqrt_drag_min)
                0244     bs   = max(vonkarm/ft, sqrt_drag_min)
                0245     qs   = max(vonkarm/fq, sqrt_drag_min)
                0246     drag_m   = us*us
                0247     drag_t   = us*bs
                0248     drag_q   = us*qs
                0249     u_star   = us*speed
                0250     b_star   = bs*delta_b
                0251   end where
                0252
                0253 end if
                0254
                0255 return
                0256 end subroutine mo_drag_1d
                0257
                0258 !=======================================================================
                0259
                0260 subroutine mo_profile_1d(zref, zref_t, z, z0, zt, zq, u_star, b_star, q_star, &
                0261                          del_m, del_t, del_q, myThid, avail)
                0262
                0263 real,    intent(in)                :: zref, zref_t
                0264 real,    intent(in) , dimension(:) :: z, z0, zt, zq, u_star, b_star, q_star
                0265 real,    intent(out), dimension(:) :: del_m, del_t, del_q
                0266 integer, intent(in)                :: myThid
                0267 logical, intent(in) , optional, dimension(:) :: avail
                0268
                0269 real, dimension(size(z)) :: zeta, zeta_0, zeta_t, zeta_q, zeta_ref, zeta_ref_t, &
                0270                             ln_z_z0, ln_z_zt, ln_z_zq, ln_z_zref, ln_z_zref_t,  &
                0271                             f_m_ref, f_m, f_t_ref, f_t, f_q_ref, f_q,           &
                0272                             mo_length_inv
                0273
                0274 logical, dimension(size(z)) :: mask
                0275
                0276 if(.not.init) call monin_obukhov_init(myThid)
                0277
                0278 mask = .true.
                0279 if(present(avail)) mask = avail
                0280
                0281 del_m = 0.0  ! zero output arrays
                0282 del_t = 0.0
                0283 del_q = 0.0
                0284
                0285 where(mask)
                0286   ln_z_z0     = log(z/z0)
                0287   ln_z_zt     = log(z/zt)
                0288   ln_z_zq     = log(z/zq)
                0289   ln_z_zref   = log(z/zref)
                0290   ln_z_zref_t = log(z/zref_t)
                0291 endwhere
                0292
                0293 if(neutral) then
                0294
                0295   where(mask)
                0296     del_m = 1.0 - ln_z_zref  /ln_z_z0
                0297     del_t = 1.0 - ln_z_zref_t/ln_z_zt
                0298     del_q = 1.0 - ln_z_zref_t/ln_z_zq
                0299   endwhere
                0300
                0301 else
                0302
                0303   where(mask .and. u_star > 0.0)
                0304     mo_length_inv = - vonkarm * b_star/(u_star*u_star)
                0305     zeta       = z     *mo_length_inv
                0306     zeta_0     = z0    *mo_length_inv
                0307     zeta_t     = zt    *mo_length_inv
                0308     zeta_q     = zq    *mo_length_inv
                0309     zeta_ref   = zref  *mo_length_inv
                0310     zeta_ref_t = zref_t*mo_length_inv
                0311   endwhere
                0312
                0313   call mo_integral_m(f_m,     zeta, zeta_0,   ln_z_z0,   mask, myThid)
                0314   call mo_integral_m(f_m_ref, zeta, zeta_ref, ln_z_zref, mask, myThid)
                0315
                0316   call mo_integral_tq(f_t, f_q, zeta, zeta_t, zeta_q, ln_z_zt, ln_z_zq, mask, myThid)
                0317   call mo_integral_tq(f_t_ref, f_q_ref, zeta, zeta_ref_t, zeta_ref_t, &
                0318                       ln_z_zref_t, ln_z_zref_t,  mask, myThid)
                0319
                0320   where(mask)
                0321     del_m = 1.0 - f_m_ref/f_m
                0322     del_t = 1.0 - f_t_ref/f_t
                0323     del_q = 1.0 - f_q_ref/f_q
                0324   endwhere
                0325
                0326 end if
                0327
                0328 return
                0329 end subroutine mo_profile_1d
                0330
                0331 !=======================================================================
                0332
                0333 subroutine stable_mix_3d(rich, mix, myThid)
                0334
                0335 real, intent(in) , dimension(:,:,:)  :: rich
                0336 real, intent(out), dimension(:,:,:)  :: mix
                0337 integer, intent(in)                  :: myThid
                0338
                0339 real, dimension(size(rich,1),size(rich,2),size(rich,3)) :: &
                0340       r, a, b, c, zeta, phi
                0341
                0342 mix = 0.0
                0343
                0344 if(stable_option == 1) then
                0345
                0346   where(rich > 0.0 .and. rich < rich_crit)
                0347     r = 1.0/rich
                0348     a = r - b_stab
                0349     b = r - (1.0 + 5.0)
                0350     c = - 1.0
                0351
                0352     zeta = (-b + sqrt(b*b - 4.0*a*c))/(2.0*a)
                0353     phi = 1.0 + b_stab*zeta + (5.0 - b_stab)*zeta/(1.0 + zeta)
                0354     mix = 1./(phi*phi)
                0355   end where
                0356
                0357 else if(stable_option == 2) then
                0358
                0359   where(rich > 0.0 .and. rich <= rich_trans)
                0360     mix = (1.0 - 5.0*rich)**2
                0361   end where
                0362   where(rich > rich_trans .and. rich < rich_crit)
                0363     mix = ((1.0 - b_stab*rich)/lambda)**2
                0364   end where
                0365
                0366 end if
                0367
                0368 return
                0369 end subroutine stable_mix_3d
                0370
                0371 !=======================================================================
                0372
                0373 subroutine mo_diff_2d_n(z, u_star, b_star, k_m, k_h, myThid)
                0374
                0375 real, intent(in),  dimension(:,:,:) :: z
                0376 real, intent(in),  dimension(:,:)   :: u_star, b_star
                0377 real, intent(out), dimension(:,:,:) :: k_m, k_h
                0378 integer, intent(in)                 :: myThid
                0379
                0380 real , dimension(size(z,1),size(z,2)) :: phi_m, phi_h, zeta, uss
                0381 integer :: j, k
                0382 INTEGER :: I
                0383 logical, dimension(size(z,1)) :: mask
                0384
                0385 if(.not.init) call monin_obukhov_init(myThid)
                0386
                0387 mask = .true.
                0388 uss = max(u_star, 1.e-10)
                0389
                0390 if(neutral) then
                0391   do k = 1, size(z,3)
                0392     k_m(:,:,k) = vonkarm *uss*z(:,:,k)
                0393     k_h(:,:,k) = k_m(:,:,k)
                0394   end do
                0395 else
                0396   do k = 1, size(z,3)
                0397     zeta = - vonkarm * b_star*z(:,:,k)/(uss*uss)
                0398
                0399     do j = 1, size(z,2)
                0400       call mo_derivative_m(phi_m(:,j), zeta(:,j), mask, myThid)
                0401       call mo_derivative_t(phi_h(:,j), zeta(:,j), mask, myThid)
                0402     enddo
                0403     k_m(:,:,k) = vonkarm * uss*z(:,:,k)/phi_m
                0404     k_h(:,:,k) = vonkarm * uss*z(:,:,k)/phi_h
                0405   end do
                0406 endif
                0407
                0408 return
                0409 end subroutine mo_diff_2d_n
                0410
                0411 !=======================================================================
                0412 ! The following routines are used by the public interfaces above
                0413 !=======================================================================
                0414
                0415 subroutine solve_zeta(rich, z, z0, zt, zq, f_m, f_t, f_q, mask, myThid)
                0416
                0417 real   , intent(in) , dimension(:) :: rich, z, z0, zt, zq
                0418 logical, intent(in) , dimension(:) :: mask
                0419 real   , intent(out), dimension(:) :: f_m, f_t, f_q
                0420 integer, intent(in)                :: myThid
                0421
                0422 real, parameter    :: error    = 1.e-04
                0423 real, parameter    :: zeta_min = 1.e-06
                0424 integer, parameter :: max_iter = 20
                0425
                0426 real    :: max_cor
                0427 integer :: iter
                0428
                0429 real, dimension(size(rich)) ::   &
                0430           d_rich, rich_1, correction, corr, z_z0, z_zt, z_zq, &
                0431           ln_z_z0, ln_z_zt, ln_z_zq, zeta,                    &
                0432           phi_m, phi_m_0, phi_t, phi_t_0, rzeta,              &
                0433           zeta_0, zeta_t, zeta_q, df_m, df_t
                0434
                0435 logical, dimension(size(rich)) :: mask_1
                0436
                0437 z_z0 = z/z0
                0438 z_zt = z/zt
                0439 z_zq = z/zq
                0440 ln_z_z0 = log(z_z0)
                0441 ln_z_zt = log(z_zt)
                0442 ln_z_zq = log(z_zq)
                0443
                0444 corr = 0.0
                0445 mask_1 = mask
                0446
                0447 ! initial guess
                0448
                0449 where(mask_1)
                0450   zeta = rich*ln_z_z0*ln_z_z0/ln_z_zt
                0451 else where
                0452   zeta = 0.0
                0453 end where
                0454
                0455 where (mask_1 .and. rich >= 0.0)
                0456   zeta = zeta/(1.0 - rich/rich_crit)
                0457 end where
                0458
                0459 iter_loop: do iter = 1, max_iter
                0460
                0461   where (mask_1 .and. abs(zeta).lt.zeta_min)
                0462     zeta = 0.0
                0463     f_m = ln_z_z0
                0464     f_t = ln_z_zt
                0465     f_q = ln_z_zq
                0466     mask_1 = .false.  ! don't do any more calculations at these pts
                0467   end where
                0468
                0469   where (mask_1)
                0470     rzeta  = 1.0/zeta
                0471     zeta_0 = zeta/z_z0
                0472     zeta_t = zeta/z_zt
                0473     zeta_q = zeta/z_zq
                0474   else where
                0475     zeta_0 = 0.0
                0476     zeta_t = 0.0
                0477     zeta_q = 0.0
                0478   end where
                0479
                0480   call mo_derivative_m(phi_m  , zeta  , mask_1, myThid)
                0481   call mo_derivative_m(phi_m_0, zeta_0, mask_1, myThid)
                0482   call mo_derivative_t(phi_t  , zeta  , mask_1, myThid)
                0483   call mo_derivative_t(phi_t_0, zeta_t, mask_1, myThid)
                0484
                0485   call mo_integral_m(f_m, zeta, zeta_0, ln_z_z0, mask_1, myThid)
                0486   call mo_integral_tq(f_t, f_q, zeta, zeta_t, zeta_q, ln_z_zt, ln_z_zq, mask_1, myThid)
                0487
                0488   where (mask_1)
                0489     df_m  = (phi_m - phi_m_0)*rzeta
                0490     df_t  = (phi_t - phi_t_0)*rzeta
                0491     rich_1 = zeta*f_t/(f_m*f_m)
                0492     d_rich = rich_1*( rzeta +  df_t/f_t - 2.0 *df_m/f_m)
                0493     correction = (rich - rich_1)/d_rich
                0494     corr = min(abs(correction),abs(correction/zeta))
                0495       ! the criterion corr < error seems to work ok, but is a bit arbitrary
                0496       !  when zeta is small the tolerance is reduced
                0497   end where
                0498
                0499   max_cor= maxval(corr)
                0500
                0501   if(max_cor > error) then
                0502     mask_1 = mask_1 .and. (corr > error)
                0503        ! change the mask so computation proceeds only on non-converged points
                0504     where(mask_1)
                0505       zeta = zeta + correction
                0506     end where
                0507     cycle iter_loop
                0508   else
                0509     return
                0510   end if
                0511
                0512 end do iter_loop
                0513
                0514 !call error_mesg ('solve_zeta in monin_obukhov_mod',  &
                0515 !                 'surface drag iteration did not converge', FATAL)
                0516
                0517 end subroutine solve_zeta
                0518
                0519 !=======================================================================
                0520
                0521 subroutine mo_derivative_m(phi_m, zeta, mask, myThid)
                0522
                0523 ! the differential similarity function for momentum
                0524
                0525 real    , intent(out),  dimension(:) :: phi_m
                0526 real    , intent(in),   dimension(:) :: zeta
                0527 logical , intent(in),   dimension(:) :: mask
                0528 integer, intent(in)                  :: myThid
                0529
                0530 logical, dimension(size(zeta)) :: stable, unstable
                0531 real   , dimension(size(zeta)) :: x
                0532
                0533 stable   = mask .and. zeta >= 0.0
                0534 unstable = mask .and. zeta <  0.0
                0535
                0536 where (unstable)
                0537   x     = (1 - 16.0*zeta  )**(-0.5)
                0538   phi_m = sqrt(x)  ! phi_m = (1 - 16.0*zeta)**(-0.25)
                0539 end where
                0540
                0541 if(stable_option == 1) then
                0542
                0543   where (stable)
                0544     phi_m = 1.0 + zeta  *(5.0 + b_stab*zeta)/(1.0 + zeta)
                0545   end where
                0546
                0547 else if(stable_option == 2) then
                0548
                0549   where (stable .and. zeta < zeta_trans)
                0550     phi_m = 1 + 5.0*zeta
                0551   end where
                0552   where (stable .and. zeta >= zeta_trans)
                0553     phi_m = lambda + b_stab*zeta
                0554   end where
                0555
                0556 endif
                0557
                0558 return
                0559 end subroutine mo_derivative_m
                0560
                0561 !=======================================================================
                0562
                0563 subroutine mo_derivative_t(phi_t, zeta, mask, myThid)
                0564
                0565 ! the differential similarity function for buoyancy and tracers
                0566
                0567 real    , intent(out),  dimension(:) :: phi_t
                0568 real    , intent(in),   dimension(:) :: zeta
                0569 logical , intent(in),   dimension(:) :: mask
                0570 integer, intent(in)                  :: myThid
                0571
                0572 logical, dimension(size(zeta)) :: stable, unstable
                0573
                0574 stable   = mask .and. zeta >= 0.0
                0575 unstable = mask .and. zeta <  0.0
                0576
                0577 where (unstable)
                0578   phi_t = (1 - 16.0*zeta)**(-0.5)
                0579 end where
                0580
                0581 if(stable_option == 1) then
                0582
                0583   where (stable)
                0584     phi_t = 1.0 + zeta*(5.0 + b_stab*zeta)/(1.0 + zeta)
                0585   end where
                0586
                0587 else if(stable_option == 2) then
                0588
                0589   where (stable .and. zeta < zeta_trans)
                0590     phi_t = 1 + 5.0*zeta
                0591   end where
                0592   where (stable .and. zeta >= zeta_trans)
                0593     phi_t = lambda + b_stab*zeta
                0594   end where
                0595
                0596 endif
                0597
                0598 return
                0599 end subroutine mo_derivative_t
                0600
                0601 !=======================================================================
                0602
                0603 subroutine mo_integral_tq (psi_t, psi_q, zeta, zeta_t, zeta_q, &
                0604                            ln_z_zt, ln_z_zq, mask, myThid)
                0605
                0606 ! the integral similarity function for moisture and tracers
                0607
                0608 real    , intent(out), dimension(:) :: psi_t, psi_q
                0609 real    , intent(in),  dimension(:) :: zeta, zeta_t, zeta_q, ln_z_zt, ln_z_zq
                0610 logical , intent(in),  dimension(:) :: mask
                0611 integer, intent(in)                 :: myThid
                0612
                0613 real, dimension(size(zeta)) :: x, x_t, x_q
                0614
                0615 logical, dimension(size(zeta)) :: stable, unstable, &
                0616                                   weakly_stable, strongly_stable
                0617
                0618 stable   = mask .and. zeta >= 0.0
                0619 unstable = mask .and. zeta <  0.0
                0620
                0621 where(unstable)
                0622
                0623   x     = sqrt(1 - 16.0*zeta)
                0624   x_t   = sqrt(1 - 16.0*zeta_t)
                0625   x_q   = sqrt(1 - 16.0*zeta_q)
                0626
                0627   psi_t = ln_z_zt - 2.0*log( (1.0 + x)/(1.0 + x_t) )
                0628   psi_q = ln_z_zq - 2.0*log( (1.0 + x)/(1.0 + x_q) )
                0629
                0630 end where
                0631
                0632 if( stable_option == 1) then
                0633
                0634   where (stable)
                0635
                0636     psi_t = ln_z_zt + (5.0 - b_stab)*log((1.0 + zeta)/(1.0 + zeta_t)) &
                0637        + b_stab*(zeta - zeta_t)
                0638     psi_q = ln_z_zq + (5.0 - b_stab)*log((1.0 + zeta)/(1.0 + zeta_q)) &
                0639        + b_stab*(zeta - zeta_q)
                0640
                0641   end where
                0642
                0643 else if (stable_option == 2) then
                0644
                0645   weakly_stable   = stable .and. zeta <= zeta_trans
                0646   strongly_stable = stable .and. zeta >  zeta_trans
                0647
                0648   where (weakly_stable)
                0649     psi_t = ln_z_zt + 5.0*(zeta - zeta_t)
                0650     psi_q = ln_z_zq + 5.0*(zeta - zeta_q)
                0651   end where
                0652
                0653   where(strongly_stable)
                0654     x = (lambda - 1.0)*log(zeta/zeta_trans) + b_stab*(zeta - zeta_trans)
                0655   endwhere
                0656
                0657   where (strongly_stable .and. zeta_t <= zeta_trans)
                0658     psi_t = ln_z_zt + x + 5.0*(zeta_trans - zeta_t)
                0659   end where
                0660   where (strongly_stable .and. zeta_t > zeta_trans)
                0661     psi_t = lambda*ln_z_zt + b_stab*(zeta  - zeta_t)
                0662   endwhere
                0663
                0664   where (strongly_stable .and. zeta_q <= zeta_trans)
                0665     psi_q = ln_z_zq + x + 5.0*(zeta_trans - zeta_q)
                0666   end where
                0667   where (strongly_stable .and. zeta_q > zeta_trans)
                0668     psi_q = lambda*ln_z_zq + b_stab*(zeta  - zeta_q)
                0669   endwhere
                0670
                0671 end if
                0672
                0673 return
                0674 end subroutine mo_integral_tq
                0675
                0676 !=======================================================================
                0677
                0678 subroutine mo_integral_m (psi_m, zeta, zeta_0, ln_z_z0, mask, myThid)
                0679
                0680 !  the integral similarity function for momentum
                0681
                0682 real    , intent(out), dimension(:) :: psi_m
                0683 real    , intent(in),  dimension(:) :: zeta, zeta_0, ln_z_z0
                0684 logical , intent(in),  dimension(:) :: mask
                0685 integer, intent(in)                 :: myThid
                0686
                0687 real, dimension(size(zeta)) :: x, x_0, x1, x1_0, num, denom, y
                0688
                0689 logical, dimension(size(zeta)) :: stable, unstable, &
                0690                                   weakly_stable, strongly_stable
                0691
                0692 stable   = mask .and. zeta >= 0.0
                0693 unstable = mask .and. zeta <  0.0
                0694
                0695 where(unstable)
                0696
                0697   x     = sqrt(1 - 16.0*zeta)
                0698   x_0   = sqrt(1 - 16.0*zeta_0)
                0699
                0700   x      = sqrt(x)
                0701   x_0    = sqrt(x_0)
                0702
                0703   x1     = 1.0 + x
                0704   x1_0   = 1.0 + x_0
                0705
                0706   num    = x1*x1*(1.0 + x*x)
                0707   denom  = x1_0*x1_0*(1.0 + x_0*x_0)
                0708   y      = atan(x) - atan(x_0)
                0709   psi_m  = ln_z_z0 - log(num/denom) + 2*y
                0710
                0711 end where
                0712
                0713 if( stable_option == 1) then
                0714
                0715   where (stable)
                0716     psi_m = ln_z_z0 + (5.0 - b_stab)*log((1.0 + zeta)/(1.0 + zeta_0)) &
                0717        + b_stab*(zeta - zeta_0)
                0718   end where
                0719
                0720 else if (stable_option == 2) then
                0721
                0722   weakly_stable   = stable .and. zeta <= zeta_trans
                0723   strongly_stable = stable .and. zeta >  zeta_trans
                0724
                0725   where (weakly_stable)
                0726     psi_m = ln_z_z0 + 5.0*(zeta - zeta_0)
                0727   end where
                0728
                0729   where(strongly_stable)
                0730     x = (lambda - 1.0)*log(zeta/zeta_trans) + b_stab*(zeta - zeta_trans)
                0731   endwhere
                0732
                0733   where (strongly_stable .and. zeta_0 <= zeta_trans)
                0734     psi_m = ln_z_z0 + x + 5.0*(zeta_trans - zeta_0)
                0735   end where
                0736   where (strongly_stable .and. zeta_0 > zeta_trans)
                0737     psi_m = lambda*ln_z_z0 + b_stab*(zeta  - zeta_0)
                0738   endwhere
                0739
                0740 end if
                0741
                0742 return
                0743 end subroutine mo_integral_m
                0744
                0745 !=======================================================================
                0746 ! The following routines allow the public interfaces to be used
                0747 ! with different dimensions of the input and output
                0748 !
                0749 !=======================================================================
                0750
                0751 subroutine mo_drag_2d &
                0752     (pt, pt0, z, z0, zt, zq, speed, drag_m, drag_t, drag_q, u_star, b_star, myThid)
                0753
                0754 real, intent(in)   , dimension(:,:) :: z, speed, pt, pt0, z0, zt, zq
                0755 real, intent(out)  , dimension(:,:) :: drag_m, drag_t, drag_q
                0756 real, intent(inout), dimension(:,:) :: u_star, b_star
                0757 integer, intent(in)                 :: myThid
                0758
                0759 integer :: j
                0760
                0761 do j = 1, size(pt,2)
                0762   call mo_drag_1d (pt(:,j), pt0(:,j), z(:,j), z0(:,j), zt(:,j), zq(:,j), &
                0763                    speed(:,j), drag_m(:,j), drag_t(:,j), drag_q(:,j), &
                0764                    u_star(:,j), b_star(:,j), myThid )
                0765 end do
                0766
                0767 return
                0768 end subroutine mo_drag_2d
                0769
                0770 !=======================================================================
                0771 subroutine mo_drag_0d &
                0772     (pt, pt0, z, z0, zt, zq, speed, drag_m, drag_t, drag_q, u_star, b_star, myThid)
                0773
                0774 real, intent(in)    :: z, speed, pt, pt0, z0, zt, zq
                0775 real, intent(out)   :: drag_m, drag_t, drag_q, u_star, b_star
                0776
                0777 real, dimension(1) :: pt_1, pt0_1, z_1, z0_1, zt_1, zq_1, speed_1, &
                0778                       drag_m_1, drag_t_1, drag_q_1, u_star_1, b_star_1
                0779 integer, intent(in):: myThid
                0780
                0781 pt_1   (1) = pt
                0782 pt0_1  (1) = pt0
                0783 z_1    (1) = z
                0784 z0_1   (1) = z0
                0785 zt_1   (1) = zt
                0786 zq_1   (1) = zq
                0787 speed_1(1) = speed
                0788
                0789 call mo_drag_1d (pt_1, pt0_1, z_1, z0_1, zt_1, zq_1, speed_1, &
                0790                  drag_m_1, drag_t_1, drag_q_1, u_star_1, b_star_1, myThid )
                0791
                0792 drag_m = drag_m_1(1)
                0793 drag_t = drag_t_1(1)
                0794 drag_q = drag_q_1(1)
                0795 u_star = u_star_1(1)
                0796 b_star = b_star_1(1)
                0797
                0798 return
                0799 end subroutine mo_drag_0d
                0800 !=======================================================================
                0801
                0802 subroutine mo_profile_2d(zref, zref_t, z, z0, zt, zq, u_star, b_star, q_star, &
                0803                          del_m, del_h, del_q, myThid)
                0804
                0805 real, intent(in)                  :: zref, zref_t
                0806 real, intent(in) , dimension(:,:) :: z, z0, zt, zq, u_star, b_star, q_star
                0807 real, intent(out), dimension(:,:) :: del_m, del_h, del_q
                0808 integer, intent(in)               :: myThid
                0809
                0810 integer :: j
                0811
                0812 do j = 1, size(z,2)
                0813   call mo_profile_1d (zref, zref_t, z(:,j), z0(:,j), zt(:,j),         &
                0814                       zq(:,j), u_star(:,j), b_star(:,j), q_star(:,j), &
                0815                       del_m(:,j), del_h (:,j), del_q (:,j), myThid )
                0816 enddo
                0817
                0818 return
                0819 end subroutine mo_profile_2d
                0820
                0821 !=======================================================================
                0822
                0823 subroutine mo_profile_0d(zref, zref_t, z, z0, zt, zq, u_star, b_star, q_star, &
                0824                          del_m, del_h, del_q, myThid)
                0825
                0826 real, intent(in)  :: zref, zref_t
                0827 real, intent(in)  :: z, z0, zt, zq, u_star, b_star, q_star
                0828 real, intent(out) :: del_m, del_h, del_q
                0829
                0830 real, dimension(1) :: z_1, z0_1, zt_1, zq_1, u_star_1, b_star_1, q_star_1, &
                0831                       del_m_1, del_h_1, del_q_1
                0832 integer, intent(in):: myThid
                0833
                0834 z_1     (1) = z
                0835 z0_1    (1) = z0
                0836 zt_1    (1) = zt
                0837 zq_1    (1) = zq
                0838 u_star_1(1) = u_star
                0839 b_star_1(1) = b_star
                0840 q_star_1(1) = q_star
                0841
                0842 call mo_profile_1d (zref, zref_t, z_1, z0_1, zt_1, zq_1, &
                0843                     u_star_1, b_star_1, q_star_1,        &
                0844                     del_m_1, del_h_1, del_q_1, myThid )
                0845
                0846 del_m = del_m_1(1)
                0847 del_h = del_h_1(1)
                0848 del_q = del_q_1(1)
                0849
                0850 return
                0851 end subroutine mo_profile_0d
                0852
                0853 !=======================================================================
                0854
                0855 subroutine mo_profile_1d_n(zref, z, z0, zt, zq, u_star, b_star, q_star, &
                0856                          del_m, del_t, del_q, myThid, avail)
                0857
                0858 real,    intent(in),  dimension(:)   :: zref
                0859 real,    intent(in) , dimension(:)   :: z, z0, zt, zq, u_star, b_star, q_star
                0860 real,    intent(out), dimension(:,:) :: del_m, del_t, del_q
                0861 integer, intent(in)                  :: myThid
                0862 logical, intent(in) , optional, dimension(:) :: avail
                0863
                0864 integer :: k
                0865
                0866 do k = 1, size(zref)
                0867   if(present(avail)) then
                0868     call mo_profile_1d (zref(k), zref(k), z, z0, zt, zq, &
                0869        u_star, b_star, q_star, del_m(:,k), del_t(:,k), del_q(:,k), myThid, avail)
                0870   else
                0871       call mo_profile_1d (zref(k), zref(k), z, z0, zt, zq, &
                0872        u_star, b_star, q_star, del_m(:,k), del_t(:,k), del_q(:,k), myThid )
                0873   endif
                0874 enddo
                0875
                0876 return
                0877 end subroutine mo_profile_1d_n
                0878
                0879 !=======================================================================
                0880
                0881 subroutine mo_profile_0d_n(zref, z, z0, zt, zq, u_star, b_star, q_star, &
                0882                          del_m, del_t, del_q, myThid)
                0883
                0884 real,    intent(in),  dimension(:) :: zref
                0885 real,    intent(in)                :: z, z0, zt, zq, u_star, b_star, q_star
                0886 real,    intent(out), dimension(:) :: del_m, del_t, del_q
                0887 integer, intent(in)                :: myThid
                0888
                0889 integer :: k
                0890
                0891 do k = 1, size(zref)
                0892   call mo_profile_0d (zref(k), zref(k), z, z0, zt, zq, &
                0893        u_star, b_star, q_star, del_m(k), del_t(k), del_q(k), myThid)
                0894 enddo
                0895
                0896 return
                0897 end subroutine mo_profile_0d_n
                0898
                0899 !=======================================================================
                0900
                0901 subroutine mo_profile_2d_n(zref, z, z0, zt, zq, u_star, b_star, q_star, &
                0902                          del_m, del_t, del_q, myThid)
                0903
                0904 real,    intent(in),  dimension(:)     :: zref
                0905 real,    intent(in),  dimension(:,:)   :: z, z0, zt, zq, u_star, b_star, q_star
                0906 real,    intent(out), dimension(:,:,:) :: del_m, del_t, del_q
                0907 integer, intent(in)                    :: myThid
                0908
                0909 integer :: k
                0910
                0911 do k = 1, size(zref)
                0912   call mo_profile_2d (zref(k), zref(k), z, z0, zt, zq, &
                0913        u_star, b_star, q_star, del_m(:,:,k), del_t(:,:,k), del_q(:,:,k), myThid)
                0914 enddo
                0915
                0916 return
                0917 end subroutine mo_profile_2d_n
                0918
                0919 !=======================================================================
                0920
                0921 subroutine mo_diff_2d_1(z, u_star, b_star, k_m, k_h, myThid)
                0922
                0923 real, intent(in),  dimension(:,:) :: z, u_star, b_star
                0924 real, intent(out), dimension(:,:) :: k_m, k_h
                0925
                0926 real   , dimension(size(z,1),size(z,2),1) :: z_n, k_m_n, k_h_n
                0927 integer, intent(in)                       :: myThid
                0928
                0929 z_n(:,:,1) = z
                0930
                0931 call mo_diff_2d_n(z_n, u_star, b_star, k_m_n, k_h_n, myThid)
                0932
                0933 k_m = k_m_n(:,:,1)
                0934 k_h = k_h_n(:,:,1)
                0935
                0936 return
                0937 end subroutine mo_diff_2d_1
                0938
                0939 !=======================================================================
                0940
                0941 subroutine mo_diff_1d_1(z, u_star, b_star, k_m, k_h, myThid)
                0942
                0943 real, intent(in),  dimension(:) :: z, u_star, b_star
                0944 real, intent(out), dimension(:) :: k_m, k_h
                0945
                0946 real, dimension(size(z),1,1) :: z_n, k_m_n, k_h_n
                0947 real, dimension(size(z),1)   :: u_star_n, b_star_n
                0948 integer, intent(in)          :: myThid
                0949
                0950 z_n   (:,1,1) = z
                0951 u_star_n(:,1) = u_star
                0952 b_star_n(:,1) = b_star
                0953
                0954 call mo_diff_2d_n(z_n, u_star_n, b_star_n, k_m_n, k_h_n, myThid)
                0955
                0956 k_m = k_m_n(:,1,1)
                0957 k_h = k_h_n(:,1,1)
                0958
                0959 return
                0960 end subroutine mo_diff_1d_1
                0961
                0962 !=======================================================================
                0963
                0964 subroutine mo_diff_1d_n(z, u_star, b_star, k_m, k_h, myThid)
                0965
                0966 real, intent(in),  dimension(:,:) :: z
                0967 real, intent(in),  dimension(:)   :: u_star, b_star
                0968 real, intent(out), dimension(:,:) :: k_m, k_h
                0969 integer, intent(in)               :: myThid
                0970
                0971 real, dimension(size(z,1),1)            :: u_star2, b_star2
                0972 real, dimension(size(z,1),size(z,2), 1) :: z2, k_m2, k_h2
                0973
                0974 z2   (:,:,1) = z
                0975 u_star2(:,1) = u_star
                0976 b_star2(:,1) = b_star
                0977
                0978 call mo_diff_2d_n(z2, u_star2, b_star2, k_m2, k_h2, myThid)
                0979
                0980 k_m = k_m2(:,:,1)
                0981 k_h = k_h2(:,:,1)
                0982
                0983 return
                0984 end subroutine mo_diff_1d_n
                0985
                0986 !=======================================================================
                0987
                0988 subroutine mo_diff_0d_1(z, u_star, b_star, k_m, k_h, myThid)
                0989
                0990 real, intent(in)   :: z, u_star, b_star
                0991 real, intent(out)  :: k_m, k_h
                0992 integer, intent(in):: myThid
                0993
                0994 real, dimension(1,1,1) :: z_n, k_m_n, k_h_n
                0995 real, dimension(1,1)   :: u_star_n, b_star_n
                0996
                0997 z_n   (1,1,1) = z
                0998 u_star_n(1,1) = u_star
                0999 b_star_n(1,1) = b_star
                1000
                1001 call mo_diff_2d_n(z_n, u_star_n, b_star_n, k_m_n, k_h_n, myThid)
                1002
                1003 k_m = k_m_n(1,1,1)
                1004 k_h = k_h_n(1,1,1)
                1005
                1006 return
                1007 end subroutine mo_diff_0d_1
                1008
                1009 !=======================================================================
                1010
                1011 subroutine mo_diff_0d_n(z, u_star, b_star, k_m, k_h, myThid)
                1012
                1013 real, intent(in),  dimension(:) :: z
                1014 real, intent(in)                :: u_star, b_star
                1015 real, intent(out), dimension(:) :: k_m, k_h
                1016 integer, intent(in)             :: myThid
                1017
                1018 real, dimension(1,1)            :: u_star2, b_star2
                1019 real, dimension(size(z,1),1, 1) :: z2, k_m2, k_h2
                1020
                1021 z2   (:,1,1) = z
                1022 u_star2(1,1) = u_star
                1023 b_star2(1,1) = b_star
                1024
                1025 call mo_diff_2d_n(z2, u_star2, b_star2, k_m2, k_h2, myThid)
                1026
                1027 k_m = k_m2(:,1,1)
                1028 k_h = k_h2(:,1,1)
                1029
                1030 return
                1031 end subroutine mo_diff_0d_n
                1032
                1033 !=======================================================================
                1034
                1035 subroutine stable_mix_2d(rich, mix, myThid)
                1036
                1037 real, intent(in) , dimension(:,:)  :: rich
                1038 real, intent(out), dimension(:,:)  :: mix
                1039 integer, intent(in)                :: myThid
                1040
                1041 real, dimension(size(rich,1),size(rich,2),1) :: rich_3d, mix_3d
                1042
                1043 rich_3d(:,:,1) = rich
                1044
                1045 call stable_mix_3d(rich_3d, mix_3d, myThid)
                1046
                1047 mix = mix_3d(:,:,1)
                1048
                1049 return
                1050 end subroutine stable_mix_2d
                1051
                1052 !=======================================================================
                1053
                1054 subroutine stable_mix_1d(rich, mix, myThid)
                1055
                1056 real, intent(in) , dimension(:)  :: rich
                1057 real, intent(out), dimension(:)  :: mix
                1058 integer, intent(in)              :: myThid
                1059
                1060 real, dimension(size(rich),1,1) :: rich_3d, mix_3d
                1061
                1062 rich_3d(:,1,1) = rich
                1063
                1064 call stable_mix_3d(rich_3d, mix_3d, myThid)
                1065
                1066 mix = mix_3d(:,1,1)
                1067
                1068 return
                1069 end subroutine stable_mix_1d
                1070
                1071 !=======================================================================
                1072
                1073 subroutine stable_mix_0d(rich, mix, myThid)
                1074
                1075 real, intent(in)    :: rich
                1076 real, intent(out)   :: mix
                1077 integer, intent(in) :: myThid
                1078
                1079 real, dimension(1,1,1) :: rich_3d, mix_3d
                1080
                1081 rich_3d(1,1,1) = rich
                1082
                1083 call stable_mix_3d(rich_3d, mix_3d, myThid)
                1084
                1085 mix = mix_3d(1,1,1)
                1086
                1087 return
                1088 end subroutine stable_mix_0d
                1089 !=======================================================================
                1090
                1091 end module monin_obukhov_mod
                1092