FMS  2024.03
Flexible Modeling System
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 
24 pure 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 
69 end subroutine monin_obukhov_diff_
70 
71 
72 pure 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 
181 end subroutine monin_obukhov_drag_1d_
182 
183 
184 pure 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 
297 end 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?
302 pure 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 
355 end subroutine monin_obukhov_derivative_t_
356 
357 
358 ! the differential similarity function for momentum
359 pure 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 
405 end subroutine monin_obukhov_derivative_m_
406 
407 
408 pure 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 
489 end subroutine monin_obukhov_profile_1d_
490 
491 
492 !> The integral similarity function for momentum
493 pure 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 
566 end subroutine monin_obukhov_integral_m_
567 
568 
569 !> The integral similarity function for moisture and tracers
570 pure 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 
594 stable = mask .and. zeta >= 0.0_lkind
595 unstable = mask .and. zeta < 0.0_lkind
596 
597 !miz: modified to include a new monin-obukhov option
598 if (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
610 else
611 
612 where(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 
621 end where
622 end if
623 !miz
624 
625 if( 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 
636 else 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 
666 end if
667 
668 end subroutine monin_obukhov_integral_tq_
669 
670 
671 pure 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 
688 mix = 0.0_lkind
689 b_stab = 1.0_lkind/rich_crit
690 rich_trans = zeta_trans/(1.0_lkind + 5.0_lkind*zeta_trans)
691 
692 if(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 
706 else 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 
717 end if
718 
719 end 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