SCALE-RM
scale_bulkflux.F90
Go to the documentation of this file.
1 !-------------------------------------------------------------------------------
10 !-------------------------------------------------------------------------------
11 #include "scalelib.h"
12 module scale_bulkflux
13  !-----------------------------------------------------------------------------
14  !
15  !++ used modules
16  !
17  use scale_precision
18  use scale_io
19  use scale_prof
20  !-----------------------------------------------------------------------------
21  implicit none
22  private
23  !-----------------------------------------------------------------------------
24  !
25  !++ Public procedure
26  !
27  public :: bulkflux_setup
28  public :: bulkflux_diagnose_scales
29  public :: bulkflux_diagnose_surface
30 
31  interface bulkflux_diagnose_scales
32  module procedure bulkflux_diagnose_scales_0d
33  module procedure bulkflux_diagnose_scales_2d
34  end interface bulkflux_diagnose_scales
35 
36  interface bulkflux_diagnose_surface
37  module procedure bulkflux_diagnose_surface_0d
38  module procedure bulkflux_diagnose_surface_2d
39  end interface bulkflux_diagnose_surface
40 
41  abstract interface
42  subroutine bc( &
43  T1, T0, &
44  P1, P0, &
45  Q1, Q0, &
46  Uabs, Z1, PBL, &
47  Z0M, Z0H, Z0E, &
48  Ustar, Tstar, Qstar, &
49  Wstar, RLmo, Ra, &
50  FracU10, FracT2, FracQ2 )
51  use scale_precision
52  implicit none
53  real(RP), intent(in) :: T1 ! tempearature at the lowest atmospheric layer [K]
54  real(RP), intent(in) :: T0 ! skin temperature [K]
55  real(RP), intent(in) :: P1 ! pressure at the lowest atmospheric layer [Pa]
56  real(RP), intent(in) :: P0 ! surface pressure [Pa]
57  real(RP), intent(in) :: Q1 ! mixing ratio at the lowest atmospheric layer [kg/kg]
58  real(RP), intent(in) :: Q0 ! surface mixing ratio [kg/kg]
59  real(RP), intent(in) :: Uabs! absolute velocity at the lowest atmospheric layer [m/s]
60  real(RP), intent(in) :: Z1 ! height at the lowest atmospheric layer [m]
61  real(RP), intent(in) :: PBL ! the top of atmospheric mixing layer [m]
62  real(RP), intent(in) :: Z0M ! roughness length of momentum [m]
63  real(RP), intent(in) :: Z0H ! roughness length of heat [m]
64  real(RP), intent(in) :: Z0E ! roughness length of moisture [m]
65 
66  real(RP), intent(out) :: Ustar ! friction velocity [m/s]
67  real(RP), intent(out) :: Tstar ! friction temperature [K]
68  real(RP), intent(out) :: Qstar ! friction mixing rate [kg/kg]
69  real(RP), intent(out) :: Wstar ! free convection velocity scale [m/s]
70  real(RP), intent(out) :: RLmo ! inversed Obukhov length [1/m]
71  real(RP), intent(out) :: Ra ! Aerodynamic resistance (=1/Ce)
72  real(RP), intent(out) :: FracU10 ! calculation parameter for U10 [-]
73  real(RP), intent(out) :: FracT2 ! calculation parameter for T2 [-]
74  real(RP), intent(out) :: FracQ2 ! calculation parameter for Q2 [-]
75  end subroutine bc
76  end interface
77 
78  procedure(bc), pointer :: bulkflux => null()
79  public :: bulkflux
80 
81  !-----------------------------------------------------------------------------
82  !
83  !++ Public parameters & variables
84  !
85  character(len=H_SHORT), public :: bulkflux_type = 'B91W01' ! 'U95', 'B91', and 'B91W01'
86 
87  !-----------------------------------------------------------------------------
88  !
89  !++ Private procedure
90  !
91  private :: bulkflux_u95
92  private :: bulkflux_b91w01
93  private :: fm_unstable
94  private :: fh_unstable
95  private :: fm_stable
96  private :: fh_stable
97 
98  !-----------------------------------------------------------------------------
99  !
100  !++ Private parameters & variables
101  !
102  real(dp), parameter, private :: eps = 1.e-16_dp
103 
104  logical, private :: bulkflux_nk2018 = .true.
105 
106  integer, private :: bulkflux_itr_sa_max = 5 ! maximum iteration number for successive approximation
107  integer, private :: bulkflux_itr_nr_max = 10 ! maximum iteration number for Newton-Raphson method
108 
109  real(rp), private :: bulkflux_wscf ! empirical scaling factor of Wstar (Beljaars 1994)
110 
111  ! limiter
112  real(rp), private :: bulkflux_uabs_min = 1.0e-2_rp ! minimum of Uabs [m/s]
113  real(rp), private :: bulkflux_wstar_min = 1.0e-4_rp ! minimum of W* [m/s]
114 
115  ! surface diagnose
116  logical, private :: bulkflux_surfdiag_neutral = .true. ! calculate surface diagnoses with neutral condition
117 
118  logical, private :: flag_w01
119 
120 contains
121 
122  !-----------------------------------------------------------------------------
123  !
124  !-----------------------------------------------------------------------------
125  subroutine bulkflux_setup( dx )
126  use scale_prc, only: &
127  prc_abort
128  implicit none
129  real(rp), intent(in) :: dx
130 
131  namelist / param_bulkflux / &
132  bulkflux_type, &
133  bulkflux_nk2018, &
134  bulkflux_itr_sa_max, &
135  bulkflux_itr_nr_max, &
136  bulkflux_wscf, &
137  bulkflux_uabs_min, &
138  bulkflux_wstar_min, &
139  bulkflux_surfdiag_neutral
140 
141  integer :: ierr
142  !---------------------------------------------------------------------------
143 
144  log_newline
145  log_info("BULKFLUX_setup",*) 'Setup'
146 
147  ! WSCF = 1.2 for dx > 1 km (Beljaars 1994)
148  ! lim_{dx->0} WSCF = 0 for LES (Fig. 6 Kitamura and Ito 2016 BLM)
149  bulkflux_wscf = 1.2_rp * min(dx*1.e-3_rp, 1.0_rp)
150 
151 
152  !--- read namelist
153  rewind(io_fid_conf)
154  read(io_fid_conf,nml=param_bulkflux,iostat=ierr)
155 
156  if( ierr < 0 ) then !--- missing
157  log_info("BULKFLUX_setup",*) 'Not found namelist. Default used.'
158  elseif( ierr > 0 ) then !--- fatal error
159  log_error("BULKFLUX_setup",*) 'Not appropriate names in namelist PARAM_BULKFLUX. Check!'
160  call prc_abort
161  endif
162  log_nml(param_bulkflux)
163 
164  log_newline
165  log_info("BULKFLUX_setup",*) 'Scheme for surface bulk flux : ', trim(bulkflux_type)
166  select case(bulkflux_type)
167  case('U95')
168  log_info_cont(*) '=> Uno et al.(1995)'
169  bulkflux => bulkflux_u95
170  case('B91W01')
171  log_info_cont(*) '=> Beljaars and Holtslag (1991) and Wilson (2001)'
172  flag_w01 = .true.
173  bulkflux => bulkflux_b91w01
174  case('B91')
175  log_info_cont(*) '=> Beljaars and Holtslag (1991)'
176  flag_w01 = .false.
177  bulkflux => bulkflux_b91w01
178  case default
179  log_error("BULKFLUX_setup",*) 'Unsupported BULKFLUX_type. STOP'
180  call prc_abort
181  end select
182 
183  return
184  end subroutine bulkflux_setup
185 
186  !-----------------------------------------------------------------------------
187  ! estimate scales from fluxes
188  !-----------------------------------------------------------------------------
189  subroutine bulkflux_diagnose_scales_2d( &
190  IA, IS, IE, JA, JS, JE, &
191  SFLX_MW, SFLX_MU, SFLX_MV, &
192  SFLX_SH, SFLX_QV, &
193  SFC_DENS, SFC_TEMP, PBL, &
194  Ustar, Tstar, Qstar, &
195  Wstar, RLmo, &
196  mask )
197  use scale_const, only: &
198  eps => const_eps, &
199  grav => const_grav, &
200  cpdry => const_cpdry, &
201  karman => const_karman, &
202  epstvap => const_epstvap
203  implicit none
204  integer, intent(in) :: IA, IS, IE
205  integer, intent(in) :: JA, JS, JE
206 
207  real(RP), intent(in) :: SFLX_MW (IA,JA)
208  real(RP), intent(in) :: SFLX_MU (IA,JA)
209  real(RP), intent(in) :: SFLX_MV (IA,JA)
210  real(RP), intent(in) :: SFLX_SH (IA,JA)
211  real(RP), intent(in) :: SFLX_QV (IA,JA)
212  real(RP), intent(in) :: SFC_DENS(IA,JA)
213  real(RP), intent(in) :: SFC_TEMP(IA,JA)
214  real(RP), intent(in) :: PBL (IA,JA)
215 
216  real(RP), intent(out) :: Ustar(IA,JA)
217  real(RP), intent(out) :: Tstar(IA,JA)
218  real(RP), intent(out) :: Qstar(IA,JA)
219  real(RP), intent(out) :: Wstar(IA,JA)
220  real(RP), intent(out) :: RLmo (IA,JA)
221 
222  logical, intent(in), optional :: mask(IA,JA)
223 
224  integer :: i, j
225 
226  if ( present(mask) ) then
227  !$omp parallel do
228  do j = js, je
229  do i = is, ie
230  if ( mask(i,j) ) then
231  call bulkflux_diagnose_scales_0d( sflx_mw(i,j), sflx_mu(i,j), sflx_mv(i,j), & ! (in)
232  sflx_sh(i,j), sflx_qv(i,j), & ! (in)
233  sfc_dens(i,j), sfc_temp(i,j), pbl(i,j), & ! (in)
234  ustar(i,j), tstar(i,j), qstar(i,j), & ! (out)
235  wstar(i,j), rlmo(i,j) ) ! (out)
236  end if
237  end do
238  end do
239  else
240  !$omp parallel do
241  do j = js, je
242  do i = is, ie
243  call bulkflux_diagnose_scales_0d( sflx_mw(i,j), sflx_mu(i,j), sflx_mv(i,j), & ! (in)
244  sflx_sh(i,j), sflx_qv(i,j), & ! (in)
245  sfc_dens(i,j), sfc_temp(i,j), pbl(i,j), & ! (in)
246  ustar(i,j), tstar(i,j), qstar(i,j), & ! (out)
247  wstar(i,j), rlmo(i,j) ) ! (out)
248  end do
249  end do
250  end if
251 
252  return
253  end subroutine bulkflux_diagnose_scales_2d
254  !-----------------------------------------------------------------------------
255  subroutine bulkflux_diagnose_scales_0d( &
256  SFLX_MW, SFLX_MU, SFLX_MV, &
257  SFLX_SH, SFLX_QV, &
258  SFC_DENS, SFC_TEMP, PBL, &
259  Ustar, Tstar, Qstar, &
260  Wstar, RLmo )
261  use scale_const, only: &
262  eps => const_eps, &
263  grav => const_grav, &
264  cpdry => const_cpdry, &
265  karman => const_karman, &
266  epstvap => const_epstvap
267  implicit none
268  real(RP), intent(in) :: SFLX_MW
269  real(RP), intent(in) :: SFLX_MU
270  real(RP), intent(in) :: SFLX_MV
271  real(RP), intent(in) :: SFLX_SH
272  real(RP), intent(in) :: SFLX_QV
273  real(RP), intent(in) :: SFC_DENS
274  real(RP), intent(in) :: SFC_TEMP
275  real(RP), intent(in) :: PBL
276 
277  real(RP), intent(out) :: Ustar
278  real(RP), intent(out) :: Tstar
279  real(RP), intent(out) :: Qstar
280  real(RP), intent(out) :: Wstar
281  real(RP), intent(out) :: RLmo
282 
283  real(RP) :: BFLX, tmp, sw
284  logical :: ws_flag
285  integer :: i, j
286 
287  ws_flag = ( bulkflux_type == 'B91W01' )
288 
289  ustar = sqrt( sqrt( sflx_mw**2 + sflx_mu**2 + sflx_mv**2 ) / sfc_dens )
290  sw = 0.5_rp - sign( 0.5_rp, ustar - eps )
291  tstar = - sflx_sh / ( sfc_dens * ustar * cpdry + sw ) * ( 1.0_rp - sw )
292  qstar = - sflx_qv / ( sfc_dens * ustar + sw ) * ( 1.0_rp - sw )
293  bflx = - ustar * tstar - epstvap * ustar * qstar * sfc_temp
294  rlmo = - karman * grav * bflx / ( ustar**3 * sfc_temp + sw ) * ( 1.0_rp - sw )
295  if ( ws_flag ) then
296  tmp = pbl * grav / sfc_temp * bflx
297  sw = 0.5_rp + sign( 0.5_rp, tmp ) ! if tmp is plus, sw = 1
298  wstar = ( tmp * sw )**( 1.0_rp / 3.0_rp )
299  else
300  wstar = 0.0_rp
301  end if
302 
303  return
304  end subroutine bulkflux_diagnose_scales_0d
305 
306  !-----------------------------------------------------------------------------
307  ! estimate surface diagnose values (U10, V10, T2, Q2)
308  !-----------------------------------------------------------------------------
309  subroutine bulkflux_diagnose_surface_2d( &
310  IA, IS, IE, JA, JS, JE, &
311  ATM_U, ATM_V, &
312  ATM_TEMP, ATM_QV, &
313  SFC_TEMP, SFC_QV, &
314  ATM_Z1, &
315  SFC_Z0M, SFC_Z0H, SFC_Z0E, &
316  U10, V10, T2, Q2, &
317  mask, &
318  FracU10, FracT2, FracQ2 )
319  implicit none
320  integer, intent(in) :: IA, IS, IE
321  integer, intent(in) :: JA, JS, JE
322 
323  real(RP), intent(in) :: ATM_U (IA,JA)
324  real(RP), intent(in) :: ATM_V (IA,JA)
325  real(RP), intent(in) :: ATM_TEMP(IA,JA)
326  real(RP), intent(in) :: ATM_QV (IA,JA)
327  real(RP), intent(in) :: SFC_TEMP(IA,JA)
328  real(RP), intent(in) :: SFC_QV (IA,JA)
329  real(RP), intent(in) :: ATM_Z1 (IA,JA)
330  real(RP), intent(in) :: SFC_Z0M (IA,JA)
331  real(RP), intent(in) :: SFC_Z0H (IA,JA)
332  real(RP), intent(in) :: SFC_Z0E (IA,JA)
333 
334  real(RP), intent(out) :: U10(IA,JA)
335  real(RP), intent(out) :: V10(IA,JA)
336  real(RP), intent(out) :: T2 (IA,JA)
337  real(RP), intent(out) :: Q2 (IA,JA)
338 
339  logical, intent(in), optional :: mask (IA,JA)
340  real(RP), intent(in), optional :: FracU10(IA,JA)
341  real(RP), intent(in), optional :: FracT2 (IA,JA)
342  real(RP), intent(in), optional :: FracQ2 (IA,JA)
343 
344  integer :: i, j
345 
346  if ( present(fracu10) ) then
347  if ( present(mask) ) then
348  !$omp parallel do
349  do j = js, je
350  do i = is, ie
351  if ( mask(i,j) ) then
352  call bulkflux_diagnose_surface_0d( atm_u(i,j), atm_v(i,j), &
353  atm_temp(i,j), atm_qv(i,j), &
354  sfc_temp(i,j), sfc_qv(i,j), &
355  atm_z1(i,j), &
356  sfc_z0m(i,j), sfc_z0h(i,j), sfc_z0e(i,j), &
357  u10(i,j), v10(i,j), t2(i,j), q2(i,j), &
358  fracu10(i,j), fract2(i,j), fracq2(i,j) )
359  end if
360  end do
361  end do
362 
363  else
364  !$omp parallel do
365  do j = js, je
366  do i = is, ie
367  call bulkflux_diagnose_surface_0d( atm_u(i,j), atm_v(i,j), &
368  atm_temp(i,j), atm_qv(i,j), &
369  sfc_temp(i,j), sfc_qv(i,j), &
370  atm_z1(i,j), &
371  sfc_z0m(i,j), sfc_z0h(i,j), sfc_z0e(i,j), &
372  u10(i,j), v10(i,j), t2(i,j), q2(i,j), &
373  fracu10(i,j), fract2(i,j), fracq2(i,j) )
374  end do
375  end do
376  end if
377  else
378  if ( present(mask) ) then
379  !$omp parallel do
380  do j = js, je
381  do i = is, ie
382  if ( mask(i,j) ) then
383  call bulkflux_diagnose_surface_0d( atm_u(i,j), atm_v(i,j), &
384  atm_temp(i,j), atm_qv(i,j), &
385  sfc_temp(i,j), sfc_qv(i,j), &
386  atm_z1(i,j), &
387  sfc_z0m(i,j), sfc_z0h(i,j), sfc_z0e(i,j), &
388  u10(i,j), v10(i,j), t2(i,j), q2(i,j) )
389  end if
390  end do
391  end do
392 
393  else
394  !$omp parallel do
395  do j = js, je
396  do i = is, ie
397  call bulkflux_diagnose_surface_0d( atm_u(i,j), atm_v(i,j), &
398  atm_temp(i,j), atm_qv(i,j), &
399  sfc_temp(i,j), sfc_qv(i,j), &
400  atm_z1(i,j), &
401  sfc_z0m(i,j), sfc_z0h(i,j), sfc_z0e(i,j), &
402  u10(i,j), v10(i,j), t2(i,j), q2(i,j) )
403  end do
404  end do
405  end if
406  end if
407 
408  return
409  end subroutine bulkflux_diagnose_surface_2d
410 
411  subroutine bulkflux_diagnose_surface_0d( &
412  ATM_U, ATM_V, &
413  ATM_TEMP, ATM_QV, &
414  SFC_TEMP, SFC_QV, &
415  ATM_Z1, &
416  SFC_Z0M, SFC_Z0H, SFC_Z0E, &
417  U10, V10, T2, Q2, &
418  FracU10, FracT2, FracQ2 )
419  implicit none
420  real(RP), intent(in) :: ATM_U
421  real(RP), intent(in) :: ATM_V
422  real(RP), intent(in) :: ATM_TEMP
423  real(RP), intent(in) :: ATM_QV
424  real(RP), intent(in) :: SFC_TEMP
425  real(RP), intent(in) :: SFC_QV
426  real(RP), intent(in) :: ATM_Z1
427  real(RP), intent(in) :: SFC_Z0M
428  real(RP), intent(in) :: SFC_Z0H
429  real(RP), intent(in) :: SFC_Z0E
430 
431  real(RP), intent(out) :: U10
432  real(RP), intent(out) :: V10
433  real(RP), intent(out) :: T2
434  real(RP), intent(out) :: Q2
435 
436  real(RP), intent(in), optional :: FracU10
437  real(RP), intent(in), optional :: FracT2
438  real(RP), intent(in), optional :: FracQ2
439 
440  real(RP) :: f10m, f2h, f2e
441 
442  if ( bulkflux_surfdiag_neutral .or. ( .not. present(fracu10) ) ) then
443  ! assume the neutral condition
444  f10m = log( 10.0_rp / sfc_z0m ) / log( atm_z1 / sfc_z0m )
445  f2h = log( 2.0_rp / sfc_z0h ) / log( atm_z1 / sfc_z0h )
446  f2e = log( 2.0_rp / sfc_z0e ) / log( atm_z1 / sfc_z0e )
447  else
448  ! take the static stability into account
449  f10m = fracu10
450  f2h = fract2
451  f2e = fracq2
452  end if
453  u10 = atm_u * f10m
454  v10 = atm_v * f10m
455  t2 = sfc_temp + ( atm_temp - sfc_temp ) * f2h
456  q2 = sfc_qv + ( atm_qv - sfc_qv ) * f2e
457 
458  return
459  end subroutine bulkflux_diagnose_surface_0d
460 
461  !-----------------------------------------------------------------------------
462  ! ref. Uno et al. (1995)
463  !-----------------------------------------------------------------------------
464  subroutine bulkflux_u95( &
465  T1, T0, &
466  P1, P0, &
467  Q1, Q0, &
468  Uabs, Z1, PBL, &
469  Z0M, Z0H, Z0E, &
470  Ustar, Tstar, Qstar, &
471  Wstar, RLmo, Ra, &
472  FracU10, FracT2, FracQ2 )
473  use scale_const, only: &
474  grav => const_grav, &
475  karman => const_karman, &
476  rdry => const_rdry, &
477  cpdry => const_cpdry, &
478  pre00 => const_pre00
479  implicit none
480 
481  ! parameter
482  real(RP), parameter :: tPrn = 0.74_rp ! turbulent Prandtl number (Businger et al. 1971)
483  real(RP), parameter :: LFb = 9.4_rp ! Louis factor b (Louis 1979)
484  real(RP), parameter :: LFbp = 4.7_rp ! Louis factor b' (Louis 1979)
485  real(RP), parameter :: LFdm = 7.4_rp ! Louis factor d for momemtum (Louis 1979)
486  real(RP), parameter :: LFdh = 5.3_rp ! Louis factor d for heat (Louis 1979)
487 
488  ! argument
489  real(RP), intent(in) :: T1 ! tempearature at the lowest atmospheric layer [K]
490  real(RP), intent(in) :: T0 ! skin temperature [K]
491  real(RP), intent(in) :: P1 ! pressure at the lowest atmospheric layer [Pa]
492  real(RP), intent(in) :: P0 ! surface pressure [Pa]
493  real(RP), intent(in) :: Q1 ! mixing ratio at the lowest atmospheric layer [kg/kg]
494  real(RP), intent(in) :: Q0 ! surface mixing ratio [kg/kg]
495  real(RP), intent(in) :: Uabs! absolute velocity at the lowest atmospheric layer [m/s]
496  real(RP), intent(in) :: Z1 ! height at the lowest atmospheric layer [m]
497  real(RP), intent(in) :: PBL ! the top of atmospheric mixing layer [m]
498  real(RP), intent(in) :: Z0M ! roughness length of momentum [m]
499  real(RP), intent(in) :: Z0H ! roughness length of heat [m]
500  real(RP), intent(in) :: Z0E ! roughness length of moisture [m]
501 
502  real(RP), intent(out) :: Ustar ! friction velocity [m/s]
503  real(RP), intent(out) :: Tstar ! friction temperature [K]
504  real(RP), intent(out) :: Qstar ! friction mixing rate [kg/kg]
505  real(RP), intent(out) :: Wstar ! free convection velocity scale [m/s]
506  real(RP), intent(out) :: RLmo ! inversed Obukhov length [1/m]
507  real(RP), intent(out) :: Ra ! Aerodynamic resistance (=1/Ce)
508  real(RP), intent(out) :: FracU10 ! calculation parameter for U10 [-]
509  real(RP), intent(out) :: FracT2 ! calculation parameter for T2 [-]
510  real(RP), intent(out) :: FracQ2 ! calculation parameter for Q2 [-]
511 
512  ! work
513  real(RP) :: UabsW
514  real(RP) :: RiB0, RiB ! bulk Richardson number [-]
515  real(RP) :: C0Z1, C010, C002 ! initial drag coefficient [-]
516  real(RP) :: CmZ1, ChZ1, CqZ1, fmZ1, fhZ1, t0thZ1, q0qeZ1
517  real(RP) :: Cm10, fm10
518  real(RP) :: Cm02, Ch02, Cq02, fm02, fh02, t0th02, q0qe02
519  real(RP) :: TH1, TH0
520  real(RP) :: logZ1Z0M, log10Z0m, log02Z0m
521  real(RP) :: logZ0MZ0E
522  real(RP) :: logZ0MZ0H
523  !---------------------------------------------------------------------------
524 
525  logz1z0m = log( z1 / z0m )
526  log10z0m = log( 10.0_rp / z0m )
527  log02z0m = log( 2.0_rp / z0m )
528 
529  logz0mz0e = max( log( z0m/z0e ), 1.0_rp )
530  logz0mz0h = max( log( z0m/z0h ), 1.0_rp )
531 
532  uabsw = max( uabs, bulkflux_uabs_min )
533  th1 = t1 * ( p0 / p1 )**( rdry / cpdry )
534  th0 = t0
535 
536  ! bulk Richardson number
537  rib0 = grav * z1 * ( th1 - th0 ) / ( th1 * uabsw**2 )
538 
539  c0z1 = ( karman / logz1z0m )**2
540  c010 = ( karman / log10z0m )**2
541  c002 = ( karman / log02z0m )**2
542 
543  rib = rib0
544 
545  if( rib0 > 0.0_rp ) then
546  ! stable condition
547  fmz1 = 1.0_rp / ( 1.0_rp + lfbp * rib )**2
548  fhz1 = 1.0_rp / ( 1.0_rp + lfbp * rib )**2
549  fm10 = 1.0_rp / ( 1.0_rp + lfbp * rib )**2
550  fm02 = 1.0_rp / ( 1.0_rp + lfbp * rib )**2
551  fh02 = 1.0_rp / ( 1.0_rp + lfbp * rib )**2
552  else
553  ! unstable condition
554  fmz1 = 1.0_rp - lfb * rib / ( 1.0_rp + lfb * lfdm * c0z1 * sqrt( z1 / z0m ) * sqrt( abs( rib ) ) )
555  fhz1 = 1.0_rp - lfb * rib / ( 1.0_rp + lfb * lfdh * c0z1 * sqrt( z1 / z0m ) * sqrt( abs( rib ) ) )
556  fm10 = 1.0_rp - lfb * rib / ( 1.0_rp + lfb * lfdm * c010 * sqrt( 10.0_rp / z0m ) * sqrt( abs( rib ) ) )
557  fm02 = 1.0_rp - lfb * rib / ( 1.0_rp + lfb * lfdm * c002 * sqrt( 2.0_rp / z0m ) * sqrt( abs( rib ) ) )
558  fh02 = 1.0_rp - lfb * rib / ( 1.0_rp + lfb * lfdh * c002 * sqrt( 2.0_rp / z0m ) * sqrt( abs( rib ) ) )
559  end if
560 
561  t0thz1 = 1.0_rp / ( 1.0_rp + logz0mz0h / logz1z0m / sqrt( fmz1 ) * fhz1 )
562  q0qez1 = 1.0_rp / ( 1.0_rp + logz0mz0e / logz1z0m / sqrt( fmz1 ) * fhz1 )
563  t0th02 = 1.0_rp / ( 1.0_rp + logz0mz0h / log02z0m / sqrt( fm02 ) * fh02 )
564  q0qe02 = 1.0_rp / ( 1.0_rp + logz0mz0e / log02z0m / sqrt( fm02 ) * fh02 )
565 
566  rib = rib * t0thz1
567 
568  if( rib0 > 0.0_rp ) then
569  ! stable condition
570  fmz1 = 1.0_rp / ( 1.0_rp + lfbp * rib )**2
571  fhz1 = 1.0_rp / ( 1.0_rp + lfbp * rib )**2
572  fm10 = 1.0_rp / ( 1.0_rp + lfbp * rib )**2
573  fm02 = 1.0_rp / ( 1.0_rp + lfbp * rib )**2
574  fh02 = 1.0_rp / ( 1.0_rp + lfbp * rib )**2
575  else
576  ! unstable condition
577  fmz1 = 1.0_rp - lfb * rib / ( 1.0_rp + lfb * lfdm * c0z1 * sqrt( z1 / z0m ) * sqrt( abs( rib ) ) )
578  fhz1 = 1.0_rp - lfb * rib / ( 1.0_rp + lfb * lfdh * c0z1 * sqrt( z1 / z0m ) * sqrt( abs( rib ) ) )
579  fm10 = 1.0_rp - lfb * rib / ( 1.0_rp + lfb * lfdm * c010 * sqrt( 10.0_rp / z0m ) * sqrt( abs( rib ) ) )
580  fm02 = 1.0_rp - lfb * rib / ( 1.0_rp + lfb * lfdm * c002 * sqrt( 2.0_rp / z0m ) * sqrt( abs( rib ) ) )
581  fh02 = 1.0_rp - lfb * rib / ( 1.0_rp + lfb * lfdh * c002 * sqrt( 2.0_rp / z0m ) * sqrt( abs( rib ) ) )
582  end if
583 
584  t0thz1 = 1.0_rp / ( 1.0_rp + logz0mz0h / logz1z0m / sqrt( fmz1 ) * fhz1 )
585  q0qez1 = 1.0_rp / ( 1.0_rp + logz0mz0e / logz1z0m / sqrt( fmz1 ) * fhz1 )
586  t0th02 = 1.0_rp / ( 1.0_rp + logz0mz0h / log02z0m / sqrt( fm02 ) * fh02 )
587  q0qe02 = 1.0_rp / ( 1.0_rp + logz0mz0e / log02z0m / sqrt( fm02 ) * fh02 )
588 
589  cmz1 = c0z1 * fmz1
590  chz1 = c0z1 * fhz1 * t0thz1 / tprn
591  cqz1 = c0z1 * fhz1 * q0qez1 / tprn
592  cm10 = c010 * fm10
593  cm02 = c002 * fm02
594  ch02 = c002 * fh02 * t0th02 / tprn
595  cq02 = c002 * fh02 * q0qe02 / tprn
596 
597  ustar = sqrt( cmz1 ) * uabsw
598  tstar = chz1 * uabsw / ustar * ( th1 - th0 )
599  qstar = cqz1 * uabsw / ustar * ( q1 - q0 )
600  wstar = 0.0_rp
601 
602  fracu10 = sqrt( cmz1 / cm10 )
603  fract2 = chz1 / ch02 * sqrt( cm02 / cmz1 )
604  fracq2 = cqz1 / cq02 * sqrt( cm02 / cmz1 )
605 
606  rlmo = rib / z1 * karman * chz1 / ( cmz1 * sqrt( cmz1 ) )
607 
608  ra = 1.0_rp / ( cqz1 * uabsw )
609 
610  return
611  end subroutine bulkflux_u95
612 
613  !-----------------------------------------------------------------------------
614  !
615  ! refs. Beljaars and Holtslag (1991) and Wilson (2001)
616  !
617  ! Iteration method: refs. JMA-NHM Description Note II, Mar 2008
618  !
619  !-----------------------------------------------------------------------------
620 !OCL SERIAL
621  subroutine bulkflux_b91w01( &
622  T1, T0, &
623  P1, P0, &
624  Q1, Q0, &
625  Uabs, Z1, PBL, &
626  Z0M, Z0H, Z0E, &
627  Ustar, Tstar, Qstar, &
628  Wstar, RLmo, Ra, &
629  FracU10, FracT2, FracQ2 )
630  use scale_const, only: &
631  grav => const_grav, &
632  karman => const_karman, &
633  rdry => const_rdry, &
634  rvap => const_rvap, &
635  cpdry => const_cpdry, &
636  cpvap => const_cpvap, &
637  epstvap => const_epstvap, &
638  pre00 => const_pre00
639  implicit none
640 
641  ! parameter
642  real(DP), parameter :: dIL = 1.0e-6_dp ! delta [1/m]
643 
644  ! argument
645  real(RP), intent(in) :: T1 ! tempearature at the lowest atmospheric layer [K]
646  real(RP), intent(in) :: T0 ! skin temperature [K]
647  real(RP), intent(in) :: P1 ! pressure at the lowest atmospheric layer [Pa]
648  real(RP), intent(in) :: P0 ! surface pressure [Pa]
649  real(RP), intent(in) :: Q1 ! mixing ratio at the lowest atmospheric layer [kg/kg]
650  real(RP), intent(in) :: Q0 ! mixing ratio at surface [kg/kg]
651  real(RP), intent(in) :: Uabs! absolute velocity at the lowest atmospheric layer [m/s]
652  real(RP), intent(in) :: Z1 ! height at the lowest atmospheric layer [m]
653  real(RP), intent(in) :: PBL ! the top of atmospheric mixing layer [m]
654  real(RP), intent(in) :: Z0M ! roughness length of momentum [m]
655  real(RP), intent(in) :: Z0H ! roughness length of heat [m]
656  real(RP), intent(in) :: Z0E ! roughness length of moisture [m]
657 
658  real(RP), intent(out) :: Ustar ! friction velocity [m/s]
659  real(RP), intent(out) :: Tstar ! friction temperature [K]
660  real(RP), intent(out) :: Qstar ! friction mixing rate [kg/kg]
661  real(RP), intent(out) :: Wstar ! free convection velocity scale [m/s]
662  real(RP), intent(out) :: RLmo ! inversed Obukhov length [1/m]
663  real(RP), intent(out) :: Ra ! Aerodynamic resistance (=1/Ce)
664  real(RP), intent(out) :: FracU10 ! calculation parameter for U10 [-]
665  real(RP), intent(out) :: FracT2 ! calculation parameter for T2 [-]
666  real(RP), intent(out) :: FracQ2 ! calculation parameter for Q2 [-]
667 
668  ! work
669  integer :: n
670 
671  real(DP) :: res ! residual
672  real(DP) :: dres ! d(residual)/dIL
673 
674  real(DP) :: RiB0 ! bulk Richardson number [no unit]
675 
676  real(DP) :: UabsC
677  real(DP) :: UstarC, dUstarC
678  real(DP) :: TstarC, dTstarC
679  real(DP) :: QstarC, dQstarC
680  real(DP) :: WstarC, dWstar
681 
682  real(DP) :: Rtot, CPtot
683  real(DP) :: TH1, TH0
684  real(DP) :: TV1, TV0
685  real(DP) :: sw
686 
687  real(DP) :: IL, IL2
688  real(DP) :: BFLX, dBFLX
689 
690  real(DP) :: DP_Z1, DP_Z0M, DP_Z0H, DP_Z0E
691  real(DP) :: log_Z1ovZ0M, log_Z1ovZ0H, log_Z1ovZ0E
692 
693  real(DP) :: RzM, RzH, RzE
694  !---------------------------------------------------------------------------
695 
696  ! convert to DP
697  if ( bulkflux_nk2018 ) then
698  dp_z1 = real( z1*2.0_rp, kind=dp )
699  else
700  dp_z1 = real( z1, kind=dp )
701  end if
702  dp_z0m = real( z0m, kind=dp )
703  dp_z0h = real( z0h, kind=dp )
704  dp_z0e = real( z0e, kind=dp )
705 
706  uabsc = max( uabs, bulkflux_uabs_min )
707 
708  rtot = rdry * ( 1.0_dp - q1 ) + rvap * q1
709  cptot = cpdry * ( 1.0_dp - q1 ) + cpvap * q1
710  th1 = t1 * ( p0 / p1 )**( rtot / cptot )
711  th0 = t0
712  tv1 = th1 * ( 1.0_dp + epstvap * q1 )
713  tv0 = th0 * ( 1.0_dp + epstvap * q0 )
714 
715  rzm = 1.0_dp - dp_z0m / dp_z1
716  rzh = 1.0_dp - dp_z0h / dp_z1
717  rze = 1.0_dp - dp_z0e / dp_z1
718 
719  ! make log constant
720  log_z1ovz0m = log( dp_z1 / dp_z0m )
721  log_z1ovz0h = log( dp_z1 / dp_z0h )
722  log_z1ovz0e = log( dp_z1 / dp_z0e )
723 
724 
725  ! bulk Richardson number at initial step
726  rib0 = grav * dp_z1 * ( tv1 - tv0 ) / ( tv0 * uabsc**2 )
727 
728  ! first guess of inversed Obukhov length under neutral condition
729  il = rib0 / dp_z1 * log_z1ovz0m**2 / log_z1ovz0h
730 
731  ! free convection velocity scale at initial step
732  wstarc = bulkflux_wstar_min
733 
734  ! Successive approximation
735  do n = 1, bulkflux_itr_sa_max
736 
737  call calc_scales_b91w01( &
738  il, uabsc, th1, th0, tv0, q1, q0, pbl, & ! (in)
739  log_z1ovz0m, log_z1ovz0h, log_z1ovz0e, & ! (in)
740  dp_z1, dp_z0m, dp_z0h, dp_z0e, & ! (in)
741  rzm, rzh, rze, & ! (in)
742  wstarc, & ! (inout)
743  ustarc, tstarc, qstarc, bflx ) ! (out)
744 
745  ! estimate the inversed Obukhov length
746  il = - karman * grav * bflx / ( ustarc**3 * tv0 )
747  end do
748 
749  ! Newton-Raphson method
750  do n = 1, bulkflux_itr_nr_max
751 
752  dwstar = wstarc
753 
754  call calc_scales_b91w01( &
755  il, uabsc, th1, th0, tv0, q1, q0, pbl, & ! (in)
756  log_z1ovz0m, log_z1ovz0h, log_z1ovz0e, & ! (in)
757  dp_z1, dp_z0m, dp_z0h, dp_z0e, & ! (in)
758  rzm, rzh, rze, & ! (in)
759  wstarc, & ! (inout)
760  ustarc, tstarc, qstarc, bflx ) ! (out)
761 
762  il2 = sign( abs(il) + dil, il )
763  call calc_scales_b91w01( &
764  il2, uabsc, th1, th0, tv0, q1, q0, pbl, & ! (in)
765  log_z1ovz0m, log_z1ovz0h, log_z1ovz0e, & ! (in)
766  dp_z1, dp_z0m, dp_z0h, dp_z0e, & ! (in)
767  rzm, rzh, rze, & ! (in)
768  dwstar, & ! (inout)
769  dustarc, dtstarc, dqstarc, dbflx ) ! (out)
770 
771  res = il + karman * grav * bflx / ( ustarc**3 * tv0 )
772 
773  ! calculate d(residual)/dIL
774  dres = 1.0_dp + karman * grav / ( tv0 * sign(dil,il) ) * ( dbflx / dustarc**3 - bflx / ustarc**3 )
775 
776  ! stop iteration to prevent numerical error
777  if( abs( dres ) < 1e-20_rp ) exit
778 
779  ! convergence test with error levels
780  if( abs( res/dres ) < abs(il * eps) ) exit
781 
782  ! avoid sign changing
783  if( il * ( il - res / dres ) < 0.0_rp ) then
784  if ( abs(il) <= eps ) exit
785  il = sign(eps, il)
786  end if
787 
788  ! update the inversed Obukhov length
789  il = il - res / dres
790  end do
791 
792  ! Successive approximation after Newton-Raphson method
793  if( .NOT. abs( res/dres ) < abs(il * eps) ) then
794 
795  call calc_scales_b91w01( &
796  il, uabsc, th1, th0, tv0, q1, q0, pbl, & ! (in)
797  log_z1ovz0m, log_z1ovz0h, log_z1ovz0e, & ! (in)
798  dp_z1, dp_z0m, dp_z0h, dp_z0e, & ! (in)
799  rzm, rzh, rze, & ! (in)
800  wstarc, & ! (inout)
801  ustarc, tstarc, qstarc, bflx ) ! (out)
802 
803  ! estimate the inversed Obukhov length
804  il = - karman * grav * bflx / ( ustarc**3 * tv0 )
805  end if
806 
807  ! calculate Ustar, Tstar, and Qstar based on IL
808 
809  call calc_scales_b91w01( &
810  il, uabsc, th1, th0, tv0, q1, q0, pbl, & ! (in)
811  log_z1ovz0m, log_z1ovz0h, log_z1ovz0e, & ! (in)
812  dp_z1, dp_z0m, dp_z0h, dp_z0e, & ! (in)
813  rzm, rzh, rze, & ! (in)
814  wstarc, & ! (inout)
815  ustarc, tstarc, qstarc, bflx, & ! (out)
816  fracu10 = fracu10, & ! (out)
817  fract2 = fract2, fracq2 = fracq2 ) ! (out)
818 
819 
820  ! revert to RP
821  ustar = real( ustarc, kind=rp )
822  tstar = real( tstarc, kind=rp )
823  qstar = real( qstarc, kind=rp )
824  wstar = real( wstarc, kind=rp )
825 
826  ra = max( ( q1 - q0 ) / real(ustarc * qstarc + eps,kind=rp), real(eps,kind=rp) )
827 
828  rlmo = real(il, kind=rp)
829 
830  return
831  end subroutine bulkflux_b91w01
832 
833  subroutine calc_scales_b91w01( &
834  IL, Uabs, TH1, TH0, TV0, Q1, Q0, PBL, &
835  log_Z1ovZ0M, log_Z1ovZ0H, log_Z1ovZ0E, &
836  DP_Z1, DP_Z0M, DP_Z0H, DP_Z0E, &
837  RzM, RzH, RzE, &
838  Wstar, &
839  Ustar, Tstar, Qstar, BFLX, &
840  FracU10, FracT2, FracQ2 )
841  use scale_const, only: &
842  grav => const_grav, &
843  karman => const_karman, &
844  epstvap => const_epstvap
845  implicit none
846  real(DP), intent(in) :: IL
847  real(DP), intent(in) :: Uabs
848  real(DP), intent(in) :: TH1
849  real(DP), intent(in) :: TH0
850  real(DP), intent(in) :: TV0
851  real(RP), intent(in) :: Q1
852  real(RP), intent(in) :: Q0
853  real(RP), intent(in) :: PBL
854  real(DP), intent(in) :: log_Z1ovZ0M
855  real(DP), intent(in) :: log_Z1ovZ0H
856  real(DP), intent(in) :: log_Z1ovZ0E
857  real(DP), intent(in) :: DP_Z1
858  real(DP), intent(in) :: DP_Z0M
859  real(DP), intent(in) :: DP_Z0H
860  real(DP), intent(in) :: DP_Z0E
861  real(DP), intent(in) :: RzM
862  real(DP), intent(in) :: RzH
863  real(DP), intent(in) :: RzE
864 
865  real(DP), intent(inout) :: Wstar
866 
867  real(DP), intent(out) :: Ustar
868  real(DP), intent(out) :: Tstar
869  real(DP), intent(out) :: Qstar
870  real(DP), intent(out) :: BFLX
871 
872  real(RP), intent(out), optional :: FracU10
873  real(RP), intent(out), optional :: FracT2
874  real(RP), intent(out), optional :: FracQ2
875 
876  real(DP) :: denoM, denoH, denoE
877 
878  real(DP) :: tmp, sw
879 
880  if ( il < 0.0_dp ) then ! unstable condition
881 
882  if ( bulkflux_nk2018 ) then
883  denom = log_z1ovz0m &
884  - fmm_unstable(dp_z1,il) + fmm_unstable(dp_z0m,il) * dp_z0m / dp_z1 &
885  + rzm * ( fm_unstable(dp_z0m,il) - 1.0_dp )
886  tmp = fhm_unstable(dp_z1,il)
887  denoh = log_z1ovz0h &
888  - tmp + fhm_unstable(dp_z0h,il) * dp_z0h / dp_z1 &
889  + rzh * ( fh_unstable(dp_z0h,il) - 1.0_dp )
890  denoe = log_z1ovz0e &
891  - tmp + fhm_unstable(dp_z0e,il) * dp_z0e / dp_z1 &
892  + rze * ( fh_unstable(dp_z0e,il) - 1.0_dp )
893  else
894  denom = log_z1ovz0m - fm_unstable(dp_z1,il) + fm_unstable(dp_z0m,il)
895  tmp = fh_unstable(dp_z1,il)
896  denoh = log_z1ovz0h - tmp + fh_unstable(dp_z0h,il)
897  denoe = log_z1ovz0e - tmp + fh_unstable(dp_z0e,il)
898  end if
899  ustar = karman / denom * sqrt( uabs**2 + (bulkflux_wscf*wstar)**2 )
900  tstar = karman / denoh * ( th1 - th0 )
901  qstar = karman / denoe * ( q1 - q0 )
902  if ( present(fracu10) ) then
903  fracu10 = ( log(10.0_dp/dp_z0m) - fm_unstable(10.0_dp,il) + fm_unstable(dp_z0m,il) ) / denom
904  tmp = fm_unstable(2.0_dp,il)
905  fract2 = ( log(2.0_dp/dp_z0h) - tmp + fh_unstable(dp_z0h,il) ) / denoh
906  fracq2 = ( log(2.0_dp/dp_z0e) - tmp + fh_unstable(dp_z0e,il) ) / denoe
907  end if
908 
909  else ! stable condition
910 
911  if ( bulkflux_nk2018 ) then
912  denom = log_z1ovz0m &
913  - fmm_stable(dp_z1,il) + fmm_stable(dp_z0m,il) * dp_z0m / dp_z1 &
914  + rzm * ( fm_stable(dp_z0m,il) - 1.0_dp )
915  tmp = fhm_stable(dp_z1,il)
916  denoh = log_z1ovz0h &
917  - tmp + fhm_stable(dp_z0h,il) * dp_z0h / dp_z1 &
918  + rzh * ( fh_stable(dp_z0h,il) - 1.0_dp )
919  denoe = log_z1ovz0e &
920  - tmp + fhm_stable(dp_z0e,il) * dp_z0e / dp_z1 &
921  + rze * ( fh_stable(dp_z0e,il) - 1.0_dp )
922  else
923  denom = log_z1ovz0m - fm_stable(dp_z1,il) + fm_stable(dp_z0m,il)
924  tmp = fh_stable(dp_z1,il)
925  denoh = log_z1ovz0h - tmp + fh_stable(dp_z0h,il)
926  denoe = log_z1ovz0e - tmp + fh_stable(dp_z0e,il)
927  end if
928  ustar = karman / denom * uabs
929  tstar = karman / denoh * ( th1 - th0 )
930  qstar = karman / denoe * ( q1 - q0 )
931  if ( present(fracu10) ) then
932  fracu10 = ( log(10.0_dp/dp_z0m) - fm_stable(10.0_dp,il) + fm_stable(dp_z0m,il) ) / denom
933  tmp = fh_stable(2.0_dp,il)
934  fract2 = ( log(2.0_dp/dp_z0h) - tmp + fh_stable(dp_z0h,il) ) / denoh
935  fracq2 = ( log(2.0_dp/dp_z0e) - tmp + fh_stable(dp_z0e,il) ) / denoe
936  end if
937 
938  end if
939 
940  ! estimate buoyancy flux
941  bflx = - ustar * tstar - epstvap * ustar * qstar * tv0
942 
943  ! update free convection velocity scale
944  tmp = pbl * grav / tv0 * bflx
945  sw = 0.5_dp + sign( 0.5_dp, tmp ) ! if tmp is plus, sw = 1
946  wstar = ( tmp * sw )**( 1.0_dp / 3.0_dp )
947 
948  ! avoid zero division with UstarC = 0
949  ustar = sign( max( abs(ustar), eps ), ustar )
950 
951  return
952  end subroutine calc_scales_b91w01
953 
954  !-----------------------------------------------------------------------------
955  ! stability function for momemtum in unstable condition
956  function fm_unstable( Z, IL )
957  use scale_const, only: &
958  pi => const_pi
959  implicit none
960 
961  ! argument
962  real(DP), intent(in) :: Z
963  real(DP), intent(in) :: IL
964 
965  ! function
966  real(DP) :: fm_unstable
967 
968  ! Wilson (2001)
969  real(DP), parameter :: gamma = 3.6_dp
970 
971  ! works
972  real(DP) :: R
973  real(DP) :: r4R
974  !---------------------------------------------------------------------------
975 
976  r = min( z * il, 0.0_dp )
977 
978  if ( flag_w01 ) then
979  ! Wilson (2001)
980  fm_unstable = 3.0_dp * log( ( 1.0_dp + sqrt( 1.0_dp + gamma * (-r)**(2.0_dp/3.0_dp) ) ) * 0.5_dp )
981  else
982  ! Beljaars and Holtslag (1991), originally Paulson (1974) and Dyer (1974)
983  r4r = sqrt( sqrt( 1.0_dp - 16.0_dp * r ) )
984  fm_unstable = log( ( 1.0_dp + r4r )**2 * ( 1.0_dp + r4r * r4r ) * 0.125_dp ) - 2.0_dp * atan( r4r ) + pi * 0.5_dp
985  end if
986 
987  return
988  end function fm_unstable
989  function fmm_unstable( Z, IL )
990  use scale_const, only: &
991  pi => const_pi
992  implicit none
993 
994  ! argument
995  real(dp), intent(in) :: z
996  real(dp), intent(in) :: il
997 
998  ! function
999  real(dp) :: fmm_unstable
1000 
1001  ! Wilson (2001)
1002  real(dp), parameter :: gamma = 3.6_dp
1003 
1004  ! works
1005  real(dp) :: r
1006  real(dp) :: f, r3
1007  real(dp) :: r4r, r2r
1008  !---------------------------------------------------------------------------
1009 
1010  r = min( z * il, 0.0_dp )
1011 
1012  if ( flag_w01 ) then
1013  ! Wilson (2001)
1014  r3 = (-r)**(1.0_dp/3.0_dp)
1015  if ( r > -eps ) then
1016  fmm_unstable = 9.0_dp / 20.0_dp * gamma * r3**2
1017  else
1018  f = sqrt( 1.0_dp + gamma * r3**2 )
1019  fmm_unstable = 3.0_dp * log( ( 1.0_dp + f ) * 0.5_dp ) &
1020  + 1.5_dp / ( sqrt(gamma)**3 * r) * asinh( sqrt(gamma) * r3 ) &
1021  + 1.5_dp * f / ( gamma * r3**2 ) &
1022  - 1.0_dp
1023  end if
1024  else
1025  ! Beljaars and Holtslag (1991), originally Paulson (1974) and Dyer (1974)
1026  if ( r > -eps ) then ! |R|<EPS, now R < 0
1027  fmm_unstable = - 2.0_dp * r
1028  else
1029  r2r = sqrt( 1.0_dp - 16.0_dp * r )
1030  r4r = sqrt( r2r )
1031  fmm_unstable = log( ( 1.0_dp + r4r )**2 * ( 1.0_dp + r2r ) * 0.125_dp ) &
1032  - 2.0_dp * atan( r4r ) &
1033  + ( 1.0_dp - r4r*r2r ) / ( 12.0_dp * r ) &
1034  + pi * 0.5_dp - 1.0_dp
1035  end if
1036  end if
1037 
1038  return
1039  end function fmm_unstable
1040 
1041  !-----------------------------------------------------------------------------
1042  ! stability function for heat/vapor in unstable condition
1043  function fh_unstable( Z, IL )
1044  implicit none
1045 
1046  ! argument
1047  real(dp), intent(in) :: z
1048  real(dp), intent(in) :: il
1049 
1050  ! function
1051  real(dp) :: fh_unstable
1052 
1053  ! Wilson (2001)
1054  real(dp), parameter :: pt = 0.95_dp ! turbulent Prandtl number
1055  real(dp), parameter :: gamma = 7.9_dp
1056 
1057  ! works
1058  real(dp) :: r
1059  !---------------------------------------------------------------------------
1060 
1061  r = min( z * il, 0.0_dp )
1062 
1063  if ( flag_w01 ) then
1064  ! Wilson (2001)
1065  fh_unstable = 3.0_dp * log( ( 1.0_dp + sqrt( 1.0_dp + gamma * (-r)**(2.0_dp/3.0_dp) ) ) * 0.5_dp ) * pt &
1066  + log( z ) * ( 1.0_dp - pt )
1067  else
1068  ! Beljaars and Holtslag (1991), originally Paulson (1974); Dyer (1974)
1069  fh_unstable = 2.0_dp * log( ( 1.0_dp + sqrt( 1.0_dp - 16.0_dp * r ) ) * 0.5_dp )
1070  end if
1071 
1072  return
1073  end function fh_unstable
1074  function fhm_unstable( Z, IL )
1075  implicit none
1076 
1077  ! argument
1078  real(dp), intent(in) :: z
1079  real(dp), intent(in) :: il
1080 
1081  ! function
1082  real(dp) :: fhm_unstable
1083 
1084  ! Wilson (2001)
1085  real(dp), parameter :: pt = 0.95_dp ! turbulent Prandtl number
1086  real(dp), parameter :: gamma = 7.9_dp
1087 
1088  ! works
1089  real(dp) :: r
1090  real(dp) :: f, r3
1091  real(dp) :: r2r
1092  !---------------------------------------------------------------------------
1093 
1094  r = min( z * il, 0.0_dp )
1095 
1096  if ( flag_w01 ) then
1097  ! Wilson (2001)
1098  r3 = (-r)**(1.0_dp/3.0_dp)
1099  if ( r > -eps ) then
1100  fhm_unstable = 9.0_dp / 20.0_dp * gamma * r3**2
1101  else
1102  f = sqrt( 1.0_dp + gamma * r3**2 )
1103  fhm_unstable = 3.0_dp * log( ( 1.0_dp + f ) * 0.5_dp ) &
1104  + 1.5_dp / ( sqrt(gamma)**3 * r) * asinh( sqrt(gamma) * r3 ) &
1105  + 1.5_dp * f / ( gamma * r3**2 ) &
1106  - 1.0_dp
1107  fhm_unstable = fhm_unstable * pt + ( log( z ) - 1.0_dp ) * ( 1.0_dp - pt )
1108  end if
1109  else
1110  ! Beljaars and Holtslag (1991), originally Paulson (1974); Dyer (1974)
1111  if ( r > -eps ) then ! |R| < EPS, now R < 0
1112  fhm_unstable = - 4.0_dp * r
1113  else
1114  r2r = sqrt( 1.0_dp - 16.0_dp * r )
1115  fhm_unstable = 2.0_dp * log( ( 1.0_dp + r2r ) * 0.5_dp ) &
1116  + ( 1.0_dp - r2r ) / ( 8.0_dp * r ) &
1117  - 1.0_dp
1118  end if
1119  end if
1120 
1121  return
1122  end function fhm_unstable
1123 
1124  !-----------------------------------------------------------------------------
1125  ! stability function for momemtum in stable condition
1126  function fm_stable( Z, IL )
1127  implicit none
1128 
1129  ! argument
1130  real(dp), intent(in) :: z
1131  real(dp), intent(in) :: il
1132 
1133  ! function
1134  real(dp) :: fm_stable
1135 
1136  ! parameters of stability functions (Beljaars and Holtslag 1991)
1137  real(dp), parameter :: a = 1.0_dp
1138  real(dp), parameter :: b = 0.667_dp
1139  real(dp), parameter :: c = 5.0_dp
1140  real(dp), parameter :: d = 0.35_dp
1141 
1142  ! works
1143  real(dp) :: r
1144  !---------------------------------------------------------------------------
1145 
1146  r = max( z * il, 0.0_dp )
1147 
1148  ! Holtslag and DeBruin (1988)
1149 #if defined(PGI) || defined(SX)
1150  fm_stable = - a*r - b*( r - c/d )*exp( -min( d*r, 1.e+3_rp ) ) - b*c/d ! apply exp limiter
1151 #else
1152  fm_stable = - a*r - b*( r - c/d )*exp( -d*r ) - b*c/d
1153 #endif
1154 
1155  return
1156  end function fm_stable
1157  function fmm_stable( Z, IL )
1158  implicit none
1159 
1160  ! argument
1161  real(dp), intent(in) :: z
1162  real(dp), intent(in) :: il
1163 
1164  ! function
1165  real(dp) :: fmm_stable
1166 
1167  ! parameters of stability functions (Beljaars and Holtslag 1991)
1168  real(dp), parameter :: a = 1.0_dp
1169  real(dp), parameter :: b = 0.667_dp
1170  real(dp), parameter :: c = 5.0_dp
1171  real(dp), parameter :: d = 0.35_dp
1172 
1173  ! works
1174  real(dp) :: r
1175  !---------------------------------------------------------------------------
1176 
1177  r = max( z * il, 0.0_dp )
1178 
1179  ! Holtslag and DeBruin (1988)
1180  if ( r < eps ) then
1181  fmm_stable = - 0.5_dp * ( a + b * c + d ) * r
1182  else
1183  fmm_stable = b * ( d*r - c + 1.0_dp ) / ( d**2 * r ) &
1184 #if defined(PGI) || defined(SX)
1185  ! apply exp limiter
1186  * exp( -min( d*r, 1.e+3_dp ) ) &
1187 #else
1188  * exp( -d*r ) &
1189 #endif
1190  - a * r * 0.5_dp &
1191  - b * ( c*d*r - c + 1.0_dp ) / ( d**2 * r )
1192  end if
1193 
1194  return
1195  end function fmm_stable
1196 
1197  !-----------------------------------------------------------------------------
1198  ! stability function for heat/vapor in stable condition
1199  function fh_stable( Z, IL )
1200  implicit none
1201 
1202  ! argument
1203  real(dp), intent(in) :: z
1204  real(dp), intent(in) :: il
1205 
1206  ! function
1207  real(dp) :: fh_stable
1208 
1209  ! parameters of stability functions (Beljaars and Holtslag 1991)
1210  real(dp), parameter :: a = 1.0_dp
1211  real(dp), parameter :: b = 0.667_dp
1212  real(dp), parameter :: c = 5.0_dp
1213  real(dp), parameter :: d = 0.35_dp
1214 
1215  ! works
1216  real(dp) :: r
1217  !---------------------------------------------------------------------------
1218 
1219  r = max( z * il, 0.0_dp )
1220 
1221  ! Beljaars and Holtslag (1991)
1222 #if defined(PGI) || defined(SX)
1223  fh_stable = 1.0_dp - ( 1.0_dp + 2.0_dp/3.0_dp * a*r )**1.5_dp - b*( r - c/d )*exp( -min( d*r, 1.e+3_rp ) ) - b*c/d ! apply exp limiter
1224 #else
1225  fh_stable = 1.0_dp - ( 1.0_dp + 2.0_dp/3.0_dp * a*r )**1.5_dp - b*( r - c/d )*exp( -d*r ) - b*c/d
1226 #endif
1227 
1228  return
1229  end function fh_stable
1230  function fhm_stable( Z, IL )
1231  implicit none
1232 
1233  ! argument
1234  real(dp), intent(in) :: z
1235  real(dp), intent(in) :: il
1236 
1237  ! function
1238  real(dp) :: fhm_stable
1239 
1240  ! parameters of stability functions (Beljaars and Holtslag 1991)
1241  real(dp), parameter :: a = 1.0_dp
1242  real(dp), parameter :: b = 0.667_dp
1243  real(dp), parameter :: c = 5.0_dp
1244  real(dp), parameter :: d = 0.35_dp
1245 
1246  ! works
1247  real(dp) :: r
1248  !---------------------------------------------------------------------------
1249 
1250  r = max( z * il, 0.0_dp )
1251 
1252  ! Beljaars and Holtslag (1991)
1253  if ( r < eps ) then
1254  fhm_stable = - 0.5_dp * ( a + b*c + b ) * r
1255  else
1256  fhm_stable = b * ( d*r - c + 1.0_dp ) / ( d**2 * r ) &
1257 #if defined(PGI) || defined(SX)
1258  * exp( -min( d*r, 1.e+3_dp) ) &
1259 #else
1260  * exp( -d*r ) &
1261 #endif
1262  - 3.0_dp * sqrt( 1.0_dp + 2.0_dp*a*r/3.0_dp )**5 / ( 5.0_dp * a * r ) &
1263  - b * ( c*d*r - c + 1.0_dp ) / ( d**2 * r ) &
1264  + 0.6_rp / ( a * r ) + 1.0_dp
1265  end if
1266 
1267  return
1268  end function fhm_stable
1269 
1270 end module scale_bulkflux
scale_const::const_grav
real(rp), public const_grav
standard acceleration of gravity [m/s2]
Definition: scale_const.F90:46
scale_bulkflux::bulkflux_diagnose_scales_0d
subroutine bulkflux_diagnose_scales_0d(SFLX_MW, SFLX_MU, SFLX_MV, SFLX_SH, SFLX_QV, SFC_DENS, SFC_TEMP, PBL, Ustar, Tstar, Qstar, Wstar, RLmo)
Definition: scale_bulkflux.F90:261
scale_prc::prc_abort
subroutine, public prc_abort
Abort Process.
Definition: scale_prc.F90:342
scale_const::const_epstvap
real(rp), public const_epstvap
1 / epsilon - 1
Definition: scale_const.F90:70
scale_bulkflux::bulkflux_diagnose_surface_2d
subroutine bulkflux_diagnose_surface_2d(IA, IS, IE, JA, JS, JE, ATM_U, ATM_V, ATM_TEMP, ATM_QV, SFC_TEMP, SFC_QV, ATM_Z1, SFC_Z0M, SFC_Z0H, SFC_Z0E, U10, V10, T2, Q2, mask, FracU10, FracT2, FracQ2)
Definition: scale_bulkflux.F90:319
scale_precision
module PRECISION
Definition: scale_precision.F90:14
scale_const::const_rvap
real(rp), parameter, public const_rvap
specific gas constant (water vapor) [J/kg/K]
Definition: scale_const.F90:63
scale_const::const_eps
real(rp), public const_eps
small number
Definition: scale_const.F90:33
scale_bulkflux
module Surface bulk flux
Definition: scale_bulkflux.F90:12
scale_bulkflux::fmm_unstable
real(dp) function fmm_unstable(Z, IL)
Definition: scale_bulkflux.F90:990
scale_const::const_cpvap
real(rp), parameter, public const_cpvap
specific heat (water vapor, constant pressure) [J/kg/K]
Definition: scale_const.F90:64
scale_prc
module PROCESS
Definition: scale_prc.F90:11
scale_const::const_pi
real(rp), public const_pi
pi
Definition: scale_const.F90:31
scale_precision::rp
integer, parameter, public rp
Definition: scale_precision.F90:41
scale_io
module STDIO
Definition: scale_io.F90:10
scale_bulkflux::bulkflux_setup
subroutine, public bulkflux_setup(dx)
Definition: scale_bulkflux.F90:126
scale_const
module CONSTANT
Definition: scale_const.F90:11
scale_const::const_cpdry
real(rp), public const_cpdry
specific heat (dry air,constant pressure) [J/kg/K]
Definition: scale_const.F90:56
scale_prof
module profiler
Definition: scale_prof.F90:11
scale_precision::dp
integer, parameter, public dp
Definition: scale_precision.F90:32
scale_bulkflux::bulkflux_type
character(len=h_short), public bulkflux_type
Definition: scale_bulkflux.F90:85
scale_const::const_karman
real(rp), parameter, public const_karman
von Karman constant
Definition: scale_const.F90:50
scale_const::const_rdry
real(rp), public const_rdry
specific gas constant (dry air) [J/kg/K]
Definition: scale_const.F90:55
scale_bulkflux::bulkflux_diagnose_scales_2d
subroutine bulkflux_diagnose_scales_2d(IA, IS, IE, JA, JS, JE, SFLX_MW, SFLX_MU, SFLX_MV, SFLX_SH, SFLX_QV, SFC_DENS, SFC_TEMP, PBL, Ustar, Tstar, Qstar, Wstar, RLmo, mask)
Definition: scale_bulkflux.F90:197
scale_bulkflux::calc_scales_b91w01
subroutine calc_scales_b91w01(IL, Uabs, TH1, TH0, TV0, Q1, Q0, PBL, log_Z1ovZ0M, log_Z1ovZ0H, log_Z1ovZ0E, DP_Z1, DP_Z0M, DP_Z0H, DP_Z0E, RzM, RzH, RzE, Wstar, Ustar, Tstar, Qstar, BFLX, FracU10, FracT2, FracQ2)
Definition: scale_bulkflux.F90:841
scale_const::const_pre00
real(rp), public const_pre00
pressure reference [Pa]
Definition: scale_const.F90:88
scale_io::io_fid_conf
integer, public io_fid_conf
Config file ID.
Definition: scale_io.F90:56
scale_bulkflux::bulkflux
procedure(bc), pointer, public bulkflux
Definition: scale_bulkflux.F90:78
Generated by   doxygen 1.8.17