FMS 2025.01-dev
Flexible Modeling System
Loading...
Searching...
No Matches
monin_obukhov_inter.inc
1!***********************************************************************
2!* GNU Lesser General Public License
3!*
4!* This file is part of the GFDL Flexible Modeling System (FMS).
5!*
6!* FMS is free software: you can redistribute it and/or modify it under
7!* the terms of the GNU Lesser General Public License as published by
8!* the Free Software Foundation, either version 3 of the License, or (at
9!* your option) any later version.
10!*
11!* FMS is distributed in the hope that it will be useful, but WITHOUT
12!* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
13!* FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
14!* for more details.
15!*
16!* You should have received a copy of the GNU Lesser General Public
17!* License along with FMS. If not, see <http://www.gnu.org/licenses/>.
18!***********************************************************************
19
20!> @addtogroup monin_obukhov_inter
21!> @{
22
23
24pure subroutine monin_obukhov_diff_(vonkarm, &
25 & ustar_min, &
26 & neutral, stable_option,new_mo_option,rich_crit, zeta_trans, &
27 & ni, nj, nk, z, u_star, b_star, k_m, k_h, ier)
28
29 real(kind=fms_mo_kind_), intent(in) :: vonkarm
30 real(kind=fms_mo_kind_), intent(in) :: ustar_min !< = 1.0E-10
31 logical, intent(in) :: neutral
32 integer, intent(in) :: stable_option
33 logical, intent(in) :: new_mo_option !miz
34 real(kind=fms_mo_kind_), intent(in) :: rich_crit, zeta_trans
35 integer, intent(in) :: ni, nj, nk
36 real(kind=fms_mo_kind_), intent(in), dimension(ni, nj, nk) :: z
37 real(kind=fms_mo_kind_), intent(in), dimension(ni, nj) :: u_star, b_star
38 real(kind=fms_mo_kind_), intent(out), dimension(ni, nj, nk) :: k_m, k_h
39 integer, intent(out) :: ier
40
41 real(kind=fms_mo_kind_), dimension(ni, nj) :: phi_m, phi_h, zeta, uss
42 integer :: j, k
43 logical, dimension(ni) :: mask
44
45 ier = 0
46
47 mask = .true.
48 uss = max(u_star, ustar_min)
49
50 if(neutral) then
51 do k = 1, size(z,3)
52 k_m(:,:,k) = vonkarm *uss*z(:,:,k)
53 k_h(:,:,k) = k_m(:,:,k)
54 end do
55 else
56 do k = 1, size(z,3)
57 zeta = - vonkarm * b_star*z(:,:,k)/(uss*uss)
58 do j = 1, size(z,2)
59 call monin_obukhov_derivative_m(stable_option, rich_crit, zeta_trans, &
60 & ni, phi_m(:,j), zeta(:,j), mask, ier)
61 call monin_obukhov_derivative_t(stable_option, new_mo_option,rich_crit, zeta_trans, &
62 & ni, phi_h(:,j), zeta(:,j), mask, ier)
63 enddo
64 k_m(:,:,k) = vonkarm * uss*z(:,:,k)/phi_m
65 k_h(:,:,k) = vonkarm * uss*z(:,:,k)/phi_h
66 end do
67 endif
68
69end subroutine monin_obukhov_diff_
70
71
72pure subroutine monin_obukhov_drag_1d_(grav, vonkarm, &
73 & error, zeta_min, max_iter, small, &
74 & neutral, stable_option, new_mo_option, rich_crit, zeta_trans,&
75 & drag_min_heat, drag_min_moist, drag_min_mom, &
76 & n, pt, pt0, z, z0, zt, zq, speed, drag_m, drag_t, &
77 & drag_q, u_star, b_star, lavail, avail, ier)
78
79 real(kind=fms_mo_kind_), intent(in) :: grav
80 real(kind=fms_mo_kind_), intent(in) :: vonkarm
81 real(kind=fms_mo_kind_), intent(in) :: error !< = 1.0E-04
82 real(kind=fms_mo_kind_), intent(in) :: zeta_min !< = 1.0E-06
83 integer, intent(in) :: max_iter !< = 20
84 real(kind=fms_mo_kind_), intent(in) :: small !< = 1.0E-04
85 logical, intent(in) :: neutral
86 integer, intent(in) :: stable_option
87 logical, intent(in) :: new_mo_option
88 real(kind=fms_mo_kind_), intent(in) :: rich_crit, zeta_trans
89 real(kind=fms_mo_kind_), intent(in) :: drag_min_heat, drag_min_moist, drag_min_mom
90 integer, intent(in) :: n
91 real(kind=fms_mo_kind_), intent(in), dimension(n) :: pt, pt0, z, z0, zt, zq, speed
92 real(kind=fms_mo_kind_), intent(inout), dimension(n) :: drag_m, drag_t, drag_q, u_star, b_star
93 logical, intent(in) :: lavail !< whether to use provided mask or not
94 logical, intent(in), dimension(n) :: avail !< provided mask
95 integer, intent(out) :: ier
96
97 real(kind=fms_mo_kind_), dimension(n) :: rich, fm, ft, fq, zz
98 logical, dimension(n) :: mask, mask_1, mask_2
99 real(kind=fms_mo_kind_), dimension(n) :: delta_b !!, us, bs, qs
100 real(kind=fms_mo_kind_) :: r_crit, sqrt_drag_min_heat
101 real(kind=fms_mo_kind_) :: sqrt_drag_min_moist, sqrt_drag_min_mom
102 real(kind=fms_mo_kind_) :: us, bs, qs
103 integer :: i
104 integer, parameter :: lkind = fms_mo_kind_
105
106 ier = 0
107 r_crit = 0.95_lkind*rich_crit ! convergence can get slow if one is
108 ! close to rich_crit
109 sqrt_drag_min_heat = 0.0_lkind
110 if(drag_min_heat.ne.0.0_lkind) sqrt_drag_min_heat = sqrt(drag_min_heat)
111 sqrt_drag_min_moist = 0.0_lkind
112 if(drag_min_moist.ne.0.0_lkind) sqrt_drag_min_moist = sqrt(drag_min_moist)
113 sqrt_drag_min_mom = 0.0_lkind
114 if(drag_min_mom.ne.0.0_lkind) sqrt_drag_min_mom = sqrt(drag_min_mom)
115
116 mask = .true.
117 if(lavail) mask = avail
118
119 where(mask)
120 delta_b = grav*(pt0 - pt)/pt0
121 rich = - z*delta_b/(speed*speed + small)
122 zz = max(z,z0,zt,zq)
123 elsewhere
124 rich = 0.0_lkind
125 end where
126
127 if(neutral) then
128
129 do i = 1, n
130 if(mask(i)) then
131 fm(i) = log(zz(i)/z0(i))
132 ft(i) = log(zz(i)/zt(i))
133 fq(i) = log(zz(i)/zq(i))
134 us = vonkarm/fm(i)
135 bs = vonkarm/ft(i)
136 qs = vonkarm/fq(i)
137 drag_m(i) = us*us
138 drag_t(i) = us*bs
139 drag_q(i) = us*qs
140 u_star(i) = us*speed(i)
141 b_star(i) = bs*delta_b(i)
142 end if
143 enddo
144
145 else
146
147 mask_1 = mask .and. rich < r_crit
148 mask_2 = mask .and. rich >= r_crit
149
150 do i = 1, n
151 if(mask_2(i)) then
152 drag_m(i) = drag_min_mom
153 drag_t(i) = drag_min_heat
154 drag_q(i) = drag_min_moist
155 us = sqrt_drag_min_mom
156 bs = sqrt_drag_min_heat
157 u_star(i) = us*speed(i)
158 b_star(i) = bs*delta_b(i)
159 end if
160 enddo
161
162 call monin_obukhov_solve_zeta (error, zeta_min, max_iter, small, &
163 & stable_option, new_mo_option, rich_crit, zeta_trans, &
164 & n, rich, zz, z0, zt, zq, fm, ft, fq, mask_1, ier)
165
166 do i = 1, n
167 if(mask_1(i)) then
168 us = max(vonkarm/fm(i), sqrt_drag_min_mom)
169 bs = max(vonkarm/ft(i), sqrt_drag_min_heat)
170 qs = max(vonkarm/fq(i), sqrt_drag_min_moist)
171 drag_m(i) = us*us
172 drag_t(i) = us*bs
173 drag_q(i) = us*qs
174 u_star(i) = us*speed(i)
175 b_star(i) = bs*delta_b(i)
176 endif
177 enddo
178
179 end if
180
181end subroutine monin_obukhov_drag_1d_
182
183
184pure subroutine monin_obukhov_solve_zeta_(error, zeta_min, max_iter, small, &
185 & stable_option, new_mo_option, rich_crit, zeta_trans, & !miz
186 & n, rich, z, z0, zt, zq, f_m, f_t, f_q, mask, ier)
187
188 real(kind=fms_mo_kind_), intent(in) :: error !< = 1.0E-04
189 real(kind=fms_mo_kind_), intent(in) :: zeta_min !< = 1.0E-06
190 integer, intent(in) :: max_iter !< = 20
191 real(kind=fms_mo_kind_), intent(in) :: small !< = 1.0E-04
192 integer, intent(in) :: stable_option
193 logical, intent(in) :: new_mo_option
194 real(kind=fms_mo_kind_), intent(in) :: rich_crit, zeta_trans
195 integer, intent(in) :: n
196 real(kind=fms_mo_kind_), intent(in), dimension(n) :: rich, z, z0, zt, zq
197 logical, intent(in), dimension(n) :: mask
198 real(kind=fms_mo_kind_), intent(out), dimension(n) :: f_m, f_t, f_q
199 integer, intent(out) :: ier
200
201 real(kind=fms_mo_kind_) :: max_cor
202 integer :: iter
203 real(kind=fms_mo_kind_), dimension(n) :: &
204 d_rich, rich_1, correction, corr, z_z0, z_zt, z_zq, &
205 ln_z_z0, ln_z_zt, ln_z_zq, zeta, &
206 phi_m, phi_m_0, phi_t, phi_t_0, rzeta, &
207 zeta_0, zeta_t, zeta_q, df_m, df_t
208 logical, dimension(n) :: mask_1
209 integer, parameter :: lkind = fms_mo_kind_
210
211 ier = 0
212
213 z_z0 = z/z0
214 z_zt = z/zt
215 z_zq = z/zq
216 ln_z_z0 = log(z_z0)
217 ln_z_zt = log(z_zt)
218 ln_z_zq = log(z_zq)
219
220 corr = 0.0_lkind
221 mask_1 = mask
222
223 ! initial guess
224
225 zeta = 0.0_lkind
226 where(mask_1)
227 zeta = rich*ln_z_z0*ln_z_z0/ln_z_zt
228 end where
229
230 where (mask_1 .and. rich >= 0.0_lkind)
231 zeta = zeta/(1.0_lkind - rich/rich_crit)
232 end where
233
234 iter_loop: do iter = 1, max_iter
235
236 where (mask_1 .and. abs(zeta).lt.zeta_min)
237 zeta = 0.0_lkind
238 f_m = ln_z_z0
239 f_t = ln_z_zt
240 f_q = ln_z_zq
241 mask_1 = .false. ! don't do any more calculations at these pts
242 end where
243
244
245 zeta_0 = 0.0_lkind
246 zeta_t = 0.0_lkind
247 zeta_q = 0.0_lkind
248 where (mask_1)
249 rzeta = 1.0_lkind/zeta
250 zeta_0 = zeta/z_z0
251 zeta_t = zeta/z_zt
252 zeta_q = zeta/z_zq
253 end where
254
255 call monin_obukhov_derivative_m(stable_option, rich_crit, zeta_trans, &
256 & n, phi_m , zeta , mask_1, ier)
257 call monin_obukhov_derivative_m(stable_option, rich_crit, zeta_trans, &
258 & n, phi_m_0, zeta_0, mask_1, ier)
259 call monin_obukhov_derivative_t(stable_option, new_mo_option,rich_crit, zeta_trans, &
260 & n, phi_t , zeta , mask_1, ier)
261 call monin_obukhov_derivative_t(stable_option, new_mo_option,rich_crit, zeta_trans, &
262 & n, phi_t_0, zeta_t, mask_1, ier)
263
264 call monin_obukhov_integral_m(stable_option, rich_crit, zeta_trans, &
265 & n, f_m, zeta, zeta_0, ln_z_z0, mask_1, ier)
266 call monin_obukhov_integral_tq(stable_option, new_mo_option, rich_crit, zeta_trans, &
267 & n, f_t, f_q, zeta, zeta_t, zeta_q, ln_z_zt, ln_z_zq, mask_1, ier)
268
269 where (mask_1)
270 df_m = (phi_m - phi_m_0)*rzeta
271 df_t = (phi_t - phi_t_0)*rzeta
272 rich_1 = zeta*f_t/(f_m*f_m)
273 d_rich = rich_1*( rzeta + df_t/f_t - 2.0_lkind *df_m/f_m)
274 correction = (rich - rich_1)/d_rich
275 corr = min(abs(correction),abs(correction/zeta))
276 ! the criterion corr < error seems to work ok, but is a bit arbitrary
277 ! when zeta is small the tolerance is reduced
278 end where
279
280 max_cor = maxval(corr)
281
282 if(max_cor > error) then
283 mask_1 = mask_1 .and. (corr > error)
284 ! change the mask so computation proceeds only on non-converged points
285 where(mask_1)
286 zeta = zeta + correction
287 end where
288 cycle iter_loop
289 else
290 return
291 end if
292
293 end do iter_loop
294
295 ier = 1 ! surface drag iteration did not converge
296
297end subroutine monin_obukhov_solve_zeta_
298
299
300!> The differential similarity function for buoyancy and tracers
301! seems to be the same as monin_obukhov_derivative_m?
302pure subroutine monin_obukhov_derivative_t_(stable_option,new_mo_option,rich_crit, zeta_trans, &
303 & n, phi_t, zeta, mask, ier)
304
305 integer, intent(in) :: stable_option
306 logical, intent(in) :: new_mo_option !miz
307 real(kind=fms_mo_kind_), intent(in) :: rich_crit, zeta_trans
308 integer, intent(in) :: n
309 real(kind=fms_mo_kind_), intent(out), dimension(n) :: phi_t
310 real(kind=fms_mo_kind_), intent(in), dimension(n) :: zeta
311 logical, intent(in), dimension(n) :: mask
312 integer, intent(out) :: ier
313
314 logical, dimension(n) :: stable, unstable
315 real(kind=fms_mo_kind_) :: b_stab, lambda
316 integer, parameter :: lkind = fms_mo_kind_
317
318 ier = 0
319 b_stab = 1.0_lkind/rich_crit
320
321 stable = mask .and. zeta >= 0.0_lkind
322 unstable = mask .and. zeta < 0.0_lkind
323
324!miz: modified to include new monin-obukhov option
325 if (new_mo_option) then
326 where (unstable)
327 phi_t = (1.0_lkind - 16.0_lkind*zeta)**(-1.0_lkind/3.0_lkind)
328 end where
329 else
330 where (unstable)
331 phi_t = (1.0_lkind - 16.0_lkind*zeta)**(-0.5_lkind)
332 end where
333 end if
334!miz
335
336 if(stable_option == 1) then
337
338 where (stable)
339 phi_t = 1.0_lkind + zeta*(5.0_lkind + b_stab*zeta)/(1.0_lkind + zeta)
340 end where
341
342 else if(stable_option == 2) then
343
344 lambda = 1.0_lkind + (5.0_lkind - b_stab)*zeta_trans
345
346 where (stable .and. zeta < zeta_trans)
347 phi_t = 1.0_lkind + 5.0_lkind*zeta
348 end where
349 where (stable .and. zeta >= zeta_trans)
350 phi_t = lambda + b_stab*zeta
351 end where
352
353 endif
354
355end subroutine monin_obukhov_derivative_t_
356
357
358! the differential similarity function for momentum
359pure subroutine monin_obukhov_derivative_m_(stable_option, rich_crit, zeta_trans, &
360 & n, phi_m, zeta, mask, ier)
361
362 integer, intent(in) :: stable_option
363 real(kind=fms_mo_kind_), intent(in) :: rich_crit, zeta_trans
364 integer, intent(in) :: n
365 real(kind=fms_mo_kind_), intent(out), dimension(n) :: phi_m
366 real(kind=fms_mo_kind_), intent(in), dimension(n) :: zeta
367 logical, intent(in), dimension(n) :: mask
368 integer, intent(out) :: ier
369
370 logical, dimension(n) :: stable, unstable
371 real(kind=fms_mo_kind_), dimension(n) :: x
372 real(kind=fms_mo_kind_) :: b_stab, lambda
373 integer, parameter :: lkind = fms_mo_kind_
374
375 ier = 0
376 b_stab = 1.0_lkind/rich_crit
377
378 stable = mask .and. zeta >= 0.0_lkind
379 unstable = mask .and. zeta < 0.0_lkind
380
381 where (unstable)
382 x = (1.0_lkind - 16.0_lkind*zeta )**(-0.5_lkind)
383 phi_m = sqrt(x) ! phi_m = (1 - 16.0*zeta)**(-0.25)
384 end where
385
386 if(stable_option == 1) then
387
388 where (stable)
389 phi_m = 1.0_lkind + zeta *(5.0_lkind + b_stab*zeta)/(1.0_lkind + zeta)
390 end where
391
392 else if(stable_option == 2) then
393
394 lambda = 1.0_lkind + (5.0_lkind - b_stab)*zeta_trans
395
396 where (stable .and. zeta < zeta_trans)
397 phi_m = 1.0_lkind + 5.0_lkind*zeta
398 end where
399 where (stable .and. zeta >= zeta_trans)
400 phi_m = lambda + b_stab*zeta
401 end where
402
403 endif
404
405end subroutine monin_obukhov_derivative_m_
406
407
408pure subroutine monin_obukhov_profile_1d_(vonkarm, &
409 & neutral, stable_option, new_mo_option, rich_crit, zeta_trans, &
410 & n, zref, zref_t, z, z0, zt, zq, u_star, b_star, q_star, &
411 & del_m, del_t, del_q, lavail, avail, ier)
412
413 real(kind=fms_mo_kind_), intent(in) :: vonkarm
414 logical, intent(in) :: neutral
415 integer, intent(in) :: stable_option
416 logical, intent(in) :: new_mo_option
417 real(kind=fms_mo_kind_), intent(in) :: rich_crit, zeta_trans
418 integer, intent(in) :: n
419 real(kind=fms_mo_kind_), intent(in) :: zref, zref_t
420 real(kind=fms_mo_kind_), intent(in), dimension(n) :: z, z0, zt, zq, u_star, b_star, q_star
421 real(kind=fms_mo_kind_), intent(out), dimension(n) :: del_m, del_t, del_q
422 logical, intent(in) :: lavail !< whether to use provided mask or not
423 logical, intent(in), dimension(n) :: avail !< provided mask
424 integer, intent(out) :: ier
425
426 real(kind=fms_mo_kind_), dimension(n) :: zeta, zeta_0, zeta_t, zeta_q, zeta_ref, zeta_ref_t, &
427 ln_z_z0, ln_z_zt, ln_z_zq, ln_z_zref, ln_z_zref_t, &
428 f_m_ref, f_m, f_t_ref, f_t, f_q_ref, f_q, &
429 mo_length_inv
430 logical, dimension(n) :: mask
431 integer, parameter :: lkind = fms_mo_kind_
432
433 ier = 0
434
435 mask = .true.
436 if(lavail) mask = avail
437
438 del_m = 0.0_lkind ! zero output arrays
439 del_t = 0.0_lkind
440 del_q = 0.0_lkind
441
442 where(mask)
443 ln_z_z0 = log(z/z0)
444 ln_z_zt = log(z/zt)
445 ln_z_zq = log(z/zq)
446 ln_z_zref = log(z/zref)
447 ln_z_zref_t = log(z/zref_t)
448 endwhere
449
450 if(neutral) then
451
452 where(mask)
453 del_m = 1.0_lkind - ln_z_zref /ln_z_z0
454 del_t = 1.0_lkind - ln_z_zref_t/ln_z_zt
455 del_q = 1.0_lkind - ln_z_zref_t/ln_z_zq
456 endwhere
457
458 else
459
460 where(mask .and. u_star > 0.0_lkind)
461 mo_length_inv = - vonkarm * b_star/(u_star*u_star)
462 zeta = z *mo_length_inv
463 zeta_0 = z0 *mo_length_inv
464 zeta_t = zt *mo_length_inv
465 zeta_q = zq *mo_length_inv
466 zeta_ref = zref *mo_length_inv
467 zeta_ref_t = zref_t*mo_length_inv
468 endwhere
469
470 call monin_obukhov_integral_m(stable_option, rich_crit, zeta_trans, &
471 & n, f_m, zeta, zeta_0, ln_z_z0, mask, ier)
472 call monin_obukhov_integral_m(stable_option, rich_crit, zeta_trans, &
473 & n, f_m_ref, zeta, zeta_ref, ln_z_zref, mask, ier)
474
475 call monin_obukhov_integral_tq(stable_option, new_mo_option, rich_crit, zeta_trans, &
476 & n, f_t, f_q, zeta, zeta_t, zeta_q, ln_z_zt, ln_z_zq, mask, ier)
477 call monin_obukhov_integral_tq(stable_option, new_mo_option, rich_crit, zeta_trans, &
478 & n, f_t_ref, f_q_ref, zeta, zeta_ref_t, zeta_ref_t, ln_z_zref_t, ln_z_zref_t, mask, ier)
479
480 where(mask)
481 del_m = 1.0_lkind - f_m_ref/f_m
482 del_t = 1.0_lkind - f_t_ref/f_t
483 del_q = 1.0_lkind - f_q_ref/f_q
484 endwhere
485
486 end if
487
488
489end subroutine monin_obukhov_profile_1d_
490
491
492!> The integral similarity function for momentum
493pure subroutine monin_obukhov_integral_m_(stable_option, rich_crit, zeta_trans, &
494 & n, psi_m, zeta, zeta_0, ln_z_z0, mask, ier)
495
496 integer, intent(in) :: stable_option
497 real(kind=fms_mo_kind_), intent(in) :: rich_crit, zeta_trans
498 integer, intent(in) :: n
499 real(kind=fms_mo_kind_), intent(inout), dimension(n) :: psi_m
500 real(kind=fms_mo_kind_), intent(in) , dimension(n) :: zeta, zeta_0, ln_z_z0
501 logical, intent(in) , dimension(n) :: mask
502 integer, intent(out) :: ier
503
504 real(kind=fms_mo_kind_) :: b_stab, lambda
505 real(kind=fms_mo_kind_), dimension(n) :: x, x_0, x1, x1_0, num, denom, y
506 logical, dimension(n) :: stable, unstable, &
507 weakly_stable, strongly_stable
508 integer, parameter :: lkind = fms_mo_kind_
509
510 ier = 0
511
512 b_stab = 1.0_lkind/rich_crit
513
514 stable = mask .and. zeta >= 0.0_lkind
515 unstable = mask .and. zeta < 0.0_lkind
516
517 where(unstable)
518
519 x = sqrt(1.0_lkind - 16.0_lkind*zeta)
520 x_0 = sqrt(1.0_lkind - 16.0_lkind*zeta_0)
521
522 x = sqrt(x)
523 x_0 = sqrt(x_0)
524
525 x1 = 1.0_lkind + x
526 x1_0 = 1.0_lkind + x_0
527
528 num = x1*x1*(1.0_lkind + x*x)
529 denom = x1_0*x1_0*(1.0_lkind + x_0*x_0)
530 y = atan(x) - atan(x_0)
531 psi_m = ln_z_z0 - log(num/denom) + 2.0_lkind*y
532
533 end where
534
535 if( stable_option == 1) then
536
537 where (stable)
538 psi_m = ln_z_z0 + (5.0_lkind - b_stab)*log((1.0_lkind + zeta)/(1.0_lkind + zeta_0)) &
539 + b_stab*(zeta - zeta_0)
540 end where
541
542 else if (stable_option == 2) then
543
544 lambda = 1.0_lkind + (5.0_lkind - b_stab)*zeta_trans
545
546 weakly_stable = stable .and. zeta <= zeta_trans
547 strongly_stable = stable .and. zeta > zeta_trans
548
549 where (weakly_stable)
550 psi_m = ln_z_z0 + 5.0_lkind*(zeta - zeta_0)
551 end where
552
553 where(strongly_stable)
554 x = (lambda - 1.0_lkind)*log(zeta/zeta_trans) + b_stab*(zeta - zeta_trans)
555 endwhere
556
557 where (strongly_stable .and. zeta_0 <= zeta_trans)
558 psi_m = ln_z_z0 + x + 5.0_lkind*(zeta_trans - zeta_0)
559 end where
560 where (strongly_stable .and. zeta_0 > zeta_trans)
561 psi_m = lambda*ln_z_z0 + b_stab*(zeta - zeta_0)
562 endwhere
563
564 end if
565
566end subroutine monin_obukhov_integral_m_
567
568
569!> The integral similarity function for moisture and tracers
570pure subroutine monin_obukhov_integral_tq_(stable_option, new_mo_option, rich_crit, zeta_trans, &
571 & n, psi_t, psi_q, zeta, zeta_t, zeta_q, &
572 & ln_z_zt, ln_z_zq, mask, ier)
573
574 integer, intent(in) :: stable_option
575 logical, intent(in) :: new_mo_option !miz
576 real(kind=fms_mo_kind_), intent(in) :: rich_crit, zeta_trans
577 integer, intent(in) :: n
578 real(kind=fms_mo_kind_), intent(inout), dimension(n) :: psi_t, psi_q
579 real(kind=fms_mo_kind_), intent(in) , dimension(n) :: zeta, zeta_t, zeta_q, ln_z_zt, ln_z_zq
580 logical, intent(in) , dimension(n) :: mask
581 integer, intent(out) :: ier
582
583 real(kind=fms_mo_kind_), dimension(n) :: x, x_t, x_q
584 logical, dimension(n) :: stable, unstable, &
585 weakly_stable, strongly_stable
586 real(kind=fms_mo_kind_) :: b_stab, lambda
587 real(kind=fms_mo_kind_) :: s3 !miz
588 integer, parameter :: lkind = fms_mo_kind_
589
590 ier = 0
591
592 b_stab = 1.0_lkind/rich_crit
593
594stable = mask .and. zeta >= 0.0_lkind
595unstable = mask .and. zeta < 0.0_lkind
596
597!miz: modified to include a new monin-obukhov option
598if (new_mo_option) then
599 s3 = sqrt(3.0_lkind)
600 where(unstable)
601 x = (1.0_lkind - 16.0_lkind*zeta)**(1.0_lkind/3.0_lkind)
602 x_t = (1.0_lkind - 16.0_lkind*zeta_t)**(1.0_lkind/3.0_lkind)
603 x_q = (1.0_lkind - 16.0_lkind*zeta_q)**(1.0_lkind/3.0_lkind)
604
605 psi_t = ln_z_zt - 1.5_lkind*log((x**2+x+1.0_lkind)/(x_t**2 + x_t + 1.0_lkind)) + s3 * &
606 (atan((2.0_lkind*x+1.0_lkind)/s3) - atan((2.0_lkind*x_t + 1.0_lkind)/s3))
607 psi_q = ln_z_zq - 1.5_lkind*log((x**2+x+1.0_lkind)/(x_q**2 + x_q + 1.0_lkind)) + s3 * &
608 (atan((2.0_lkind*x+1.0_lkind)/s3) - atan((2.0_lkind*x_q + 1.0_lkind)/s3))
609 end where
610else
611
612where(unstable)
613
614 x = sqrt(1.0_lkind - 16.0_lkind*zeta)
615 x_t = sqrt(1.0_lkind - 16.0_lkind*zeta_t)
616 x_q = sqrt(1.0_lkind - 16.0_lkind*zeta_q)
617
618 psi_t = ln_z_zt - 2.0_lkind*log( (1.0_lkind + x)/(1.0_lkind + x_t) )
619 psi_q = ln_z_zq - 2.0_lkind*log( (1.0_lkind + x)/(1.0_lkind + x_q) )
620
621end where
622end if
623!miz
624
625if( stable_option == 1) then
626
627 where (stable)
628
629 psi_t = ln_z_zt + (5.0_lkind - b_stab)*log((1.0_lkind + zeta)/(1.0_lkind + zeta_t)) &
630 + b_stab*(zeta - zeta_t)
631 psi_q = ln_z_zq + (5.0_lkind - b_stab)*log((1.0_lkind + zeta)/(1.0_lkind + zeta_q)) &
632 + b_stab*(zeta - zeta_q)
633
634 end where
635
636else if (stable_option == 2) then
637
638 lambda = 1.0_lkind + (5.0_lkind - b_stab)*zeta_trans
639
640 weakly_stable = stable .and. zeta <= zeta_trans
641 strongly_stable = stable .and. zeta > zeta_trans
642
643 where (weakly_stable)
644 psi_t = ln_z_zt + 5.0_lkind*(zeta - zeta_t)
645 psi_q = ln_z_zq + 5.0_lkind*(zeta - zeta_q)
646 end where
647
648 where(strongly_stable)
649 x = (lambda - 1.0_lkind)*log(zeta/zeta_trans) + b_stab*(zeta - zeta_trans)
650 endwhere
651
652 where (strongly_stable .and. zeta_t <= zeta_trans)
653 psi_t = ln_z_zt + x + 5.0_lkind*(zeta_trans - zeta_t)
654 end where
655 where (strongly_stable .and. zeta_t > zeta_trans)
656 psi_t = lambda*ln_z_zt + b_stab*(zeta - zeta_t)
657 endwhere
658
659 where (strongly_stable .and. zeta_q <= zeta_trans)
660 psi_q = ln_z_zq + x + 5.0_lkind*(zeta_trans - zeta_q)
661 end where
662 where (strongly_stable .and. zeta_q > zeta_trans)
663 psi_q = lambda*ln_z_zq + b_stab*(zeta - zeta_q)
664 endwhere
665
666end if
667
668end subroutine monin_obukhov_integral_tq_
669
670
671pure subroutine monin_obukhov_stable_mix_(stable_option, rich_crit, zeta_trans, &
672 & n, rich, mix, ier)
673
674 integer, intent(in) :: stable_option
675 real(kind=fms_mo_kind_), intent(in) :: rich_crit, zeta_trans
676 integer, intent(in) :: n
677 real(kind=fms_mo_kind_), intent(in), dimension(n) :: rich
678 real(kind=fms_mo_kind_), intent(out), dimension(n) :: mix
679 integer, intent(out) :: ier
680
681 real(kind=fms_mo_kind_) :: r, a, b, c, zeta, phi
682 real(kind=fms_mo_kind_) :: b_stab, rich_trans, lambda
683 integer :: i
684 integer, parameter :: lkind = fms_mo_kind_
685
686 ier = 0
687
688mix = 0.0_lkind
689b_stab = 1.0_lkind/rich_crit
690rich_trans = zeta_trans/(1.0_lkind + 5.0_lkind*zeta_trans)
691
692if(stable_option == 1) then
693
694 c = - 1.0_lkind
695 do i = 1, n
696 if(rich(i) > 0.0_lkind .and. rich(i) < rich_crit) then
697 r = 1.0_lkind/rich(i)
698 a = r - b_stab
699 b = r - (1.0_lkind + 5.0_lkind)
700 zeta = (-b + sqrt(b*b - 4.0_lkind*a*c))/(2.0_lkind*a)
701 phi = 1.0_lkind + b_stab*zeta + (5.0_lkind - b_stab)*zeta/(1.0_lkind + zeta)
702 mix(i) = 1.0_lkind/(phi*phi)
703 endif
704 end do
705
706else if(stable_option == 2) then
707
708 lambda = 1.0_lkind + (5.0_lkind - b_stab)*zeta_trans
709
710 where(rich > 0.0_lkind .and. rich <= rich_trans)
711 mix = (1.0_lkind - 5.0_lkind*rich)**2
712 end where
713 where(rich > rich_trans .and. rich < rich_crit)
714 mix = ((1.0_lkind - b_stab*rich)/lambda)**2
715 end where
716
717end if
718
719end subroutine monin_obukhov_stable_mix_
720
721!> @}
722! close documentation grouping
monin_obukhov_solve_zeta_
pure subroutine monin_obukhov_solve_zeta_(error, zeta_min, max_iter, small, stable_option, new_mo_option, rich_crit, zeta_trans, n, rich, z, z0, zt, zq, f_m, f_t, f_q, mask, ier)
Definition monin_obukhov_inter.inc:187
monin_obukhov_integral_m_
pure subroutine monin_obukhov_integral_m_(stable_option, rich_crit, zeta_trans, n, psi_m, zeta, zeta_0, ln_z_z0, mask, ier)
The integral similarity function for momentum.
Definition monin_obukhov_inter.inc:495
monin_obukhov_integral_tq_
pure subroutine monin_obukhov_integral_tq_(stable_option, new_mo_option, rich_crit, zeta_trans, n, psi_t, psi_q, zeta, zeta_t, zeta_q, ln_z_zt, ln_z_zq, mask, ier)
The integral similarity function for moisture and tracers.
Definition monin_obukhov_inter.inc:573
monin_obukhov_diff_
pure subroutine monin_obukhov_diff_(vonkarm, ustar_min, neutral, stable_option, new_mo_option, rich_crit, zeta_trans, ni, nj, nk, z, u_star, b_star, k_m, k_h, ier)
Definition monin_obukhov_inter.inc:28
monin_obukhov_drag_1d_
pure subroutine monin_obukhov_drag_1d_(grav, vonkarm, error, zeta_min, max_iter, small, neutral, stable_option, new_mo_option, rich_crit, zeta_trans, drag_min_heat, drag_min_moist, drag_min_mom, n, pt, pt0, z, z0, zt, zq, speed, drag_m, drag_t, drag_q, u_star, b_star, lavail, avail, ier)
Definition monin_obukhov_inter.inc:78
monin_obukhov_profile_1d_
pure subroutine monin_obukhov_profile_1d_(vonkarm, neutral, stable_option, new_mo_option, rich_crit, zeta_trans, n, zref, zref_t, z, z0, zt, zq, u_star, b_star, q_star, del_m, del_t, del_q, lavail, avail, ier)
Definition monin_obukhov_inter.inc:412
monin_obukhov_derivative_t_
pure subroutine monin_obukhov_derivative_t_(stable_option, new_mo_option, rich_crit, zeta_trans, n, phi_t, zeta, mask, ier)
The differential similarity function for buoyancy and tracers.
Definition monin_obukhov_inter.inc:304