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