FV3 Bundle
monin_obukhov_kernel.F90
Go to the documentation of this file.
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 ! -*-F90-*-
21 ! $Id$
22 
23 !==============================================================================
24 ! Kernel routine interface
25 !==============================================================================
26 
27 module monin_obukhov_inter
28 #include <fms_platform.h>
29 
30 ! explicit interface to all kernel routines
31 #include "monin_obukhov_interfaces.h"
32 
33 end module monin_obukhov_inter
34 
35 !==============================================================================
36 ! Kernel routines
37 !==============================================================================
38 
39 _pure subroutine monin_obukhov_diff(vonkarm, &
40  & ustar_min, &
41  & neutral, stable_option,new_mo_option,rich_crit, zeta_trans, &
42  & ni, nj, nk, z, u_star, b_star, k_m, k_h, ier)
43 
44  implicit none
45 
46  real , intent(in ) :: vonkarm
47  real , intent(in ) :: ustar_min ! = 1.e-10
48  logical, intent(in ) :: neutral
49  integer, intent(in ) :: stable_option
50  logical, intent(in ) :: new_mo_option !miz
51  real , intent(in ) :: rich_crit, zeta_trans
52  integer, intent(in ) :: ni, nj, nk
53  real , intent(in ), dimension(ni, nj, nk) :: z
54  real , intent(in ), dimension(ni, nj) :: u_star, b_star
55  real , intent( out), dimension(ni, nj, nk) :: k_m, k_h
56  integer, intent( out) :: ier
57 
58  real , dimension(ni, nj) :: phi_m, phi_h, zeta, uss
59  integer :: j, k
60 
61  logical, dimension(ni) :: mask
62 
63  interface
64  _pure subroutine monin_obukhov_derivative_t(stable_option,new_mo_option,rich_crit, zeta_trans, &
65  & n, phi_t, zeta, mask, ier)
66 
67  ! the differential similarity function for buoyancy and tracers
68  ! Note: seems to be the same as monin_obukhov_derivative_m?
69 
70  integer, intent(in ) :: stable_option
71  logical, intent(in ) :: new_mo_option !miz
72  real , intent(in ) :: rich_crit, zeta_trans
73  integer, intent(in ) :: n
74  real , intent( out), dimension(n) :: phi_t
75  real , intent(in ), dimension(n) :: zeta
76  logical, intent(in ), dimension(n) :: mask
77  integer, intent( out) :: ier
78  end subroutine monin_obukhov_derivative_t
79 
80  _pure subroutine monin_obukhov_derivative_m(stable_option, rich_crit, zeta_trans, &
81  & n, phi_m, zeta, mask, ier)
82 
83  ! the differential similarity function for momentum
84 
85  integer, intent(in ) :: stable_option
86  real , intent(in ) :: rich_crit, zeta_trans
87  integer, intent(in ) :: n
88  real , intent( out), dimension(n) :: phi_m
89  real , intent(in ), dimension(n) :: zeta
90  logical, intent(in ), dimension(n) :: mask
91  integer, intent(out ) :: ier
92  end subroutine monin_obukhov_derivative_m
93  end interface
94 
95  ier = 0
96 
97  mask = .true.
98  uss = max(u_star, ustar_min)
99 
100  if(neutral) then
101  do k = 1, size(z,3)
102  k_m(:,:,k) = vonkarm *uss*z(:,:,k)
103  k_h(:,:,k) = k_m(:,:,k)
104  end do
105  else
106  do k = 1, size(z,3)
107  zeta = - vonkarm * b_star*z(:,:,k)/(uss*uss)
108  do j = 1, size(z,2)
109  call monin_obukhov_derivative_m(stable_option, rich_crit, zeta_trans, &
110  & ni, phi_m(:,j), zeta(:,j), mask, ier)
111  call monin_obukhov_derivative_t(stable_option, new_mo_option,rich_crit, zeta_trans, &
112  & ni, phi_h(:,j), zeta(:,j), mask, ier)
113  enddo
114  k_m(:,:,k) = vonkarm * uss*z(:,:,k)/phi_m
115  k_h(:,:,k) = vonkarm * uss*z(:,:,k)/phi_h
116  end do
117  endif
118 
119 end subroutine monin_obukhov_diff
120 !==============================================================================
121 _pure subroutine monin_obukhov_drag_1d(grav, vonkarm, &
122  & error, zeta_min, max_iter, small, &
123  & neutral, stable_option, new_mo_option, rich_crit, zeta_trans,&
124  & drag_min_heat, drag_min_moist, drag_min_mom, &
125  & n, pt, pt0, z, z0, zt, zq, speed, drag_m, drag_t, &
126  & drag_q, u_star, b_star, lavail, avail, ier)
127 
128  implicit none
129 
130  real , intent(in ) :: grav
131  real , intent(in ) :: vonkarm
132  real , intent(in ) :: error ! = 1.e-04
133  real , intent(in ) :: zeta_min ! = 1.e-06
134  integer, intent(in ) :: max_iter ! = 20
135  real , intent(in ) :: small ! = 1.e-04
136  logical, intent(in ) :: neutral
137  integer, intent(in ) :: stable_option
138  logical, intent(in ) :: new_mo_option
139  real , intent(in ) :: rich_crit, zeta_trans
140  real , intent(in ) :: drag_min_heat, drag_min_moist, drag_min_mom
141  integer, intent(in ) :: n
142  real , intent(in ), dimension(n) :: pt, pt0, z, z0, zt, zq, speed
143  real , intent(inout), dimension(n) :: drag_m, drag_t, drag_q, u_star, b_star
144  logical, intent(in ) :: lavail ! whether to use provided mask or not
145  logical, intent(in ), dimension(n) :: avail ! provided mask
146  integer, intent(out ) :: ier
147 
148  real , dimension(n) :: rich, fm, ft, fq, zz
149  logical, dimension(n) :: mask, mask_1, mask_2
150  real , dimension(n) :: delta_b !!, us, bs, qs
151  real :: r_crit, sqrt_drag_min_heat
152  real :: sqrt_drag_min_moist, sqrt_drag_min_mom
153  real :: us, bs, qs
154  integer :: i
155 
156  interface
157  _pure subroutine monin_obukhov_solve_zeta(error, zeta_min, max_iter, small, &
158  & stable_option, new_mo_option, rich_crit, zeta_trans, &
159  & n, rich, z, z0, zt, zq, f_m, f_t, f_q, mask, ier)
160 
161  real , intent(in ) :: error ! = 1.e-04
162  real , intent(in ) :: zeta_min ! = 1.e-06
163  integer, intent(in ) :: max_iter ! = 20
164  real , intent(in ) :: small ! = 1.e-04
165  integer, intent(in ) :: stable_option
166  logical, intent(in ) :: new_mo_option
167  real , intent(in ) :: rich_crit, zeta_trans
168  integer, intent(in ) :: n
169  real , intent(in ), dimension(n) :: rich, z, z0, zt, zq
170  logical, intent(in ), dimension(n) :: mask
171  real , intent( out), dimension(n) :: f_m, f_t, f_q
172  integer, intent( out) :: ier
173  end subroutine monin_obukhov_solve_zeta
174  end interface
175 
176  ier = 0
177  r_crit = 0.95*rich_crit ! convergence can get slow if one is
178  ! close to rich_crit
179  sqrt_drag_min_heat = 0.0
180  if(drag_min_heat.ne.0.0) sqrt_drag_min_heat = sqrt(drag_min_heat)
181  sqrt_drag_min_moist = 0.0
182  if(drag_min_moist.ne.0.0) sqrt_drag_min_moist = sqrt(drag_min_moist)
183  sqrt_drag_min_mom = 0.0
184  if(drag_min_mom.ne.0.0) sqrt_drag_min_mom = sqrt(drag_min_mom)
185 
186  mask = .true.
187  if(lavail) mask = avail
188 
189  where(mask)
190  delta_b = grav*(pt0 - pt)/pt0
191  rich = - z*delta_b/(speed*speed + small)
192  zz = max(z,z0,zt,zq)
193  elsewhere
194  rich = 0.0
195  end where
196 
197  if(neutral) then
198 
199  do i = 1, n
200  if(mask(i)) then
201  fm(i) = log(zz(i)/z0(i))
202  ft(i) = log(zz(i)/zt(i))
203  fq(i) = log(zz(i)/zq(i))
204  us = vonkarm/fm(i)
205  bs = vonkarm/ft(i)
206  qs = vonkarm/fq(i)
207  drag_m(i) = us*us
208  drag_t(i) = us*bs
209  drag_q(i) = us*qs
210  u_star(i) = us*speed(i)
211  b_star(i) = bs*delta_b(i)
212  end if
213  enddo
214 
215  else
216 
217  mask_1 = mask .and. rich < r_crit
218  mask_2 = mask .and. rich >= r_crit
219 
220  do i = 1, n
221  if(mask_2(i)) then
222  drag_m(i) = drag_min_mom
223  drag_t(i) = drag_min_heat
224  drag_q(i) = drag_min_moist
225  us = sqrt_drag_min_mom
226  bs = sqrt_drag_min_heat
227  u_star(i) = us*speed(i)
228  b_star(i) = bs*delta_b(i)
229  end if
230  enddo
231 
232  call monin_obukhov_solve_zeta (error, zeta_min, max_iter, small, &
233  & stable_option, new_mo_option, rich_crit, zeta_trans, &
234  & n, rich, zz, z0, zt, zq, fm, ft, fq, mask_1, ier)
235 
236  do i = 1, n
237  if(mask_1(i)) then
238  us = max(vonkarm/fm(i), sqrt_drag_min_mom)
239  bs = max(vonkarm/ft(i), sqrt_drag_min_heat)
240  qs = max(vonkarm/fq(i), sqrt_drag_min_moist)
241  drag_m(i) = us*us
242  drag_t(i) = us*bs
243  drag_q(i) = us*qs
244  u_star(i) = us*speed(i)
245  b_star(i) = bs*delta_b(i)
246  endif
247  enddo
248 
249  end if
250 
251 end subroutine monin_obukhov_drag_1d
252 !==============================================================================
253 _pure subroutine monin_obukhov_solve_zeta(error, zeta_min, max_iter, small, &
254  & stable_option, new_mo_option, rich_crit, zeta_trans, & !miz
255  & n, rich, z, z0, zt, zq, f_m, f_t, f_q, mask, ier)
256 
257  implicit none
258 
259  real , intent(in ) :: error ! = 1.e-04
260  real , intent(in ) :: zeta_min ! = 1.e-06
261  integer, intent(in ) :: max_iter ! = 20
262  real , intent(in ) :: small ! = 1.e-04
263  integer, intent(in ) :: stable_option
264  logical, intent(in ) :: new_mo_option
265  real , intent(in ) :: rich_crit, zeta_trans
266  integer, intent(in ) :: n
267  real , intent(in ), dimension(n) :: rich, z, z0, zt, zq
268  logical, intent(in ), dimension(n) :: mask
269  real , intent( out), dimension(n) :: f_m, f_t, f_q
270  integer, intent( out) :: ier
271 
272 
273  real :: max_cor
274  integer :: iter
275 
276  real, dimension(n) :: &
277  d_rich, rich_1, correction, corr, z_z0, z_zt, z_zq, &
278  ln_z_z0, ln_z_zt, ln_z_zq, zeta, &
279  phi_m, phi_m_0, phi_t, phi_t_0, rzeta, &
280  zeta_0, zeta_t, zeta_q, df_m, df_t
281 
282  logical, dimension(n) :: mask_1
283 
284  interface
285  _pure subroutine monin_obukhov_derivative_t(stable_option,new_mo_option,rich_crit, zeta_trans, &
286  & n, phi_t, zeta, mask, ier)
287 
288  ! the differential similarity function for buoyancy and tracers
289  ! Note: seems to be the same as monin_obukhov_derivative_m?
290 
291  integer, intent(in ) :: stable_option
292  logical, intent(in ) :: new_mo_option !miz
293  real , intent(in ) :: rich_crit, zeta_trans
294  integer, intent(in ) :: n
295  real , intent( out), dimension(n) :: phi_t
296  real , intent(in ), dimension(n) :: zeta
297  logical, intent(in ), dimension(n) :: mask
298  integer, intent( out) :: ier
299  end subroutine monin_obukhov_derivative_t
300  _pure subroutine monin_obukhov_derivative_m(stable_option, rich_crit, zeta_trans, &
301  & n, phi_m, zeta, mask, ier)
302 
303  ! the differential similarity function for momentum
304 
305  integer, intent(in ) :: stable_option
306  real , intent(in ) :: rich_crit, zeta_trans
307  integer, intent(in ) :: n
308  real , intent( out), dimension(n) :: phi_m
309  real , intent(in ), dimension(n) :: zeta
310  logical, intent(in ), dimension(n) :: mask
311  integer, intent(out ) :: ier
312  end subroutine monin_obukhov_derivative_m
313  _pure subroutine monin_obukhov_integral_tq(stable_option, new_mo_option,rich_crit, zeta_trans, &
314  & n, psi_t, psi_q, zeta, zeta_t, zeta_q, &
315  & ln_z_zt, ln_z_zq, mask, ier)
316 
317  ! the integral similarity function for moisture and tracers
318 
319  integer, intent(in ) :: stable_option
320  logical, intent(in ) :: new_mo_option
321  real, intent(in ) :: rich_crit, zeta_trans
322  integer, intent(in ) :: n
323  real , intent(inout), dimension(n) :: psi_t, psi_q
324  real , intent(in) , dimension(n) :: zeta, zeta_t, zeta_q, ln_z_zt, ln_z_zq
325  logical, intent(in) , dimension(n) :: mask
326  integer, intent( out) :: ier
327  end subroutine monin_obukhov_integral_tq
328 
329  _pure subroutine monin_obukhov_integral_m(stable_option, rich_crit, zeta_trans, &
330  & n, psi_m, zeta, zeta_0, ln_z_z0, mask, ier)
331 
332  ! the integral similarity function for momentum
333 
334  integer, intent(in ) :: stable_option
335  real , intent(in ) :: rich_crit, zeta_trans
336  integer, intent(in ) :: n
337  real , intent(inout), dimension(n) :: psi_m
338  real , intent(in) , dimension(n) :: zeta, zeta_0, ln_z_z0
339  logical, intent(in) , dimension(n) :: mask
340  integer, intent(out) :: ier
341  end subroutine monin_obukhov_integral_m
342 
343  end interface
344 
345 
346  ier = 0
347 
348  z_z0 = z/z0
349  z_zt = z/zt
350  z_zq = z/zq
351  ln_z_z0 = log(z_z0)
352  ln_z_zt = log(z_zt)
353  ln_z_zq = log(z_zq)
354 
355  corr = 0.0
356  mask_1 = mask
357 
358  ! initial guess
359 
360  zeta = 0.0
361  where(mask_1)
362  zeta = rich*ln_z_z0*ln_z_z0/ln_z_zt
363  end where
364 
365  where (mask_1 .and. rich >= 0.0)
366  zeta = zeta/(1.0 - rich/rich_crit)
367  end where
368 
369  iter_loop: do iter = 1, max_iter
370 
371  where (mask_1 .and. abs(zeta).lt.zeta_min)
372  zeta = 0.0
373  f_m = ln_z_z0
374  f_t = ln_z_zt
375  f_q = ln_z_zq
376  mask_1 = .false. ! don't do any more calculations at these pts
377  end where
378 
379 
380  zeta_0 = 0.0
381  zeta_t = 0.0
382  zeta_q = 0.0
383  where (mask_1)
384  rzeta = 1.0/zeta
385  zeta_0 = zeta/z_z0
386  zeta_t = zeta/z_zt
387  zeta_q = zeta/z_zq
388  end where
389 
390  call monin_obukhov_derivative_m(stable_option, rich_crit, zeta_trans, &
391  & n, phi_m , zeta , mask_1, ier)
392  call monin_obukhov_derivative_m(stable_option, rich_crit, zeta_trans, &
393  & n, phi_m_0, zeta_0, mask_1, ier)
394  call monin_obukhov_derivative_t(stable_option, new_mo_option,rich_crit, zeta_trans, &
395  & n, phi_t , zeta , mask_1, ier)
396  call monin_obukhov_derivative_t(stable_option, new_mo_option,rich_crit, zeta_trans, &
397  & n, phi_t_0, zeta_t, mask_1, ier)
398 
399  call monin_obukhov_integral_m(stable_option, rich_crit, zeta_trans, &
400  & n, f_m, zeta, zeta_0, ln_z_z0, mask_1, ier)
401  call monin_obukhov_integral_tq(stable_option, new_mo_option, rich_crit, zeta_trans, &
402  & n, f_t, f_q, zeta, zeta_t, zeta_q, ln_z_zt, ln_z_zq, mask_1, ier)
403 
404  where (mask_1)
405  df_m = (phi_m - phi_m_0)*rzeta
406  df_t = (phi_t - phi_t_0)*rzeta
407  rich_1 = zeta*f_t/(f_m*f_m)
408  d_rich = rich_1*( rzeta + df_t/f_t - 2.0 *df_m/f_m)
409  correction = (rich - rich_1)/d_rich
410  corr = min(abs(correction),abs(correction/zeta))
411  ! the criterion corr < error seems to work ok, but is a bit arbitrary
412  ! when zeta is small the tolerance is reduced
413  end where
414 
415  max_cor= maxval(corr)
416 
417  if(max_cor > error) then
418  mask_1 = mask_1 .and. (corr > error)
419  ! change the mask so computation proceeds only on non-converged points
420  where(mask_1)
421  zeta = zeta + correction
422  end where
423  cycle iter_loop
424  else
425  return
426  end if
427 
428  end do iter_loop
429 
430  ier = 1 ! surface drag iteration did not converge
431 
432 end subroutine monin_obukhov_solve_zeta
433 !==============================================================================
434 _pure subroutine monin_obukhov_derivative_t(stable_option,new_mo_option,rich_crit, zeta_trans, &
435  & n, phi_t, zeta, mask, ier)
436 
437  ! the differential similarity function for buoyancy and tracers
438  ! Note: seems to be the same as monin_obukhov_derivative_m?
439 
440  implicit none
441 
442  integer, intent(in ) :: stable_option
443  logical, intent(in ) :: new_mo_option !miz
444  real , intent(in ) :: rich_crit, zeta_trans
445  integer, intent(in ) :: n
446  real , intent( out), dimension(n) :: phi_t
447  real , intent(in ), dimension(n) :: zeta
448  logical, intent(in ), dimension(n) :: mask
449  integer, intent( out) :: ier
450 
451  logical, dimension(n) :: stable, unstable
452  real :: b_stab, lambda
453 
454  ier = 0
455  b_stab = 1.0/rich_crit
456 
457  stable = mask .and. zeta >= 0.0
458  unstable = mask .and. zeta < 0.0
459 
460 !miz: modified to include new monin-obukhov option
461  if (new_mo_option) then
462  where (unstable)
463  phi_t = (1 - 16.0*zeta)**(-1./3.)
464  end where
465  else
466  where (unstable)
467  phi_t = (1 - 16.0*zeta)**(-0.5)
468  end where
469  end if
470 !miz
471 
472  if(stable_option == 1) then
473 
474  where (stable)
475  phi_t = 1.0 + zeta*(5.0 + b_stab*zeta)/(1.0 + zeta)
476  end where
477 
478  else if(stable_option == 2) then
479 
480  lambda = 1.0 + (5.0 - b_stab)*zeta_trans
481 
482  where (stable .and. zeta < zeta_trans)
483  phi_t = 1 + 5.0*zeta
484  end where
485  where (stable .and. zeta >= zeta_trans)
486  phi_t = lambda + b_stab*zeta
487  end where
488 
489  endif
490 
491 end subroutine monin_obukhov_derivative_t
492 !==============================================================================
493 _pure subroutine monin_obukhov_derivative_m(stable_option, rich_crit, zeta_trans, &
494  & n, phi_m, zeta, mask, ier)
495 
496  ! the differential similarity function for momentum
497 
498  implicit none
499 
500  integer, intent(in ) :: stable_option
501  real , intent(in ) :: rich_crit, zeta_trans
502  integer, intent(in ) :: n
503  real , intent( out), dimension(n) :: phi_m
504  real , intent(in ), dimension(n) :: zeta
505  logical, intent(in ), dimension(n) :: mask
506  integer, intent(out ) :: ier
507 
508  logical, dimension(n) :: stable, unstable
509  real , dimension(n) :: x
510  real :: b_stab, lambda
511 
512 
513  ier = 0
514  b_stab = 1.0/rich_crit
515 
516  stable = mask .and. zeta >= 0.0
517  unstable = mask .and. zeta < 0.0
518 
519  where (unstable)
520  x = (1 - 16.0*zeta )**(-0.5)
521  phi_m = sqrt(x) ! phi_m = (1 - 16.0*zeta)**(-0.25)
522  end where
523 
524  if(stable_option == 1) then
525 
526  where (stable)
527  phi_m = 1.0 + zeta *(5.0 + b_stab*zeta)/(1.0 + zeta)
528  end where
529 
530  else if(stable_option == 2) then
531 
532  lambda = 1.0 + (5.0 - b_stab)*zeta_trans
533 
534  where (stable .and. zeta < zeta_trans)
535  phi_m = 1 + 5.0*zeta
536  end where
537  where (stable .and. zeta >= zeta_trans)
538  phi_m = lambda + b_stab*zeta
539  end where
540 
541  endif
542 
543 end subroutine monin_obukhov_derivative_m
544 !==============================================================================
545 _pure subroutine monin_obukhov_profile_1d(vonkarm, &
546  & neutral, stable_option, new_mo_option, rich_crit, zeta_trans, &
547  & n, zref, zref_t, z, z0, zt, zq, u_star, b_star, q_star, &
548  & del_m, del_t, del_q, lavail, avail, ier)
549 
550  implicit none
551 
552  real , intent(in ) :: vonkarm
553  logical, intent(in ) :: neutral
554  integer, intent(in ) :: stable_option
555  logical, intent(in ) :: new_mo_option
556  real , intent(in ) :: rich_crit, zeta_trans
557  integer, intent(in ) :: n
558  real, intent(in ) :: zref, zref_t
559  real, intent(in ), dimension(n) :: z, z0, zt, zq, u_star, b_star, q_star
560  real, intent( out), dimension(n) :: del_m, del_t, del_q
561  logical, intent(in ) :: lavail ! whether to use provided mask or not
562  logical, intent(in ), dimension(n) :: avail ! provided mask
563  integer, intent(out ) :: ier
564 
565  real, dimension(n) :: zeta, zeta_0, zeta_t, zeta_q, zeta_ref, zeta_ref_t, &
566  ln_z_z0, ln_z_zt, ln_z_zq, ln_z_zref, ln_z_zref_t, &
567  f_m_ref, f_m, f_t_ref, f_t, f_q_ref, f_q, &
568  mo_length_inv
569 
570  logical, dimension(n) :: mask
571 
572  interface
573  _pure subroutine monin_obukhov_integral_tq(stable_option, new_mo_option,rich_crit, zeta_trans, &
574  & n, psi_t, psi_q, zeta, zeta_t, zeta_q, &
575  & ln_z_zt, ln_z_zq, mask, ier)
576 
577  ! the integral similarity function for moisture and tracers
578 
579  integer, intent(in ) :: stable_option
580  logical, intent(in ) :: new_mo_option
581  real, intent(in ) :: rich_crit, zeta_trans
582  integer, intent(in ) :: n
583  real , intent(inout), dimension(n) :: psi_t, psi_q
584  real , intent(in) , dimension(n) :: zeta, zeta_t, zeta_q, ln_z_zt, ln_z_zq
585  logical, intent(in) , dimension(n) :: mask
586  integer, intent( out) :: ier
587  end subroutine monin_obukhov_integral_tq
588 
589  _pure subroutine monin_obukhov_integral_m(stable_option, rich_crit, zeta_trans, &
590  & n, psi_m, zeta, zeta_0, ln_z_z0, mask, ier)
591 
592  ! the integral similarity function for momentum
593 
594  integer, intent(in ) :: stable_option
595  real , intent(in ) :: rich_crit, zeta_trans
596  integer, intent(in ) :: n
597  real , intent(inout), dimension(n) :: psi_m
598  real , intent(in) , dimension(n) :: zeta, zeta_0, ln_z_z0
599  logical, intent(in) , dimension(n) :: mask
600  integer, intent(out) :: ier
601  end subroutine monin_obukhov_integral_m
602  end interface
603 
604  ier = 0
605 
606  mask = .true.
607  if(lavail) mask = avail
608 
609  del_m = 0.0 ! zero output arrays
610  del_t = 0.0
611  del_q = 0.0
612 
613  where(mask)
614  ln_z_z0 = log(z/z0)
615  ln_z_zt = log(z/zt)
616  ln_z_zq = log(z/zq)
617  ln_z_zref = log(z/zref)
618  ln_z_zref_t = log(z/zref_t)
619  endwhere
620 
621  if(neutral) then
622 
623  where(mask)
624  del_m = 1.0 - ln_z_zref /ln_z_z0
625  del_t = 1.0 - ln_z_zref_t/ln_z_zt
626  del_q = 1.0 - ln_z_zref_t/ln_z_zq
627  endwhere
628 
629  else
630 
631  where(mask .and. u_star > 0.0)
632  mo_length_inv = - vonkarm * b_star/(u_star*u_star)
633  zeta = z *mo_length_inv
634  zeta_0 = z0 *mo_length_inv
635  zeta_t = zt *mo_length_inv
636  zeta_q = zq *mo_length_inv
637  zeta_ref = zref *mo_length_inv
638  zeta_ref_t = zref_t*mo_length_inv
639  endwhere
640 
641  call monin_obukhov_integral_m(stable_option, rich_crit, zeta_trans, &
642  & n, f_m, zeta, zeta_0, ln_z_z0, mask, ier)
643  call monin_obukhov_integral_m(stable_option, rich_crit, zeta_trans, &
644  & n, f_m_ref, zeta, zeta_ref, ln_z_zref, mask, ier)
645 
646  call monin_obukhov_integral_tq(stable_option, new_mo_option, rich_crit, zeta_trans, &
647  & n, f_t, f_q, zeta, zeta_t, zeta_q, ln_z_zt, ln_z_zq, mask, ier)
648  call monin_obukhov_integral_tq(stable_option, new_mo_option, rich_crit, zeta_trans, &
649  & n, f_t_ref, f_q_ref, zeta, zeta_ref_t, zeta_ref_t, ln_z_zref_t, ln_z_zref_t, mask, ier)
650 
651  where(mask)
652  del_m = 1.0 - f_m_ref/f_m
653  del_t = 1.0 - f_t_ref/f_t
654  del_q = 1.0 - f_q_ref/f_q
655  endwhere
656 
657  end if
658 
659 
660 end subroutine monin_obukhov_profile_1d
661 !==============================================================================
662 _pure subroutine monin_obukhov_integral_m(stable_option, rich_crit, zeta_trans, &
663  & n, psi_m, zeta, zeta_0, ln_z_z0, mask, ier)
664 
665  ! the integral similarity function for momentum
666 
667  implicit none
668 
669  integer, intent(in ) :: stable_option
670  real , intent(in ) :: rich_crit, zeta_trans
671  integer, intent(in ) :: n
672  real , intent(inout), dimension(n) :: psi_m
673  real , intent(in) , dimension(n) :: zeta, zeta_0, ln_z_z0
674  logical, intent(in) , dimension(n) :: mask
675  integer, intent(out) :: ier
676 
677  real :: b_stab, lambda
678 
679  real, dimension(n) :: x, x_0, x1, x1_0, num, denom, y
680  logical, dimension(n) :: stable, unstable, &
681  weakly_stable, strongly_stable
682 
683  ier = 0
684 
685  b_stab = 1.0/rich_crit
686 
687  stable = mask .and. zeta >= 0.0
688  unstable = mask .and. zeta < 0.0
689 
690  where(unstable)
691 
692  x = sqrt(1 - 16.0*zeta)
693  x_0 = sqrt(1 - 16.0*zeta_0)
694 
695  x = sqrt(x)
696  x_0 = sqrt(x_0)
697 
698  x1 = 1.0 + x
699  x1_0 = 1.0 + x_0
700 
701  num = x1*x1*(1.0 + x*x)
702  denom = x1_0*x1_0*(1.0 + x_0*x_0)
703  y = atan(x) - atan(x_0)
704  psi_m = ln_z_z0 - log(num/denom) + 2*y
705 
706  end where
707 
708  if( stable_option == 1) then
709 
710  where (stable)
711  psi_m = ln_z_z0 + (5.0 - b_stab)*log((1.0 + zeta)/(1.0 + zeta_0)) &
712  + b_stab*(zeta - zeta_0)
713  end where
714 
715  else if (stable_option == 2) then
716 
717  lambda = 1.0 + (5.0 - b_stab)*zeta_trans
718 
719  weakly_stable = stable .and. zeta <= zeta_trans
720  strongly_stable = stable .and. zeta > zeta_trans
721 
722  where (weakly_stable)
723  psi_m = ln_z_z0 + 5.0*(zeta - zeta_0)
724  end where
725 
726  where(strongly_stable)
727  x = (lambda - 1.0)*log(zeta/zeta_trans) + b_stab*(zeta - zeta_trans)
728  endwhere
729 
730  where (strongly_stable .and. zeta_0 <= zeta_trans)
731  psi_m = ln_z_z0 + x + 5.0*(zeta_trans - zeta_0)
732  end where
733  where (strongly_stable .and. zeta_0 > zeta_trans)
734  psi_m = lambda*ln_z_z0 + b_stab*(zeta - zeta_0)
735  endwhere
736 
737  end if
738 
739 end subroutine monin_obukhov_integral_m
740 !==============================================================================
741 _pure subroutine monin_obukhov_integral_tq(stable_option, new_mo_option, rich_crit, zeta_trans, &
742  & n, psi_t, psi_q, zeta, zeta_t, zeta_q, &
743  & ln_z_zt, ln_z_zq, mask, ier)
744 
745  ! the integral similarity function for moisture and tracers
746 
747  implicit none
748 
749  integer, intent(in ) :: stable_option
750  logical, intent(in ) :: new_mo_option !miz
751  real, intent(in ) :: rich_crit, zeta_trans
752  integer, intent(in ) :: n
753  real , intent(inout), dimension(n) :: psi_t, psi_q
754  real , intent(in) , dimension(n) :: zeta, zeta_t, zeta_q, ln_z_zt, ln_z_zq
755  logical, intent(in) , dimension(n) :: mask
756  integer, intent( out) :: ier
757 
758  real, dimension(n) :: x, x_t, x_q
759  logical, dimension(n) :: stable, unstable, &
760  weakly_stable, strongly_stable
761  real :: b_stab, lambda
762  real :: s3 !miz
763 
764  ier = 0
765 
766  b_stab = 1.0/rich_crit
767 
768 stable = mask .and. zeta >= 0.0
769 unstable = mask .and. zeta < 0.0
770 
771 !miz: modified to include a new monin-obukhov option
772 if (new_mo_option) then
773  s3 = sqrt(3.0)
774  where(unstable)
775  x = (1 - 16.0*zeta)**(1./3.)
776  x_t = (1 - 16.0*zeta_t)**(1./3.)
777  x_q = (1 - 16.0*zeta_q)**(1./3.)
778 
779  psi_t = ln_z_zt - 1.5*log((x**2+x+1)/(x_t**2 + x_t + 1)) + s3*(atan((2*x+1)/s3) - atan((2*x_t + 1)/s3))
780  psi_q = ln_z_zq - 1.5*log((x**2+x+1)/(x_q**2 + x_q + 1)) + s3*(atan((2*x+1)/s3) - atan((2*x_q + 1)/s3))
781  end where
782 else
783 
784 where(unstable)
785 
786  x = sqrt(1 - 16.0*zeta)
787  x_t = sqrt(1 - 16.0*zeta_t)
788  x_q = sqrt(1 - 16.0*zeta_q)
789 
790  psi_t = ln_z_zt - 2.0*log( (1.0 + x)/(1.0 + x_t) )
791  psi_q = ln_z_zq - 2.0*log( (1.0 + x)/(1.0 + x_q) )
792 
793 end where
794 end if
795 !miz
796 
797 if( stable_option == 1) then
798 
799  where (stable)
800 
801  psi_t = ln_z_zt + (5.0 - b_stab)*log((1.0 + zeta)/(1.0 + zeta_t)) &
802  + b_stab*(zeta - zeta_t)
803  psi_q = ln_z_zq + (5.0 - b_stab)*log((1.0 + zeta)/(1.0 + zeta_q)) &
804  + b_stab*(zeta - zeta_q)
805 
806  end where
807 
808 else if (stable_option == 2) then
809 
810  lambda = 1.0 + (5.0 - b_stab)*zeta_trans
811 
812  weakly_stable = stable .and. zeta <= zeta_trans
813  strongly_stable = stable .and. zeta > zeta_trans
814 
815  where (weakly_stable)
816  psi_t = ln_z_zt + 5.0*(zeta - zeta_t)
817  psi_q = ln_z_zq + 5.0*(zeta - zeta_q)
818  end where
819 
820  where(strongly_stable)
821  x = (lambda - 1.0)*log(zeta/zeta_trans) + b_stab*(zeta - zeta_trans)
822  endwhere
823 
824  where (strongly_stable .and. zeta_t <= zeta_trans)
825  psi_t = ln_z_zt + x + 5.0*(zeta_trans - zeta_t)
826  end where
827  where (strongly_stable .and. zeta_t > zeta_trans)
828  psi_t = lambda*ln_z_zt + b_stab*(zeta - zeta_t)
829  endwhere
830 
831  where (strongly_stable .and. zeta_q <= zeta_trans)
832  psi_q = ln_z_zq + x + 5.0*(zeta_trans - zeta_q)
833  end where
834  where (strongly_stable .and. zeta_q > zeta_trans)
835  psi_q = lambda*ln_z_zq + b_stab*(zeta - zeta_q)
836  endwhere
837 
838 end if
839 
840 end subroutine monin_obukhov_integral_tq
841 !==============================================================================
842 _pure subroutine monin_obukhov_stable_mix(stable_option, rich_crit, zeta_trans, &
843  & n, rich, mix, ier)
844 
845  implicit none
846 
847  integer, intent(in ) :: stable_option
848  real , intent(in ) :: rich_crit, zeta_trans
849  integer, intent(in ) :: n
850  real , intent(in ), dimension(n) :: rich
851  real , intent( out), dimension(n) :: mix
852  integer, intent( out) :: ier
853 
854  real :: r, a, b, c, zeta, phi
855  real :: b_stab, rich_trans, lambda
856  integer :: i
857 
858  ier = 0
859 
860 mix = 0.0
861 b_stab = 1.0/rich_crit
862 rich_trans = zeta_trans/(1.0 + 5.0*zeta_trans)
863 
864 if(stable_option == 1) then
865 
866  c = - 1.0
867  do i = 1, n
868  if(rich(i) > 0.0 .and. rich(i) < rich_crit) then
869  r = 1.0/rich(i)
870  a = r - b_stab
871  b = r - (1.0 + 5.0)
872  zeta = (-b + sqrt(b*b - 4.0*a*c))/(2.0*a)
873  phi = 1.0 + b_stab*zeta + (5.0 - b_stab)*zeta/(1.0 + zeta)
874  mix(i) = 1./(phi*phi)
875  endif
876  end do
877 
878 else if(stable_option == 2) then
879 
880  lambda = 1.0 + (5.0 - b_stab)*zeta_trans
881 
882  where(rich > 0.0 .and. rich <= rich_trans)
883  mix = (1.0 - 5.0*rich)**2
884  end where
885  where(rich > rich_trans .and. rich < rich_crit)
886  mix = ((1.0 - b_stab*rich)/lambda)**2
887  end where
888 
889 end if
890 
891 end subroutine monin_obukhov_stable_mix
892 
893 
894 #ifdef _TEST_MONIN_OBUKHOV
895 !==============================================================================
896 ! Unit test
897 !==============================================================================
898 
899 program test
900 
901  use monin_obukhov_inter
902 
903  implicit none
904  integer, parameter :: i8 = selected_int_kind(18)
905  integer(i8) :: ier_tot, ier
906 
907  real :: grav, vonkarm, error, zeta_min, small, ustar_min
908  real :: zref, zref_t
909  integer :: max_iter
910 
911  real :: rich_crit, zeta_trans
912  real :: drag_min_heat, drag_min_moist, drag_min_mom
913  logical :: neutral
914  integer :: stable_option
915  logical :: new_mo_option
916 
917  grav = 9.80
918  vonkarm = 0.4
919  error = 1.0e-4
920  zeta_min = 1.0e-6
921  max_iter = 20
922  small = 1.0e-4
923  neutral = .false.
924  stable_option = 1
925  new_mo_option = .false.
926  rich_crit =10.0
927  zeta_trans = 0.5
928  drag_min_heat = 1.0e-5
929  drag_min_moist= 1.0e-5
930  drag_min_mom = 1.0e-5
931  ustar_min = 1.e-10
932 
933 
934  zref = 10.
935  zref_t = 2.
936 
937 
938  ier_tot = 0
939  call test_drag
940  print *,'test_drag ier = ', ier
941  ier_tot = ier_tot + ier
942 
943  call test_stable_mix
944  print *,'test_stable_mix ier = ', ier
945  ier_tot = ier_tot + ier
946 
947  call test_diff
948  print *,'test_diff ier = ', ier
949  ier_tot = ier_tot + ier
950 
951  call test_profile
952  print *,'test_profile ier = ', ier
953  ier_tot = ier_tot + ier
954 
955  if(ier_tot/=0) then
956  print *, ier_tot, '***ERRORS detected***'
957  else
958  print *,'No error detected.'
959  endif
960 
961  CONTAINS
962 
963  subroutine test_drag
964 
965  integer(i8) :: w
966 
967  integer :: i, ier_l
968  integer, parameter :: n = 5
969  logical :: avail(n), lavail
970 
971  real, dimension(n) :: pt, pt0, z, z0, zt, zq, speed, drag_m, drag_t, drag_q, u_star, b_star
972 
973  ! potential temperature
974  pt = (/ 268.559120403867, 269.799228886728, 277.443023238556, 295.79192777341, 293.268717243262 /)
975  pt0 = (/ 273.42369841804 , 272.551410044203, 278.638168565727, 298.133068766049, 292.898163706587/)
976  z = (/ 29.432779269303, 30.0497139076724, 31.6880000418153, 34.1873479240475, 33.2184943356517/)
977  z0 = (/ 5.86144925739178e-05, 0.0001, 0.000641655193293549, 3.23383768877187e-05, 0.07/)
978  zt = (/ 3.69403636275411e-05, 0.0001, 1.01735489109205e-05, 7.63933834969505e-05, 0.00947346982656289/)
979  zq = (/ 5.72575636226887e-05, 0.0001, 5.72575636226887e-05, 5.72575636226887e-05, 5.72575636226887e-05/)
980  speed = (/ 2.9693638452068, 2.43308757772094, 5.69418282305367, 9.5608693754561, 4.35302260074334/)
981  lavail = .true.
982  avail = (/.true., .true., .true., .true., .true. /)
983 
984  drag_m = 0
985  drag_t = 0
986  drag_q = 0
987  u_star = 0
988  b_star = 0
989 
990  call monin_obukhov_drag_1d(grav, vonkarm, &
991  & error, zeta_min, max_iter, small, &
992  & neutral, stable_option, new_mo_option, rich_crit, zeta_trans,&
993  & drag_min_heat, drag_min_moist, drag_min_mom, &
994  & n, pt, pt0, z, z0, zt, zq, speed, drag_m, drag_t, &
995  & drag_q, u_star, b_star, lavail, avail, ier_l)
996 
997  ! check sum results
998  w = 0
999  w = w + transfer(sum(drag_m), w)
1000  w = w + transfer(sum(drag_t), w)
1001  w = w + transfer(sum(drag_q), w)
1002  w = w + transfer(sum(u_star), w)
1003  w = w + transfer(sum(b_star), w)
1004 
1005  ! plug in check sum here>>>
1006 #if defined(__INTEL_COMPILER) || defined(_LF95)
1007 #define CHKSUM_DRAG 4466746452959549648
1008 #endif
1009 #if defined(_PGF95)
1010 #define CHKSUM_DRAG 4466746452959549650
1011 #endif
1012 
1013 
1014  print *,'chksum test_drag : ', w, ' ref ', chksum_drag
1015  ier = chksum_drag - w
1016 
1017  end subroutine test_drag
1018 
1019  subroutine test_stable_mix
1020 
1021  integer(i8) :: w
1022 
1023  integer, parameter :: n = 5
1024  real , dimension(n) :: rich
1025  real , dimension(n) :: mix
1026  integer :: ier_l
1027 
1028 
1029  stable_option = 1
1030  rich_crit = 10.0
1031  zeta_trans = 0.5
1032 
1033  rich = (/1650.92431853365, 1650.9256285137, 77.7636819036559, 1.92806556391324, 0.414767442012442/)
1034 
1035 
1036  call monin_obukhov_stable_mix(stable_option, rich_crit, zeta_trans, &
1037  & n, rich, mix, ier_l)
1038 
1039  w = transfer( sum(mix) , w)
1040 
1041  ! plug in check sum here>>>
1042 #if defined(__INTEL_COMPILER) || defined(_LF95)
1043 #define CHKSUM_STABLE_MIX 4590035772608644256
1044 #endif
1045 #if defined(_PGF95)
1046 #define CHKSUM_STABLE_MIX 4590035772608644258
1047 #endif
1048 
1049  print *,'chksum test_stable_mix: ', w, ' ref ', chksum_stable_mix
1050  ier = chksum_stable_mix - w
1051 
1052  end subroutine test_stable_mix
1053 
1054  !========================================================================
1055 
1056  subroutine test_diff
1057 
1058  integer(i8) :: w
1059 
1060  integer, parameter :: ni=1, nj=1, nk=1
1061  real , dimension(ni, nj, nk) :: z
1062  real , dimension(ni, nj) :: u_star, b_star
1063  real , dimension(ni, nj, nk) :: k_m, k_h
1064  integer :: ier_l
1065 
1066  z = 19.9982554527751
1067  u_star = 0.129638955971075
1068  b_star = 0.000991799765557209
1069 
1070  call monin_obukhov_diff(vonkarm, &
1071  & ustar_min, &
1072  & neutral, stable_option, new_mo_option, rich_crit, zeta_trans, &!miz
1073  & ni, nj, nk, z, u_star, b_star, k_m, k_h, ier_l)
1074 
1075  w = 0
1076  w = w + transfer( sum(k_m) , w)
1077  w = w + transfer( sum(k_h) , w)
1078 
1079  ! plug check sum in here>>>
1080 #if defined(__INTEL_COMPILER) || defined(_LF95) || defined(_PGF95)
1081 #define CHKSUM_DIFF -9222066590093362639
1082 #endif
1083 
1084  print *,'chksum test_diff : ', w, ' ref ', chksum_diff
1085 
1086  ier = chksum_diff - w
1087 
1088  end subroutine test_diff
1089 
1090  !========================================================================
1091 
1092  subroutine test_profile
1093 
1094  integer(i8) :: w
1095 
1096  integer, parameter :: n = 5
1097  integer :: ier_l
1098 
1099  logical :: avail(n)
1100 
1101  real, dimension(n) :: z, z0, zt, zq, u_star, b_star, q_star
1102  real, dimension(n) :: del_m, del_t, del_q
1103 
1104  z = (/ 29.432779269303, 30.0497139076724, 31.6880000418153, 34.1873479240475, 33.2184943356517 /)
1105  z0 = (/ 5.86144925739178e-05, 0.0001, 0.000641655193293549, 3.23383768877187e-05, 0.07/)
1106  zt = (/ 3.69403636275411e-05, 0.0001, 1.01735489109205e-05, 7.63933834969505e-05, 0.00947346982656289/)
1107  zq = (/ 5.72575636226887e-05, 0.0001, 5.72575636226887e-05, 5.72575636226887e-05, 5.72575636226887e-05/)
1108  u_star = (/ 0.109462510724615, 0.0932942802513508, 0.223232887323184, 0.290918439028557, 0.260087579361467/)
1109  b_star = (/ 0.00690834676781433, 0.00428178089592372, 0.00121229800895103, 0.00262353784027441, -0.000570314880866852/)
1110  q_star = (/ 0.000110861442197537, 9.44983279664197e-05, 4.17643828631936e-05, 0.000133135421415819, 9.36317815993945e-06/)
1111 
1112  avail = (/ .true., .true.,.true.,.true.,.true. /)
1113 
1114  call monin_obukhov_profile_1d(vonkarm, &
1115  & neutral, stable_option, new_mo_option, rich_crit, zeta_trans, &
1116  & n, zref, zref_t, z, z0, zt, zq, u_star, b_star, q_star, &
1117  & del_m, del_t, del_q, .true., avail, ier_l)
1118 
1119  ! check sum results
1120  w = 0
1121  w = w + transfer(sum(del_m), w)
1122  w = w + transfer(sum(del_t), w)
1123  w = w + transfer(sum(del_q), w)
1124 
1125  ! plug check sum in here>>>
1126 #if defined(__INTEL_COMPILER) || defined(_LF95)
1127 #define CHKSUM_PROFILE -4596910845317820786
1128 #endif
1129 #if defined(_PGF95)
1130 #define CHKSUM_PROFILE -4596910845317820785
1131 #endif
1132 
1133  print *,'chksum test_profile : ', w, ' ref ', chksum_profile
1134 
1135  ier = chksum_profile - w
1136 
1137  end subroutine test_profile
1138 
1139 
1140 end program test
1141 
1142 !==============================================================================
1143 
1144 #endif
1145 ! _TEST_MONIN_OBUKHOV
monin_obukhov_inter
Definition: monin_obukhov_kernel.F90:27
test_profile
subroutine test_profile
Definition: monin_obukhov_kernel.F90:1093
test_drag
subroutine test_drag
Definition: monin_obukhov_kernel.F90:964
monin_obukhov_diff
*f90 *! $Id interface _PURE subroutine monin_obukhov_diff(vonkarm, &&ustar_min, &&neutral, stable_option, new_mo_option, rich_crit, zeta_trans, &!miz &ni, nj, nk, z, u_star, b_star, k_m, k_h, ier) real
test_stable_mix
subroutine test_stable_mix
Definition: monin_obukhov_kernel.F90:1020
max
#define max(a, b)
Definition: mosaic_util.h:33
min
#define min(a, b)
Definition: mosaic_util.h:32
test_diff
subroutine test_diff
Definition: monin_obukhov_kernel.F90:1057