File indexing completed on 2018-03-02 18:37:44 UTC
view on githubraw file Latest commit add29e06 on 2018-01-31 20:35:05 UTC
b2ea1d2979 Jean*0001 module monin_obukhov_mod
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
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
0018
0019
0020
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
0053
15388b052e Jean*0054
b2ea1d2979 Jean*0055
0056
0057
0058
0059
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
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
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
0089
0090
0091
0092
0093
0094
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
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
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
0117
0118
0119
0120
0121
0122
0123
0124
0125
0126
0127
0128
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
0136
0137
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
0143
0144
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
0150
0151
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
0157
0158
0159
0160 b_stab = 1.0/rich_crit
0161 r_crit = 0.95*rich_crit
0162
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
0168 rich_trans = zeta_trans/(1.0 + 5.0*zeta_trans)
0169
0170 init = .true.
0171
0172
0173
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
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
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
0445 mask_1 = mask
0446
0447
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.
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
0496
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
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
0515
0516
0517 end subroutine solve_zeta
0518
0519
0520
0521 subroutine mo_derivative_m(phi_m, zeta, mask, myThid)
0522
0523
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)
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
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
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
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
0747
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
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
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
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
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
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
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
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
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
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
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