SCALE-RM
scale_bulkflux.F90
Go to the documentation of this file.
1 !-------------------------------------------------------------------------------
11 !-------------------------------------------------------------------------------
12 module scale_bulkflux
13  !-----------------------------------------------------------------------------
14  !
15  !++ used modules
16  !
17  use scale_precision
18  use scale_stdio
19  !-----------------------------------------------------------------------------
20  implicit none
21  private
22  !-----------------------------------------------------------------------------
23  !
24  !++ Public procedure
25  !
26  public :: bulkflux_setup
27 
28  abstract interface
29  subroutine bc( &
30  Ustar, & ! (out)
31  tstar, & ! (out)
32  qstar, & ! (out)
33  uabs, & ! (out)
34  t1, & ! (in)
35  t0, & ! (in)
36  p1, & ! (in)
37  p0, & ! (in)
38  q1, & ! (in)
39  q0, & ! (in)
40  u1, & ! (in)
41  v1, & ! (in)
42  z1, & ! (in)
43  pbl, & ! (in)
44  z0m, & ! (in)
45  z0h, & ! (in)
46  z0e ) ! (in)
47  use scale_precision
48  implicit none
49 
50  real(RP), intent(out) :: ustar ! friction velocity [m/s]
51  real(RP), intent(out) :: tstar ! friction temperature [K]
52  real(RP), intent(out) :: qstar ! friction mixing rate [kg/kg]
53  real(RP), intent(out) :: uabs ! modified absolute velocity [m/s]
54 
55  real(RP), intent(in) :: t1 ! tempearature at the lowest atmospheric layer [K]
56  real(RP), intent(in) :: t0 ! skin temperature [K]
57  real(RP), intent(in) :: p1 ! pressure at the lowest atmospheric layer [Pa]
58  real(RP), intent(in) :: p0 ! surface pressure [Pa]
59  real(RP), intent(in) :: q1 ! mixing ratio at the lowest atmospheric layer [kg/kg]
60  real(RP), intent(in) :: q0 ! surface mixing ratio [kg/kg]
61  real(RP), intent(in) :: u1 ! zonal wind at the lowest atmospheric layer [m/s]
62  real(RP), intent(in) :: v1 ! meridional wind 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  end subroutine bc
69  end interface
70 
71  procedure(bc), pointer :: bulkflux => null()
72  public :: bulkflux
73 
74  !-----------------------------------------------------------------------------
75  !
76  !++ Public parameters & variables
77  !
78  !-----------------------------------------------------------------------------
79  !
80  !++ Private procedure
81  !
82  private :: bulkflux_u95
83  private :: bulkflux_b91w01
84  private :: fm_unstable
85  private :: fh_unstable
86  private :: fm_stable
87  private :: fh_stable
88 
89  !-----------------------------------------------------------------------------
90  !
91  !++ Private parameters & variables
92  !
93  character(len=H_SHORT), private :: bulkflux_type = 'B91W01'
94 
95  integer, private :: bulkflux_itr_max = 100 ! maximum iteration number
96 
97  real(RP), private :: bulkflux_res_min = 1.0e+0_rp ! minimum value of residual
98  real(RP), private :: bulkflux_err_min = 1.0e-2_rp ! minimum value of error
99 
100  real(RP), private :: bulkflux_wscf = 1.2_rp ! empirical scaling factor of Wstar (Beljaars 1994)
101 
102  ! limiter
103  real(RP), private :: bulkflux_uabs_min = 1.0e-2_rp ! minimum of Uabs [m/s]
104  real(RP), private :: bulkflux_rib_min = 1.0e-4_rp ! minimum of RiB [no unit]
105  real(RP), private :: bulkflux_wstar_min = 1.0e-4_rp ! minimum of W* [m/s]
106 
107 contains
108 
109  !-----------------------------------------------------------------------------
110  !
111  !-----------------------------------------------------------------------------
112  subroutine bulkflux_setup
113  use scale_process, only: &
114  prc_mpistop
115  implicit none
116 
117  namelist / param_bulkflux / &
118  bulkflux_type, &
119  bulkflux_itr_max, &
120  bulkflux_res_min, &
121  bulkflux_err_min, &
122  bulkflux_wscf, &
123  bulkflux_uabs_min, &
124  bulkflux_rib_min, &
125  bulkflux_wstar_min
126 
127  integer :: ierr
128  !---------------------------------------------------------------------------
129 
130  if( io_l ) write(io_fid_log,*)
131  if( io_l ) write(io_fid_log,*) '++++++ Module[BULKFLUX] / Categ[COUPLER] / Origin[SCALElib]'
132 
133  !--- read namelist
134  rewind(io_fid_conf)
135  read(io_fid_conf,nml=param_bulkflux,iostat=ierr)
136 
137  if( ierr < 0 ) then !--- missing
138  if( io_l ) write(io_fid_log,*) '*** Not found namelist. Default used.'
139  elseif( ierr > 0 ) then !--- fatal error
140  write(*,*) 'xxx Not appropriate names in namelist PARAM_BULKFLUX. Check!'
141  call prc_mpistop
142  endif
143  if( io_lnml ) write(io_fid_log,nml=param_bulkflux)
144 
145  select case( bulkflux_type )
146  case ( 'U95' )
147  if( io_l ) write(io_fid_log,*) '*** Scheme for surface bulk flux : Uno et al.(1995)'
148  bulkflux => bulkflux_u95
149  case ( 'B91W01' )
150  if( io_l ) write(io_fid_log,*) '*** Scheme for surface bulk flux : Beljaars (1991) and Wilson (2001)'
151  bulkflux => bulkflux_b91w01
152  case default
153  write(*,*) ' xxx Unsupported TYPE. STOP'
154  call prc_mpistop
155  end select
156 
157  return
158  end subroutine bulkflux_setup
159 
160  !-----------------------------------------------------------------------------
161  ! ref. Uno et al. (1995)
162  !-----------------------------------------------------------------------------
163  subroutine bulkflux_u95( &
164  Ustar, & ! (out)
165  tstar, & ! (out)
166  qstar, & ! (out)
167  uabs, & ! (out)
168  t1, & ! (in)
169  t0, & ! (in)
170  p1, & ! (in)
171  p0, & ! (in)
172  q1, & ! (in)
173  q0, & ! (in)
174  u1, & ! (in)
175  v1, & ! (in)
176  z1, & ! (in)
177  pbl, & ! (in)
178  z0m, & ! (in)
179  z0h, & ! (in)
180  z0e ) ! (in)
181  use scale_const, only: &
182  grav => const_grav, &
183  karman => const_karman, &
184  rdry => const_rdry, &
185  cpdry => const_cpdry, &
186  pre00 => const_pre00
187  implicit none
188 
189  ! parameter
190  real(RP), parameter :: tPrn = 0.74_rp ! turbulent Prandtl number (Businger et al. 1971)
191  real(RP), parameter :: LFb = 9.4_rp ! Louis factor b (Louis 1979)
192  real(RP), parameter :: LFbp = 4.7_rp ! Louis factor b' (Louis 1979)
193  real(RP), parameter :: LFdm = 7.4_rp ! Louis factor d for momemtum (Louis 1979)
194  real(RP), parameter :: LFdh = 5.3_rp ! Louis factor d for heat (Louis 1979)
195 
196  ! argument
197  real(RP), intent(out) :: Ustar ! friction velocity [m/s]
198  real(RP), intent(out) :: Tstar ! friction temperature [K]
199  real(RP), intent(out) :: Qstar ! friction mixing rate [kg/kg]
200  real(RP), intent(out) :: Uabs ! modified absolute velocity [m/s]
201 
202  real(RP), intent(in) :: T1 ! tempearature at the lowest atmospheric layer [K]
203  real(RP), intent(in) :: T0 ! skin temperature [K]
204  real(RP), intent(in) :: P1 ! pressure at the lowest atmospheric layer [Pa]
205  real(RP), intent(in) :: P0 ! surface pressure [Pa]
206  real(RP), intent(in) :: Q1 ! mixing ratio at the lowest atmospheric layer [kg/kg]
207  real(RP), intent(in) :: Q0 ! surface mixing ratio [kg/kg]
208  real(RP), intent(in) :: U1 ! zonal wind at the lowest atmospheric layer [m/s]
209  real(RP), intent(in) :: V1 ! meridional wind at the lowest atmospheric layer [m/s]
210  real(RP), intent(in) :: Z1 ! height at the lowest atmospheric layer [m]
211  real(RP), intent(in) :: PBL ! the top of atmospheric mixing layer [m]
212  real(RP), intent(in) :: Z0M ! roughness length of momentum [m]
213  real(RP), intent(in) :: Z0H ! roughness length of heat [m]
214  real(RP), intent(in) :: Z0E ! roughness length of moisture [m]
215 
216  ! work
217  real(RP) :: RiB0, RiB ! bulk Richardson number [no unit]
218  real(RP) :: C0 ! initial drag coefficient [no unit]
219  real(RP) :: fm, fh, t0th, q0qe
220  real(RP) :: TH1, TH0
221  real(RP) :: logZ1Z0M
222  real(RP) :: logZ0MZ0E
223  real(RP) :: logZ0MZ0H
224  !---------------------------------------------------------------------------
225 
226  logz1z0m = log( z1/z0m )
227  logz0mz0e = max( log( z0m/z0e ), 1.0_rp )
228  logz0mz0h = max( log( z0m/z0h ), 1.0_rp )
229 
230  uabs = max( sqrt( u1**2 + v1**2 ), bulkflux_uabs_min )
231  th1 = t1 * ( pre00 / p1 )**( rdry / cpdry )
232  th0 = t0 * ( pre00 / p0 )**( rdry / cpdry )
233 
234  rib0 = grav * z1 * ( th1 - th0 ) / ( th1 * uabs**2 )
235  if( abs( rib0 ) < bulkflux_rib_min ) then
236  rib0 = sign( bulkflux_rib_min, rib0 )
237  end if
238 
239  c0 = ( karman / logz1z0m )**2
240  rib = rib0
241 
242  if( rib0 > 0.0_rp ) then
243  ! stable condition
244  fm = 1.0_rp / ( 1.0_rp + lfbp * rib )**2
245  fh = fm
246  else
247  ! unstable condition
248  fm = 1.0_rp - lfb * rib / ( 1.0_rp + lfb * lfdm * c0 * sqrt( z1/z0m ) * sqrt( abs( rib ) ) )
249  fh = 1.0_rp - lfb * rib / ( 1.0_rp + lfb * lfdh * c0 * sqrt( z1/z0m ) * sqrt( abs( rib ) ) )
250  end if
251 
252  t0th = 1.0_rp / ( 1.0_rp + logz0mz0h / logz1z0m / sqrt( fm ) * fh )
253  q0qe = 1.0_rp / ( 1.0_rp + logz0mz0e / logz1z0m / sqrt( fm ) * fh )
254  rib = rib * t0th
255 
256  if( rib0 > 0.0_rp ) then
257  ! stable condition
258  fm = 1.0_rp / ( 1.0_rp + lfbp * rib )**2
259  fh = fm
260  else
261  ! unstable condition
262  fm = 1.0_rp - lfb * rib / ( 1.0_rp + lfb * lfdm * c0 * sqrt( z1/z0m ) * sqrt( abs( rib ) ) )
263  fh = 1.0_rp - lfb * rib / ( 1.0_rp + lfb * lfdh * c0 * sqrt( z1/z0m ) * sqrt( abs( rib ) ) )
264  end if
265 
266  t0th = 1.0_rp / ( 1.0_rp + logz0mz0h / logz1z0m / sqrt( fm ) * fh )
267  q0qe = 1.0_rp / ( 1.0_rp + logz0mz0e / logz1z0m / sqrt( fm ) * fh )
268 
269  ustar = sqrt( c0 * fm ) * uabs
270  tstar = c0 * fh * t0th / tprn * uabs / ustar * ( th1 - th0 )
271  qstar = c0 * fh * q0qe / tprn * uabs / ustar * ( q1 - q0 )
272 
273  return
274  end subroutine bulkflux_u95
275 
276  !-----------------------------------------------------------------------------
277  !
278  ! refs. Beljaars (1991) and Wilson (2001)
279  !
280  ! If you want to run with the original Beljaars scheme (Beljaars and Holtslag 1994),
281  ! you should fix the stability functions (fm_unstable, fh_unstable, fm_stable, and fh_stable).
282  !
283  !-----------------------------------------------------------------------------
284  subroutine bulkflux_b91w01( &
285  Ustar, & ! (out)
286  tstar, & ! (out)
287  qstar, & ! (out)
288  uabs, & ! (out)
289  t1, & ! (in)
290  t0, & ! (in)
291  p1, & ! (in)
292  p0, & ! (in)
293  q1, & ! (in)
294  q0, & ! (in)
295  u1, & ! (in)
296  v1, & ! (in)
297  z1, & ! (in)
298  pbl, & ! (in)
299  z0m, & ! (in)
300  z0h, & ! (in)
301  z0e ) ! (in)
302  use scale_const, only: &
303  grav => const_grav, &
304  karman => const_karman, &
305  rdry => const_rdry, &
306  cpdry => const_cpdry, &
307  eps => const_eps, &
308  pre00 => const_pre00
309  implicit none
310 
311  ! parameter
312  real(DP), parameter :: dL = 1.0e-6_dp ! delta Obukhov length [m]
313 
314  real(DP), parameter :: Pt = 0.95_dp ! turbulent Prandtl number
315 
316  ! argument
317  real(RP), intent(out) :: Ustar ! friction velocity [m/s]
318  real(RP), intent(out) :: Tstar ! friction temperature [K]
319  real(RP), intent(out) :: Qstar ! friction mixing rate [kg/kg]
320  real(RP), intent(out) :: Uabs ! modified absolute velocity [m/s]
321 
322  real(RP), intent(in) :: T1 ! tempearature at the lowest atmospheric layer [K]
323  real(RP), intent(in) :: T0 ! skin temperature [K]
324  real(RP), intent(in) :: P1 ! pressure at the lowest atmospheric layer [Pa]
325  real(RP), intent(in) :: P0 ! surface pressure [Pa]
326  real(RP), intent(in) :: Q1 ! mixing ratio at the lowest atmospheric layer [kg/kg]
327  real(RP), intent(in) :: Q0 ! surface mixing ratio [kg/kg]
328  real(RP), intent(in) :: U1 ! zonal wind at the lowest atmospheric layer [m/s]
329  real(RP), intent(in) :: V1 ! meridional wind at the lowest atmospheric layer [m/s]
330  real(RP), intent(in) :: Z1 ! height at the lowest atmospheric layer [m]
331  real(RP), intent(in) :: PBL ! the top of atmospheric mixing layer [m]
332  real(RP), intent(in) :: Z0M ! roughness length of momentum [m]
333  real(RP), intent(in) :: Z0H ! roughness length of heat [m]
334  real(RP), intent(in) :: Z0E ! roughness length of moisture [m]
335 
336  ! work
337  integer :: n
338 
339  real(DP) :: L ! Obukhov length [m]
340  real(DP) :: res ! residual
341  real(DP) :: dres ! d(residual)/dL
342 
343  real(DP) :: RiB0 ! bulk Richardson number [no unit]
344  real(DP) :: Wstar, dWstar ! free convection velocity scale [m/s]
345 
346  real(DP) :: UabsUS, UabsS, UabsC
347  real(DP) :: dUabsUS, dUabsS
348 
349  real(DP) :: UstarUS, UstarS, UstarC
350  real(DP) :: TstarUS, TstarS, TstarC
351  real(DP) :: QstarUS, QstarS, QstarC
352 
353  real(DP) :: dUstarUS, dUstarS, dUstarC
354  real(DP) :: dTstarUS, dTstarS, dTstarC
355  real(DP) :: dQstarUS, dQstarS, dQstarC
356 
357  real(DP) :: TH1, TH0
358  real(DP) :: sw, tmp
359 
360  real(DP) :: DP_Z1, DP_Z0M, DP_Z0H, DP_Z0E
361  real(DP) :: log_Z1ovZ0M, log_Z1ovZ0H, log_Z1ovZ0E
362  !---------------------------------------------------------------------------
363 
364  ! convert to DP
365  dp_z1 = real( z1, kind=dp )
366  dp_z0m = real( z0m, kind=dp )
367  dp_z0h = real( z0h, kind=dp )
368  dp_z0e = real( z0e, kind=dp )
369 
370  uabsc = max( sqrt( u1**2 + v1**2 ), bulkflux_uabs_min )
371  th1 = t1 * ( pre00 / p1 )**( rdry / cpdry )
372  th0 = t0 * ( pre00 / p0 )**( rdry / cpdry )
373 
374  ! make log constant
375  log_z1ovz0m = log( dp_z1 / dp_z0m )
376  log_z1ovz0h = log( dp_z1 / dp_z0h )
377  log_z1ovz0e = log( dp_z1 / dp_z0e )
378 
379  ! initial bulk Richardson number
380  rib0 = grav * dp_z1 * ( th1 - th0 ) / ( th1 * uabsc**2 )
381  if( abs( rib0 ) < bulkflux_rib_min ) then
382  rib0 = sign( bulkflux_rib_min, rib0 )
383  end if
384 
385  ! initial Obukhov length assumed by neutral condition
386  l = dp_z1 / rib0 * log_z1ovz0h / log_z1ovz0m**2
387 
388  ! initial free convection velocity scale
389  wstar = bulkflux_wstar_min
390  dwstar = bulkflux_wstar_min
391 
392  do n = 1, bulkflux_itr_max
393  ! unstable condition
394  uabsus = max( sqrt( u1**2 + v1**2 + (bulkflux_wscf*wstar)**2 ), bulkflux_uabs_min )
395  ustarus = karman / ( log_z1ovz0m - fm_unstable(dp_z1,l) + fm_unstable(dp_z0m,l) ) * uabsus
396  tstarus = karman / ( log_z1ovz0h - fh_unstable(dp_z1,l) + fh_unstable(dp_z0h,l) ) / pt * ( th1 - th0 )
397  qstarus = karman / ( log_z1ovz0e - fh_unstable(dp_z1,l) + fh_unstable(dp_z0e,l) ) / pt * ( q1 - q0 )
398 
399  ! stable condition
400  uabss = max( sqrt( u1**2 + v1**2 ), bulkflux_uabs_min )
401  ustars = karman / ( log_z1ovz0m - fm_stable(dp_z1,l) + fm_stable(dp_z0m,l) ) * uabss
402  tstars = karman / ( log_z1ovz0h - fh_stable(dp_z1,l) + fh_stable(dp_z0h,l) ) / pt * ( th1 - th0 )
403  qstars = karman / ( log_z1ovz0e - fh_stable(dp_z1,l) + fh_stable(dp_z0e,l) ) / pt * ( q1 - q0 )
404 
405  sw = 0.5_dp - sign( 0.5_dp, l ) ! if unstable, sw = 1
406 
407  uabsc = ( sw ) * uabsus + ( 1.0_dp-sw ) * uabss
408  ustarc = ( sw ) * ustarus + ( 1.0_dp-sw ) * ustars
409  tstarc = ( sw ) * tstarus + ( 1.0_dp-sw ) * tstars
410  qstarc = ( sw ) * qstarus + ( 1.0_dp-sw ) * qstars
411 
412  ! avoid zero division with TstarC = 0
413  sw = 0.5_dp + sign( 0.5_dp, abs(tstarc) - eps )
414  tstarc = ( sw ) * tstarc + ( 1.0_dp-sw ) * eps
415 
416  ! update free convection velocity scale
417  tmp = -pbl * grav / t1 * ustarc * tstarc
418  sw = 0.5_dp + sign( 0.5_dp, tmp ) ! if tmp is plus, sw = 1
419  wstar = ( tmp * sw )**( 1.0_dp / 3.0_dp )
420 
421  ! calculate residual
422  res = l - ustarc**2 * t1 / ( karman * grav * tstarc )
423 
424  ! unstable condition
425  duabsus = max( sqrt( u1**2 + v1**2 + (bulkflux_wscf*dwstar)**2 ), bulkflux_uabs_min )
426  dustarus = karman / ( log_z1ovz0m - fm_unstable(dp_z1,l+dl) + fm_unstable(dp_z0m,l+dl) ) * duabsus
427  dtstarus = karman / ( log_z1ovz0h - fh_unstable(dp_z1,l+dl) + fh_unstable(dp_z0h,l+dl) ) / pt * ( th1 - th0 )
428  dqstarus = karman / ( log_z1ovz0e - fh_unstable(dp_z1,l+dl) + fh_unstable(dp_z0e,l+dl) ) / pt * ( q1 - q0 )
429  ! stable condition
430  duabss = max( sqrt( u1**2 + v1**2 ), bulkflux_uabs_min )
431  dustars = karman / ( log_z1ovz0m - fm_stable(dp_z1,l+dl) + fm_stable(dp_z0m,l+dl) ) * duabss
432  dtstars = karman / ( log_z1ovz0h - fh_stable(dp_z1,l+dl) + fh_stable(dp_z0h,l+dl) ) / pt * ( th1 - th0 )
433  dqstars = karman / ( log_z1ovz0e - fh_stable(dp_z1,l+dl) + fh_stable(dp_z0e,l+dl) ) / pt * ( q1 - q0 )
434 
435  sw = 0.5_dp - sign( 0.5_dp, l+dl ) ! if unstable, sw = 1
436 
437  dustarc = ( sw ) * dustarus + ( 1.0_dp-sw ) * dustars
438  dtstarc = ( sw ) * dtstarus + ( 1.0_dp-sw ) * dtstars
439  dqstarc = ( sw ) * dqstarus + ( 1.0_dp-sw ) * dqstars
440 
441  ! avoid zero division with dTstarC = 0
442  sw = 0.5_dp + sign( 0.5_dp, abs(dtstarc) - eps )
443  dtstarc = ( sw ) * dtstarc + ( 1.0_dp-sw ) * eps
444 
445  ! update d(free convection velocity scale)
446  tmp = -pbl * grav / t1 * dustarc * dtstarc
447  sw = 0.5_dp + sign( 0.5_dp, tmp ) ! if tmp is plus, sw = 1
448  dwstar = ( tmp * sw )**( 1.0_dp / 3.0_dp )
449 
450  ! calculate d(residual)/dL
451  dres = 1.0_dp - t1 / ( karman * grav * dl ) * ( dustarc**2 / dtstarc - ustarc**2 / tstarc )
452 
453  ! convergence test with residual and error levels
454  if( abs( res ) < bulkflux_res_min .OR. &
455  abs( res/dres ) < bulkflux_err_min ) then
456  exit
457  end if
458 
459  ! stop iteration to prevent numerical error
460  if( abs( dres ) < eps ) then
461  exit
462  end if
463 
464  ! update Obukhov length
465  l = l - res / dres
466 
467  end do
468 
469  ! revert to RP
470  ustar = real( ustarc, kind=rp )
471  tstar = real( tstarc, kind=rp )
472  qstar = real( qstarc, kind=rp )
473  uabs = real( uabsc, kind=rp )
474 
475  if( n > bulkflux_itr_max ) then
476  if( io_l ) write(io_fid_log,'(A)' ) 'Warning: reach maximum iteration in the function of BULKFLUX_B91W01.'
477  if( io_l ) write(io_fid_log,'(A)' ) ''
478  if( io_l ) write(io_fid_log,'(A,I32)' ) 'DEBUG --- LOOP number [no unit] :', n
479  if( io_l ) write(io_fid_log,'(A,F32.16)') 'DEBUG --- Residual [m] :', res
480  if( io_l ) write(io_fid_log,'(A,F32.16)') 'DEBUG --- delta Residual [m] :', dres
481  if( io_l ) write(io_fid_log,'(A)' ) ''
482  if( io_l ) write(io_fid_log,'(A,F32.16)') 'DEBUG --- air tempearature [K] :', t1
483  if( io_l ) write(io_fid_log,'(A,F32.16)') 'DEBUG --- surface temperature [K] :', t0
484  if( io_l ) write(io_fid_log,'(A,F32.16)') 'DEBUG --- pressure [Pa] :', p1
485  if( io_l ) write(io_fid_log,'(A,F32.16)') 'DEBUG --- surface pressure [Pa] :', p0
486  if( io_l ) write(io_fid_log,'(A,F32.16)') 'DEBUG --- water vapor mass ratio [kg/kg] :', q1
487  if( io_l ) write(io_fid_log,'(A,F32.16)') 'DEBUG --- surface water vapor mass ratio [kg/kg] :', q0
488  if( io_l ) write(io_fid_log,'(A,F32.16)') 'DEBUG --- zonal wind [m/s] :', u1
489  if( io_l ) write(io_fid_log,'(A,F32.16)') 'DEBUG --- meridional wind [m/s] :', v1
490  if( io_l ) write(io_fid_log,'(A,F32.16)') 'DEBUG --- cell center height [m] :', z1
491  if( io_l ) write(io_fid_log,'(A,F32.16)') 'DEBUG --- atmospheric mixing layer height [m] :', pbl
492  if( io_l ) write(io_fid_log,'(A,F32.16)') 'DEBUG --- roughness length of momentum [m] :', z0m
493  if( io_l ) write(io_fid_log,'(A,F32.16)') 'DEBUG --- roughness length of heat [m] :', z0h
494  if( io_l ) write(io_fid_log,'(A,F32.16)') 'DEBUG --- roughness length of moisture [m] :', z0e
495  if( io_l ) write(io_fid_log,'(A)' ) ''
496  if( io_l ) write(io_fid_log,'(A,F32.16)') 'DEBUG --- friction velocity [m] :', ustar
497  if( io_l ) write(io_fid_log,'(A,F32.16)') 'DEBUG --- friction potential temperature [K] :', tstar
498  if( io_l ) write(io_fid_log,'(A,F32.16)') 'DEBUG --- friction water vapor mass ratio [kg/kg] :', qstar
499  if( io_l ) write(io_fid_log,'(A)' ) ''
500  if( io_l ) write(io_fid_log,'(A,F32.16)') 'DEBUG --- Obukhov length [m] :', l
501  if( io_l ) write(io_fid_log,'(A,F32.16)') 'DEBUG --- bulk Richardson number [no unit] :', rib0
502  if( io_l ) write(io_fid_log,'(A,F32.16)') 'DEBUG --- free convection velocity scale [m/s] :', wstar
503  end if
504 
505  return
506  end subroutine bulkflux_b91w01
507 
508  !-----------------------------------------------------------------------------
509  ! stability function for momemtum in unstable condition
510  function fm_unstable( Z, L )
511 ! use scale_const, only: &
512 ! PI => CONST_PI
513  implicit none
514 
515  ! argument
516  real(DP), intent(in) :: z
517  real(DP), intent(in) :: L
518 
519  ! function
520  real(DP) :: fm_unstable
521 
522  ! works
523  real(DP) :: R
524 ! real(DP) :: r4R
525  !---------------------------------------------------------------------------
526 
527  r = min( z/l, 0.0_dp )
528 
529  ! Wilson (2001)
530  fm_unstable = 3.0_dp * log( ( 1.0_dp + sqrt( 1.0_dp + 3.6_dp * (-r)**(2.0_dp/3.0_dp) ) ) * 0.5_dp )
531 
532  ! If you want to run with the original Beljaars scheme (Beljaars and Holtslag 1994),
533  ! you should comment out the above line (Wilson 2001) and uncomment the below lines (Paulson 1974; Dyer 1974).
534  !
535  !! Paulson (1974); Dyer (1974)
536  !r4R = ( 1.0_DP - 16.0_DP * R )**0.25_DP
537  !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
538 
539  return
540  end function fm_unstable
541 
542  !-----------------------------------------------------------------------------
543  ! stability function for heat/vapor in unstable condition
544  function fh_unstable( Z, L )
545  implicit none
546 
547  ! argument
548  real(DP), intent(in) :: Z
549  real(DP), intent(in) :: L
550 
551  ! function
552  real(DP) :: fh_unstable
553 
554  ! works
555  real(DP) :: R
556  !---------------------------------------------------------------------------
557 
558  r = min( z/l, 0.0_dp )
559 
560  ! Wilson (2001)
561  fh_unstable = 3.0_dp * log( ( 1.0_dp + sqrt( 1.0_dp + 7.9_dp * (-r)**(2.0_dp/3.0_dp) ) ) * 0.5_dp )
562 
563  ! If you want to run with the original Beljaars scheme (Beljaars and Holtslag 1994),
564  ! you should comment out the above line (Wilson 2001) and uncomment the below lines (Paulson 1974; Dyer 1974).
565  !
566  !! Paulson (1974); Dyer (1974)
567  !fh_unstable = 2.0_DP * log( ( 1.0_DP + sqrt( 1.0_DP - 16.0_DP * R ) ) * 0.5_DP )
568 
569  return
570  end function fh_unstable
571 
572  !-----------------------------------------------------------------------------
573  ! stability function for momemtum in stable condition
574  function fm_stable( Z, L )
575  implicit none
576 
577  ! argument
578  real(DP), intent(in) :: Z
579  real(DP), intent(in) :: L
580 
581  ! function
582  real(DP) :: fm_stable
583 
584  ! parameters of stability functions (Beljaars and Holtslag 1991)
585  real(DP), parameter :: a = 1.0_dp
586  real(DP), parameter :: b = 0.667_dp
587  real(DP), parameter :: c = 5.0_dp
588  real(DP), parameter :: d = 0.35_dp
589 
590  ! works
591  real(DP) :: R
592  !---------------------------------------------------------------------------
593 
594  r = max( z/l, 0.0_dp )
595 
596  ! Holtslag and DeBruin (1988)
597  fm_stable = - a*r - b*( r - c/d )*exp( -d*r ) - b*c/d
598 
599  return
600  end function fm_stable
601 
602  !-----------------------------------------------------------------------------
603  ! stability function for heat/vapor in stable condition
604  function fh_stable( Z, L )
605  implicit none
606 
607  ! argument
608  real(DP), intent(in) :: Z
609  real(DP), intent(in) :: L
610 
611  ! function
612  real(DP) :: fh_stable
613 
614  ! parameters of stability functions (Beljaars and Holtslag 1991)
615  real(DP), parameter :: a = 1.0_dp
616  real(DP), parameter :: b = 0.667_dp
617  real(DP), parameter :: c = 5.0_dp
618  real(DP), parameter :: d = 0.35_dp
619 
620  ! works
621  real(DP) :: R
622  !---------------------------------------------------------------------------
623 
624  r = max( z/l, 0.0_dp )
625 
626  ! Beljaars and Holtslag (1991)
627  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
628 
629  return
630  end function fh_stable
631 
632 end module scale_bulkflux
scale_const::const_cpdry
real(rp), public const_cpdry
specific heat (dry air,constant pressure) [J/kg/K]
Definition: scale_const.F90:58
scale_process::prc_mpistop
subroutine, public prc_mpistop
Abort MPI.
Definition: scale_process.F90:234
scale_stdio::io_l
logical, public io_l
output log or not? (this process)
Definition: scale_stdio.F90:59
scale_stdio
module STDIO
Definition: scale_stdio.F90:12
scale_const::const_karman
real(rp), parameter, public const_karman
von Karman constant
Definition: scale_const.F90:52
scale_const::const_rdry
real(rp), public const_rdry
specific gas constant (dry air) [J/kg/K]
Definition: scale_const.F90:57
scale_bulkflux::bulkflux_setup
subroutine, public bulkflux_setup
Definition: scale_bulkflux.F90:113
dc_types::dp
integer, parameter, public dp
Definition: dc_types.f90:27
scale_const::const_pre00
real(rp), public const_pre00
pressure reference [Pa]
Definition: scale_const.F90:93
scale_bulkflux::bulkflux
procedure(bc), pointer, public bulkflux
Definition: scale_bulkflux.F90:71
scale_const::const_grav
real(rp), public const_grav
standard acceleration of gravity [m/s2]
Definition: scale_const.F90:48
scale_process
module PROCESS
Definition: scale_process.F90:14
dc_log::log
subroutine, public log(type, message)
Definition: dc_log.f90:133
scale_const
module CONSTANT
Definition: scale_const.F90:14
scale_bulkflux
module Surface bulk flux
Definition: scale_bulkflux.F90:12
scale_const::const_eps
real(rp), public const_eps
small number
Definition: scale_const.F90:36
scale_stdio::io_lnml
logical, public io_lnml
output log or not? (for namelist, this process)
Definition: scale_stdio.F90:60
scale_precision
module PRECISION
Definition: scale_precision.F90:13
scale_stdio::io_fid_conf
integer, public io_fid_conf
Config file ID.
Definition: scale_stdio.F90:55
scale_stdio::io_fid_log
integer, public io_fid_log
Log file ID.
Definition: scale_stdio.F90:56
scale_precision::rp
integer, parameter, public rp
Definition: scale_precision.F90:43
Generated by   doxygen 1.8.13