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