SCALE-RM
scale_atmos_phy_sf_bulkcoef.F90
Go to the documentation of this file.
1 !-------------------------------------------------------------------------------
13 !-------------------------------------------------------------------------------
14 #include "inc_openmp.h"
15 module scale_atmos_phy_sf_bulkcoef
16  !-----------------------------------------------------------------------------
17  !
18  !++ used modules
19  !
20  use scale_precision
21  use scale_stdio
22  use scale_prof
23  use scale_grid_index
24  use scale_tracer
25  !-----------------------------------------------------------------------------
26  implicit none
27  private
28  !-----------------------------------------------------------------------------
29  !
30  !++ Public procedure
31  !
32  public :: atmos_phy_sf_bulkcoef_setup
33 
34  abstract interface
35  subroutine bc( &
36  ATM_Uabs, ATM_POTT, Z1, &
37  SFC_TEMP, &
38  Z0M, Z0H, Z0E, &
39  Cm, Ch, Ce , &
40  R10M, R2H, R2E )
41  use scale_precision
42  use scale_grid_index
43  use scale_tracer
44  implicit none
45 
46  real(RP), intent(in) :: atm_uabs(ia,ja) ! absolute velocity at Z1 [m/s]
47  real(RP), intent(in) :: atm_pott(ia,ja) ! potential temperature at z1, based on the local surface pressure [K]
48  real(RP), intent(in) :: z1 (ia,ja) ! height of lowermost atmosphere grid (cell center) [m]
49  real(RP), intent(in) :: sfc_temp(ia,ja) ! potential temperature at surface skin [K]
50  real(RP), intent(in) :: z0m (ia,ja) ! roughness length for momentum [m]
51  real(RP), intent(in) :: z0h (ia,ja) ! roughness length for heat [m]
52  real(RP), intent(in) :: z0e (ia,ja) ! roughness length for moisture [m]
53  real(RP), intent(out) :: cm (ia,ja) ! bulk coefficient for momentum
54  real(RP), intent(out) :: ch (ia,ja) ! bulk coefficient for heat
55  real(RP), intent(out) :: ce (ia,ja) ! bulk coefficient for moisture
56  real(RP), intent(out) :: r10m (ia,ja)
57  real(RP), intent(out) :: r2h (ia,ja)
58  real(RP), intent(out) :: r2e (ia,ja)
59  end subroutine bc
60  end interface
61 
62  procedure(bc), pointer :: atmos_phy_sf_bulkcoef => null()
63 
64  public :: atmos_phy_sf_bulkcoef
65 
66  !-----------------------------------------------------------------------------
67  !
68  !++ Public parameters & variables
69  !
70  !-----------------------------------------------------------------------------
71  !
72  !++ Private procedure
73  !
74  private :: atmos_phy_sf_bulkcoef_uno
75  private :: atmos_phy_sf_bulkcoef_beljaars
76  private :: fm_unstable
77  private :: fh_unstable
78  private :: fm_stable
79  private :: fh_stable
80 
81  !-----------------------------------------------------------------------------
82  !
83  !++ Private parameters & variables
84  !
85  character(len=H_SHORT), private :: atmos_phy_sf_bulkcoef_type = 'BH91'
86 
87  real(RP), private :: atmos_phy_sf_bulkcoef_cm_max = 2.5e-3_rp ! maximum limit of bulk coefficient for momentum [NIL]
88  real(RP), private :: atmos_phy_sf_bulkcoef_ch_max = 1.0e+0_rp ! maximum limit of bulk coefficient for heat [NIL]
89  real(RP), private :: atmos_phy_sf_bulkcoef_ce_max = 1.0e+0_rp ! maximum limit of bulk coefficient for moisture [NIL]
90  real(RP), private :: atmos_phy_sf_bulkcoef_cm_min = 1.0e-5_rp ! minimum limit of bulk coefficient for momentum [NIL]
91  real(RP), private :: atmos_phy_sf_bulkcoef_ch_min = 1.0e-5_rp ! minimum limit of bulk coefficient for heat [NIL]
92  real(RP), private :: atmos_phy_sf_bulkcoef_ce_min = 1.0e-5_rp ! minimum limit of bulk coefficient for moisture [NIL]
93 
94  !-----------------------------------------------------------------------------
95 contains
96  !-----------------------------------------------------------------------------
97  subroutine atmos_phy_sf_bulkcoef_setup
98  use scale_process, only: &
99  prc_mpistop
100  implicit none
101 
102  namelist / param_atmos_phy_sf_bulkcoef / &
103  atmos_phy_sf_bulkcoef_type, &
104  atmos_phy_sf_bulkcoef_cm_max, &
105  atmos_phy_sf_bulkcoef_ch_max, &
106  atmos_phy_sf_bulkcoef_ce_max, &
107  atmos_phy_sf_bulkcoef_cm_min, &
108  atmos_phy_sf_bulkcoef_ch_min, &
109  atmos_phy_sf_bulkcoef_ce_min
110 
111  integer :: ierr
112  !---------------------------------------------------------------------------
113 
114  if( io_l ) write(io_fid_log,*)
115  if( io_l ) write(io_fid_log,*) '*** Bulk coefficient parameter'
116 
117  !--- read namelist
118  rewind(io_fid_conf)
119  read(io_fid_conf,nml=param_atmos_phy_sf_bulkcoef,iostat=ierr)
120  if( ierr < 0 ) then !--- missing
121  if( io_l ) write(io_fid_log,*) '*** Not found namelist. Default used.'
122  elseif( ierr > 0 ) then !--- fatal error
123  write(*,*) 'xxx Not appropriate names in namelist PARAM_ATMOS_PHY_SF_BULKCOEF. Check!'
124  call prc_mpistop
125  endif
126  if( io_lnml ) write(io_fid_log,nml=param_atmos_phy_sf_bulkcoef)
127 
128  select case( atmos_phy_sf_bulkcoef_type )
129  case ('U95')
130  atmos_phy_sf_bulkcoef => atmos_phy_sf_bulkcoef_uno
131  case ('BH91')
132  atmos_phy_sf_bulkcoef => atmos_phy_sf_bulkcoef_beljaars
133  case default
134  write(*,*) 'xxx invalid bulk scheme (', trim(atmos_phy_sf_bulkcoef_type), '). CHECK!'
135  call prc_mpistop
136  end select
137 
138  return
139  end subroutine atmos_phy_sf_bulkcoef_setup
140 
141  !-----------------------------------------------------------------------------
142  subroutine atmos_phy_sf_bulkcoef_uno( &
143  ATM_Uabs, ATM_POTT, Z1, &
144  SFC_TEMP, &
145  Z0M, Z0H, Z0E, &
146  Cm, Ch, Ce , &
147  R10M, R2H, R2E )
148  use scale_const, only: &
149  eps => const_eps, &
150  grav => const_grav, &
151  karman => const_karman
152  implicit none
153 
154  real(RP), intent(in) :: ATM_Uabs(ia,ja) ! absolute velocity at Z1 [m/s]
155  real(RP), intent(in) :: ATM_POTT(ia,ja) ! potential temperature at z1, based on the local surface pressure [K]
156  real(RP), intent(in) :: Z1 (ia,ja) ! height of lowermost atmosphere grid (cell center) [m]
157  real(RP), intent(in) :: SFC_TEMP(ia,ja) ! potential temperature at surface skin [K]
158  real(RP), intent(in) :: Z0M (ia,ja) ! roughness length for momentum [m]
159  real(RP), intent(in) :: Z0H (ia,ja) ! roughness length for heat [m]
160  real(RP), intent(in) :: Z0E (ia,ja) ! roughness length for moisture [m]
161  real(RP), intent(out) :: Cm (ia,ja) ! bulk coefficient for momentum
162  real(RP), intent(out) :: Ch (ia,ja) ! bulk coefficient for heat
163  real(RP), intent(out) :: Ce (ia,ja) ! bulk coefficient for moisture
164  real(RP), intent(out) :: R10M (ia,ja)
165  real(RP), intent(out) :: R2H (ia,ja)
166  real(RP), intent(out) :: R2E (ia,ja)
167 
168  ! specific parameters
169  real(RP), parameter :: tPrn = 0.74_rp ! turbulent Prandtl number (Businger et al. 1971)
170  real(RP), parameter :: LFb = 9.4_rp ! Louis factor b (Louis 1979)
171  real(RP), parameter :: LFbp = 4.7_rp ! Louis factor b' (Louis 1979)
172  real(RP), parameter :: LFdm = 7.4_rp ! Louis factor d for momemtum (Louis 1979)
173  real(RP), parameter :: LFdh = 5.3_rp ! Louis factor d for heat (Louis 1979)
174 
175  real(RP) :: RiB (ia,ja) ! bulk Richardson number [no unit]
176  real(RP) :: RiBT(ia,ja)
177  real(RP) :: C0 (ia,ja) ! initial drag coefficient [no unit]
178  real(RP) :: fm (ia,ja)
179  real(RP) :: fh (ia,ja)
180  real(RP) :: Psi (ia,ja)
181  real(RP) :: t0th(ia,ja)
182  real(RP) :: q0qe(ia,ja)
183 
184  real(RP) :: C0_10m(ia,ja)
185  real(RP) :: fm_10m(ia,ja)
186  real(RP) :: C0_2m (ia,ja)
187  real(RP) :: fm_2m (ia,ja)
188  real(RP) :: fh_2m (ia,ja)
189  real(RP) :: Psi_2m(ia,ja)
190 
191  real(RP) :: tmp1, tmp2, tmp3, Psi1
192  integer :: i, j
193  !---------------------------------------------------------------------------
194 
195  do j = js, je
196  do i = is, ie
197  c0(i,j) = ( karman / log( z1(i,j)/z0m(i,j) ) )**2
198  c0_10m(i,j) = ( karman / log( 10.0_rp/z0m(i,j) ) )**2
199  c0_2m(i,j) = ( karman / log( 2.0_rp/z0m(i,j) ) )**2
200 
201  ribt(i,j) = grav * z1(i,j) * ( atm_pott(i,j) - sfc_temp(i,j) ) &
202  / ( atm_pott(i,j) * atm_uabs(i,j)**2 + eps )
203 
204  rib(i,j) = ribt(i,j)
205  enddo
206  enddo
207 
208  do j = js, je
209  do i = is, ie
210  if ( ribt(i,j) < 0.0_rp ) then ! unstable condition
211  tmp1 = c0(i,j) * lfb * sqrt( z1(i,j)/z0m(i,j) ) * sqrt( abs(rib(i,j)) )
212 
213  fm(i,j) = 1.0_rp - lfb * rib(i,j) / ( 1.0_rp + lfdm * tmp1 )
214  fh(i,j) = 1.0_rp - lfb * rib(i,j) / ( 1.0_rp + lfdh * tmp1 )
215  else ! stable condition
216  fm(i,j) = 1.0_rp / ( 1.0_rp + lfbp * rib(i,j) )**2
217  fh(i,j) = fm(i,j)
218  endif
219  enddo
220  enddo
221 
222  do j = js, je
223  do i = is, ie
224  psi(i,j) = tprn * sqrt(fm(i,j)) / fh(i,j) * log( z1(i,j)/z0m(i,j) )
225 
226  t0th(i,j) = 1.0_rp / ( 1.0_rp + tprn * log( z0m(i,j)/z0h(i,j) ) / psi(i,j) )
227  q0qe(i,j) = 1.0_rp / ( 1.0_rp + tprn * log( z0m(i,j)/z0e(i,j) ) / psi(i,j) )
228 
229  rib(i,j) = rib(i,j) * t0th(i,j)
230  enddo
231  enddo
232 
233  do j = js, je
234  do i = is, ie
235  if ( ribt(i,j) < 0.0_rp ) then ! unstable condition
236  tmp1 = c0(i,j) * lfb * sqrt( z1(i,j)/z0m(i,j) ) * sqrt( abs(rib(i,j)) )
237  tmp2 = c0_10m(i,j) * lfb * sqrt( 10.0_rp/z0m(i,j) ) * sqrt( abs(rib(i,j)) )
238  tmp3 = c0_2m(i,j) * lfb * sqrt( 2.0_rp/z0m(i,j) ) * sqrt( abs(rib(i,j)) )
239 
240  fm(i,j) = 1.0_rp - lfb * rib(i,j) / ( 1.0_rp + lfdm * tmp1 )
241  fh(i,j) = 1.0_rp - lfb * rib(i,j) / ( 1.0_rp + lfdh * tmp1 )
242 
243  fm_10m(i,j) = 1.0_rp - lfb * rib(i,j) / ( 1.0_rp + lfdm * tmp2 )
244  fm_2m(i,j) = 1.0_rp - lfb * rib(i,j) / ( 1.0_rp + lfdm * tmp3 )
245  fh_2m(i,j) = 1.0_rp - lfb * rib(i,j) / ( 1.0_rp + lfdh * tmp3 )
246  else ! stable condition
247  fm(i,j) = 1.0_rp / ( 1.0_rp + lfbp * rib(i,j) )**2
248  fh(i,j) = fm(i,j)
249 
250  fm_10m(i,j) = fm(i,j)
251  fm_2m(i,j) = fm(i,j)
252  fh_2m(i,j) = fm(i,j)
253  endif
254  enddo
255  enddo
256 
257  do j = js, je
258  do i = is, ie
259  psi1 = tprn * sqrt(fm(i,j)) / fh(i,j) * log( z1(i,j)/z0m(i,j) )
260  psi_2m(i,j) = tprn * sqrt(fm_2m(i,j)) / fh_2m(i,j) * log( 2.0_rp/z0m(i,j) )
261 
262  t0th(i,j) = 1.0_rp / ( 1.0_rp + tprn * log( z0m(i,j)/z0h(i,j) ) / psi1 )
263  q0qe(i,j) = 1.0_rp / ( 1.0_rp + tprn * log( z0m(i,j)/z0e(i,j) ) / psi1 )
264  enddo
265  enddo
266 
267  do j = js, je
268  do i = is, ie
269  cm(i,j) = c0(i,j) * fm(i,j)
270  ch(i,j) = c0(i,j) * fh(i,j) * t0th(i,j) / tprn
271  ce(i,j) = c0(i,j) * fh(i,j) * q0qe(i,j) / tprn
272 
273  cm(i,j) = min( max( cm(i,j), atmos_phy_sf_bulkcoef_cm_min ), atmos_phy_sf_bulkcoef_cm_max )
274  ch(i,j) = min( max( ch(i,j), atmos_phy_sf_bulkcoef_ch_min ), atmos_phy_sf_bulkcoef_ch_max )
275  ce(i,j) = min( max( ce(i,j), atmos_phy_sf_bulkcoef_ce_min ), atmos_phy_sf_bulkcoef_ce_max )
276  enddo
277  enddo
278 
279  do j = js, je
280  do i = is, ie
281  r10m(i,j) = log( 10.0_rp/z0m(i,j) ) / log( z1(i,j)/z0m(i,j) ) &
282  / sqrt( fm_10m(i,j) / fm(i,j) )
283 
284  tmp1 = log( z0m(i,j)/z0h(i,j) ) / log( z1(i,j)/z0m(i,j) )
285 
286  r2h(i,j) = ( tmp1 + psi_2m(i,j) / tprn ) &
287  / ( tmp1 + psi(i,j) / tprn )
288 
289  tmp2 = log( z0m(i,j)/z0e(i,j) ) / log( z1(i,j)/z0m(i,j) )
290 
291  r2e(i,j) = ( tmp2 + psi_2m(i,j) / tprn ) &
292  / ( tmp2 + psi(i,j) / tprn )
293 
294  r10m(i,j) = min( max( r10m(i,j), 0.0_rp ), 1.0_rp )
295  r2h(i,j) = min( max( r2h(i,j), 0.0_rp ), 1.0_rp )
296  r2e(i,j) = min( max( r2e(i,j), 0.0_rp ), 1.0_rp )
297  enddo
298  enddo
299 
300  return
301  end subroutine atmos_phy_sf_bulkcoef_uno
302 
303  !-----------------------------------------------------------------------------
304  subroutine atmos_phy_sf_bulkcoef_beljaars( &
305  ATM_Uabs, ATM_POTT, Z1, &
306  SFC_TEMP, &
307  Z0M, Z0H, Z0E, &
308  Cm, Ch, Ce , &
309  R10M, R2H, R2E )
310  use scale_const, only: &
311  eps => const_eps, &
312  grav => const_grav, &
313  karman => const_karman
314  implicit none
315 
316  real(RP), intent(in) :: ATM_Uabs(ia,ja) ! absolute velocity at Z1 [m/s]
317  real(RP), intent(in) :: ATM_POTT(ia,ja) ! potential temperature at z1, based on the local surface pressure [K]
318  real(RP), intent(in) :: Z1 (ia,ja) ! height of lowermost atmosphere grid (cell center) [m]
319  real(RP), intent(in) :: SFC_TEMP(ia,ja) ! potential temperature at surface skin [K]
320  real(RP), intent(in) :: Z0M (ia,ja) ! roughness length for momentum [m]
321  real(RP), intent(in) :: Z0H (ia,ja) ! roughness length for heat [m]
322  real(RP), intent(in) :: Z0E (ia,ja) ! roughness length for moisture [m]
323  real(RP), intent(out) :: Cm (ia,ja) ! bulk coefficient for momentum
324  real(RP), intent(out) :: Ch (ia,ja) ! bulk coefficient for heat
325  real(RP), intent(out) :: Ce (ia,ja) ! bulk coefficient for moisture
326  real(RP), intent(out) :: R10M (ia,ja)
327  real(RP), intent(out) :: R2H (ia,ja)
328  real(RP), intent(out) :: R2E (ia,ja)
329 
330  ! specific parameters
331  real(RP), parameter :: RiB_min = 1.e-4_rp
332  real(RP), parameter :: res_criteria = 1.e-6_rp
333  real(RP), parameter :: dL = 1.e-10_rp ! delta Obukhov length [m]
334 
335  real(RP) :: RiB0(ia,ja) ! bulk Richardson number [no unit]
336  real(RP) :: L (ia,ja) ! Obukhov length [m]
337  real(RP) :: Lpls(ia,ja) ! Obukhov length [m]
338 
339  real(RP) :: LogZ1Z0M(ia,ja)
340  real(RP) :: LogZ1Z0H(ia,ja)
341  real(RP) :: LogZ1Z0E(ia,ja)
342 
343  real(RP) :: CmUS(ia,ja), CmS(ia,ja)
344  real(RP) :: ChUS(ia,ja), ChS(ia,ja)
345  real(RP) :: CeUS(ia,ja), CeS(ia,ja)
346  real(RP) :: dCmUS(ia,ja), dCmS(ia,ja), dCm(ia,ja)
347  real(RP) :: dChUS(ia,ja), dChS(ia,ja), dCh(ia,ja)
348  real(RP) :: dCeUS(ia,ja), dCeS(ia,ja), dCe(ia,ja)
349  real(RP) :: res (ia,ja), dres
350 
351  real(RP) :: fmUS_Z1L, fmUS_Z0ML, fhUS_Z1L, fhUS_Z0HL, fhUS_Z0EL
352  real(RP) :: fmS_Z1L, fmS_Z0ML, fhS_Z1L, fhS_Z0HL, fhS_Z0EL
353 
354  real(RP) :: fmUS_10mL, fmUS_2mL, fhUS_2mL
355  real(RP) :: fmS_10mL, fmS_2mL, fhS_2mL
356  real(RP) :: fm (ia,ja)
357  real(RP) :: fh (ia,ja)
358  real(RP) :: fm_10m(ia,ja)
359  real(RP) :: fm_2m (ia,ja)
360  real(RP) :: fh_2m (ia,ja)
361  real(RP) :: Psi (ia,ja)
362  real(RP) :: Psi_2m(ia,ja)
363 
364  integer, parameter :: itelim = 100 ! maximum iteration number
365 
366  real(RP) :: KARMAN2
367  real(RP) :: tmp1, tmp2, sw
368  integer :: i, j, ite
369  !---------------------------------------------------------------------------
370 
371  karman2 = karman * karman
372 
373  do j = js, je
374  do i = is, ie
375  rib0(i,j) = grav * z1(i,j) * ( atm_pott(i,j) - sfc_temp(i,j) ) &
376  / ( atm_pott(i,j) * atm_uabs(i,j)**2 + eps )
377  rib0(i,j) = max( abs(rib0(i,j)), rib_min )
378  enddo
379  enddo
380 
381  ! The initial Obukhov length is assumed by bulk Richardson number.
382  do j = js, je
383  do i = is, ie
384  l(i,j) = z1(i,j) / rib0(i,j)
385  enddo
386  enddo
387 
388  do j = js, je
389  do i = is, ie
390  logz1z0m(i,j) = log( z1(i,j) / z0m(i,j) )
391  logz1z0h(i,j) = log( z1(i,j) / z0h(i,j) )
392  logz1z0e(i,j) = log( z1(i,j) / z0e(i,j) )
393  enddo
394  enddo
395 
396  do ite = 1, itelim
397 
398  !---< calc coefficient >---
399 
400  ! unstable condition
401  do j = js, je
402  do i = is, ie
403  fmus_z1l = fm_unstable( z1(i,j),l(i,j) )
404  fmus_z0ml = fm_unstable( z0m(i,j),l(i,j) )
405 
406  fhus_z1l = fh_unstable( z1(i,j),l(i,j) )
407  fhus_z0hl = fh_unstable( z0h(i,j),l(i,j) )
408  fhus_z0el = fh_unstable( z0e(i,j),l(i,j) )
409 
410  cmus(i,j) = karman2 / ( logz1z0m(i,j) - fmus_z1l + fmus_z0ml )**2
411  chus(i,j) = karman2 / ( logz1z0m(i,j) - fmus_z1l + fmus_z0ml ) / ( logz1z0h(i,j) - fhus_z1l + fhus_z0hl )
412  ceus(i,j) = karman2 / ( logz1z0m(i,j) - fmus_z1l + fmus_z0ml ) / ( logz1z0e(i,j) - fhus_z1l + fhus_z0el )
413  enddo
414  enddo
415 
416  ! stable condition
417  do j = js, je
418  do i = is, ie
419  fms_z1l = fm_stable( z1(i,j),l(i,j) )
420  fms_z0ml = fm_stable( z0m(i,j),l(i,j) )
421 
422  fhs_z1l = fh_stable( z1(i,j),l(i,j) )
423  fhs_z0hl = fh_stable( z0h(i,j),l(i,j) )
424  fhs_z0el = fh_stable( z0e(i,j),l(i,j) )
425 
426  cms(i,j) = karman2 / ( logz1z0m(i,j) - fms_z1l + fms_z0ml )**2
427  chs(i,j) = karman2 / ( logz1z0m(i,j) - fms_z1l + fms_z0ml ) / ( logz1z0h(i,j) - fhs_z1l + fhs_z0hl )
428  ces(i,j) = karman2 / ( logz1z0m(i,j) - fms_z1l + fms_z0ml ) / ( logz1z0e(i,j) - fhs_z1l + fhs_z0el )
429  enddo
430  enddo
431 
432  ! select unstable / stable
433  do j = js, je
434  do i = is, ie
435  sw = 0.5_rp - sign(0.5_rp,l(i,j)) ! if unstable, sw = 1
436 
437  cm(i,j) = ( sw ) * cmus(i,j) &
438  + ( 1.0_rp-sw ) * cms(i,j)
439  ch(i,j) = ( sw ) * chus(i,j) &
440  + ( 1.0_rp-sw ) * chs(i,j)
441  ce(i,j) = ( sw ) * ceus(i,j) &
442  + ( 1.0_rp-sw ) * ces(i,j)
443  enddo
444  enddo
445 
446  !---< calc derivative >---
447 
448  do j = js, je
449  do i = is, ie
450  lpls(i,j) = l(i,j) + dl
451  enddo
452  enddo
453 
454  ! unstable condition
455  do j = js, je
456  do i = is, ie
457  fmus_z1l = fm_unstable( z1(i,j),lpls(i,j) )
458  fmus_z0ml = fm_unstable( z0m(i,j),lpls(i,j) )
459 
460  fhus_z1l = fh_unstable( z1(i,j),lpls(i,j) )
461  fhus_z0hl = fh_unstable( z0h(i,j),lpls(i,j) )
462  fhus_z0el = fh_unstable( z0e(i,j),lpls(i,j) )
463 
464  dcmus(i,j) = karman2 / ( logz1z0m(i,j) - fmus_z1l + fmus_z0ml )**2
465  dchus(i,j) = karman2 / ( logz1z0m(i,j) - fmus_z1l + fmus_z0ml ) / ( logz1z0h(i,j) - fhus_z1l + fhus_z0hl )
466  dceus(i,j) = karman2 / ( logz1z0m(i,j) - fmus_z1l + fmus_z0ml ) / ( logz1z0e(i,j) - fhus_z1l + fhus_z0el )
467  enddo
468  enddo
469 
470  ! stable condition
471  do j = js, je
472  do i = is, ie
473  fms_z1l = fm_stable( z1(i,j),lpls(i,j) )
474  fms_z0ml = fm_stable( z0m(i,j),lpls(i,j) )
475 
476  fhs_z1l = fh_stable( z1(i,j),lpls(i,j) )
477  fhs_z0hl = fh_stable( z0h(i,j),lpls(i,j) )
478  fhs_z0el = fh_stable( z0e(i,j),lpls(i,j) )
479 
480  dcms(i,j) = karman2 / ( logz1z0m(i,j) - fms_z1l + fms_z0ml )**2
481  dchs(i,j) = karman2 / ( logz1z0m(i,j) - fms_z1l + fms_z0ml ) / ( logz1z0h(i,j) - fhs_z1l + fhs_z0hl )
482  dces(i,j) = karman2 / ( logz1z0m(i,j) - fms_z1l + fms_z0ml ) / ( logz1z0e(i,j) - fhs_z1l + fhs_z0el )
483  enddo
484  enddo
485 
486  do j = js, je
487  do i = is, ie
488  sw = 0.5_rp - sign(0.5_rp,lpls(i,j)) ! if unstable, sw = 1
489 
490  dcm(i,j) = ( sw ) * dcmus(i,j) &
491  + ( 1.0_rp-sw ) * dcms(i,j)
492  dch(i,j) = ( sw ) * dchus(i,j) &
493  + ( 1.0_rp-sw ) * dchs(i,j)
494  dce(i,j) = ( sw ) * dceus(i,j) &
495  + ( 1.0_rp-sw ) * dces(i,j)
496  enddo
497  enddo
498 
499  ! update Obukhov length
500  do j = js, je
501  do i = is, ie
502  ! residual
503  res(i,j) = l(i,j) - z1(i,j) * cm(i,j)**1.5_rp / ( karman * ch(i,j) * rib0(i,j) )
504  ! d(residual)
505  dres = lpls(i,j) - z1(i,j) * dcm(i,j)**1.5_rp / ( karman * dch(i,j) * rib0(i,j) )
506 
507  l(i,j) = l(i,j) - res(i,j) / ( dres - res(i,j) ) * dl
508  enddo
509  enddo
510 
511  if( maxval( abs( res(is:ie,js:je) ) ) < res_criteria ) exit
512  enddo
513 
514  ! result
515 ! RiB(:,:) = Z1(:,:) * Cm(:,:)**1.5_RP / ( L(:,:) * KARMAN * Ch(:,:) )
516 
517  do j = js, je
518  do i = is, ie
519  ! unstable condition
520  fmus_z1l = fm_unstable( z1(i,j),l(i,j) )
521  fhus_z1l = fh_unstable( z1(i,j),l(i,j) )
522  fmus_10ml = fm_unstable( 10.0_rp,l(i,j) )
523  fmus_2ml = fm_unstable( 2.0_rp,l(i,j) )
524  fhus_2ml = fh_unstable( 2.0_rp,l(i,j) )
525  ! stable condition
526  fms_z1l = fm_stable( z1(i,j),l(i,j) )
527  fhs_z1l = fh_stable( z1(i,j),l(i,j) )
528  fms_10ml = fm_stable( 10.0_rp,l(i,j) )
529  fms_2ml = fm_stable( 2.0_rp,l(i,j) )
530  fhs_2ml = fh_stable( 2.0_rp,l(i,j) )
531 
532  ! select unstable / stable
533  sw = 0.5_rp - sign(0.5_rp,l(i,j)) ! if unstable, sw = 1
534 
535  fm(i,j) = max(eps, &
536  ( sw ) * fmus_z1l &
537  + ( 1.0_rp-sw ) * fms_z1l )
538  fh(i,j) = ( sw ) * fhus_z1l &
539  + ( 1.0_rp-sw ) * fhs_z1l
540  fm_10m(i,j) = max(eps, &
541  ( sw ) * fmus_10ml &
542  + ( 1.0_rp-sw ) * fms_10ml )
543  fm_2m(i,j) = max(eps, &
544  ( sw ) * fmus_2ml &
545  + ( 1.0_rp-sw ) * fms_2ml )
546  fh_2m(i,j) = ( sw ) * fhus_2ml &
547  + ( 1.0_rp-sw ) * fhs_2ml
548  enddo
549  enddo
550 
551  do j = js, je
552  do i = is, ie
553  r10m(i,j) = log( 10.0_rp/z0m(i,j) ) / logz1z0m(i,j) &
554  / sqrt( fm_10m(i,j) / fm(i,j) )
555  enddo
556  enddo
557 
558  do j = js, je
559  do i = is, ie
560  psi(i,j) = sqrt(fm(i,j)) / fh(i,j) * logz1z0m(i,j)
561  psi_2m(i,j) = sqrt(fm_2m(i,j)) / fh_2m(i,j) * log( 2.0_rp/z0m(i,j) )
562 
563  tmp1 = log( z0m(i,j)/z0h(i,j) ) / logz1z0m(i,j)
564 
565  r2h(i,j) = ( tmp1 + psi_2m(i,j) ) &
566  / ( tmp1 + psi(i,j) )
567 
568  tmp2 = log( z0m(i,j)/z0e(i,j) ) / logz1z0m(i,j)
569 
570  r2e(i,j) = ( tmp2 + psi_2m(i,j) ) &
571  / ( tmp2 + psi(i,j) )
572  enddo
573  enddo
574 
575 
576  if ( ite > itelim ) then
577  if( io_l ) write(io_fid_log,*) 'Error: reach maximum iteration in the function of ATMOS_PHY_SF_bulkcoef_beljaars.'
578  endif
579 
580  do j = js, je
581  do i = is, ie
582  cm(i,j) = min( max( cm(i,j), atmos_phy_sf_bulkcoef_cm_min ), atmos_phy_sf_bulkcoef_cm_max )
583  ch(i,j) = min( max( ch(i,j), atmos_phy_sf_bulkcoef_ch_min ), atmos_phy_sf_bulkcoef_ch_max )
584  ce(i,j) = min( max( ce(i,j), atmos_phy_sf_bulkcoef_ce_min ), atmos_phy_sf_bulkcoef_ce_max )
585  enddo
586  enddo
587 
588  do j = js, je
589  do i = is, ie
590  r10m(i,j) = min( max( r10m(i,j), 0.0_rp ), 1.0_rp )
591  r2h(i,j) = min( max( r2h(i,j), 0.0_rp ), 1.0_rp )
592  r2e(i,j) = min( max( r2e(i,j), 0.0_rp ), 1.0_rp )
593  enddo
594  enddo
595 
596  return
597  end subroutine atmos_phy_sf_bulkcoef_beljaars
598 
599  !-----------------------------------------------------------------------------
600  ! stability function for momemtum in unstable condition
601  function fm_unstable( Z, L ) result(fmUS)
602  use scale_const, only: &
603  pi => const_pi, &
604  eps => const_eps
605  implicit none
606 
607  real(RP), intent(in) :: Z
608  real(RP), intent(in) :: L
609  real(RP) :: fmUS
610 
611  real(RP) :: R, sqsqR
612  !---------------------------------------------------------------------------
613 
614  r = z / min(l,-eps) ! should be negative
615  sqsqr = ( 1.0_rp - 16.0_rp * r )**0.25_rp
616 
617  ! Paulson (1974); Dyer (1974)
618  fmus = log( ( 1.0_rp+sqsqr )**2 * ( 1.0_rp+sqsqr*sqsqr ) * 0.125_rp ) &
619  - 2.0_rp * atan( sqsqr ) &
620  + 0.5_rp * pi
621 
622  return
623  end function fm_unstable
624 
625  !-----------------------------------------------------------------------------
626  ! stability function for heat/vapor in unstable condition
627  function fh_unstable( Z, L ) result(fhUS)
628  use scale_const, only: &
629  eps => const_eps
630  implicit none
631 
632  real(RP), intent(in) :: Z
633  real(RP), intent(in) :: L
634  real(RP) :: fhUS
635 
636  real(RP) :: R, sqR
637  !---------------------------------------------------------------------------
638 
639  r = z / min(l,-eps) ! should be negative
640  sqr = sqrt( 1.0_rp - 16.0_rp * r )
641 
642  ! Paulson (1974); Dyer (1974)
643  fhus = 2.0_rp * log( ( 1.0_rp+sqr ) * 0.5_rp )
644 
645  return
646  end function fh_unstable
647 
648  !-----------------------------------------------------------------------------
649  ! stability function for momemtum in stable condition
650  function fm_stable( Z, L ) result(fmS)
651  use scale_const, only: &
652  eps => const_eps
653  implicit none
654 
655  real(RP), intent(in) :: Z
656  real(RP), intent(in) :: L
657  real(RP) :: fmS
658 
659  ! parameters of stability functions (Beljaars and Holtslag 1991)
660  real(RP), parameter :: a = 1.0_rp
661  real(RP), parameter :: b = 0.667_rp
662  real(RP), parameter :: c = 5.0_rp
663  real(RP), parameter :: d = 0.35_rp
664 
665  real(RP) :: R
666  !---------------------------------------------------------------------------
667 
668  r = z / max(l,eps) ! should be positive
669 
670  ! Holtslag and DeBruin (1988)
671  fms = - a*r - b*( r - c/d )*exp( -d*r ) - b*c/d
672 
673  return
674  end function fm_stable
675 
676  !-----------------------------------------------------------------------------
677  ! stability function for heat/vapor in stable condition
678  function fh_stable( Z, L ) result(fhS)
679  use scale_const, only: &
680  eps => const_eps
681  implicit none
682 
683  real(RP), intent(in) :: Z
684  real(RP), intent(in) :: L
685  real(RP) :: fhS
686 
687  ! parameters of stability functions (Beljaars and Holtslag 1991)
688  real(RP), parameter :: a = 1.0_rp
689  real(RP), parameter :: b = 0.667_rp
690  real(RP), parameter :: c = 5.0_rp
691  real(RP), parameter :: d = 0.35_rp
692 
693  real(RP) :: R
694  !---------------------------------------------------------------------------
695 
696  r = z / max(l,eps) ! should be positive
697 
698  ! Beljaars and Holtslag (1991)
699  fhs = 1.0_rp - ( 1.0_rp + 2.0_rp/3.0_rp * a*r )**1.5_rp - b*( r - c/d )*exp( -d*r ) - b*c/d
700 
701  return
702  end function fh_stable
703 
704 end module scale_atmos_phy_sf_bulkcoef
scale_grid_index::is
integer, public is
start point of inner domain: x, local
Definition: scale_grid_index.F90:60
scale_grid_index::je
integer, public je
end point of inner domain: y, local
Definition: scale_grid_index.F90:63
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_atmos_phy_sf_bulkcoef::atmos_phy_sf_bulkcoef_setup
subroutine, public atmos_phy_sf_bulkcoef_setup
Definition: scale_atmos_phy_sf_bulkcoef.F90:98
scale_const::const_karman
real(rp), parameter, public const_karman
von Karman constant
Definition: scale_const.F90:52
scale_grid_index
module grid index
Definition: scale_grid_index.F90:14
scale_tracer
module TRACER
Definition: scale_tracer.F90:18
scale_grid_index::ia
integer, public ia
of x whole cells (local, with HALO)
Definition: scale_grid_index.F90:55
scale_const::const_grav
real(rp), public const_grav
standard acceleration of gravity [m/s2]
Definition: scale_const.F90:48
scale_grid_index::js
integer, public js
start point of inner domain: y, local
Definition: scale_grid_index.F90:62
scale_atmos_phy_sf_bulkcoef::atmos_phy_sf_bulkcoef
procedure(bc), pointer, public atmos_phy_sf_bulkcoef
Definition: scale_atmos_phy_sf_bulkcoef.F90:62
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_prof
module profiler
Definition: scale_prof.F90:10
scale_grid_index::ie
integer, public ie
end point of inner domain: x, local
Definition: scale_grid_index.F90:61
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_const::const_pi
real(rp), public const_pi
pi
Definition: scale_const.F90:34
scale_stdio::io_fid_conf
integer, public io_fid_conf
Config file ID.
Definition: scale_stdio.F90:55
scale_atmos_phy_sf_bulkcoef
module ATMOSPHERE / Physics Surface bulk coefficient
Definition: scale_atmos_phy_sf_bulkcoef.F90:15
scale_stdio::io_fid_log
integer, public io_fid_log
Log file ID.
Definition: scale_stdio.F90:56
scale_grid_index::ja
integer, public ja
of y whole cells (local, with HALO)
Definition: scale_grid_index.F90:56
Generated by   doxygen 1.8.13