SCALE-RM
scale_urban_dyn_kusaka01.F90
Go to the documentation of this file.
1 !-------------------------------------------------------------------------------
10 !-------------------------------------------------------------------------------
11 #include "scalelib.h"
12 module scale_urban_dyn_kusaka01
13  !-----------------------------------------------------------------------------
14  !
15  !++ used modules
16  !
17  use scale_precision
18  use scale_io
19  use scale_prof
20  use scale_cpl_sfc_index
21  !-----------------------------------------------------------------------------
22  implicit none
23  private
24  !-----------------------------------------------------------------------------
25  !
26  !++ Public procedure
27  !
28  public :: urban_dyn_kusaka01_setup
29  public :: urban_dyn_kusaka01
30 
31  !-----------------------------------------------------------------------------
32  !
33  !++ Public parameters & variables
34  !
35  !-----------------------------------------------------------------------------
36  !
37  !++ Private procedure
38  !
39  private :: slc_main
40  private :: canopy_wind
41  private :: cal_beta
42  private :: cal_psi
43  private :: mos
44  private :: multi_layer
45  private :: urban_param_setup
46 
47  !-----------------------------------------------------------------------------
48  !
49  !++ Private parameters & variables
50  !
51  !-----------------------------------------------------------------------------
52  !
53  ! from namelist
54  real(RP), private :: dts_max = 0.1_rp ! maximum dT during one minute [K/sec]
55  ! 0.1 [K/sec] = 6.0 [K/min]
56  real(RP), private :: zr = 10.0_rp ! roof level ( building height) [m]
57  real(RP), private :: roof_width = 9.0_rp ! roof level ( building height) [m]
58  real(RP), private :: road_width = 11.0_rp ! roof level ( building height) [m]
59  real(RP), private :: sigma_zed = 1.0_rp ! Standard deviation of roof height [m]
60  real(RP), private :: ah = 17.5_rp ! Sensible Anthropogenic heat [W/m^2]
61  real(RP), private :: alh = 0.0_rp ! Latent Anthropogenic heat [W/m^2]
62  real(RP), private :: betr = 0.0_rp ! Evaporation efficiency of roof [-]
63  real(RP), private :: betb = 0.0_rp ! of building [-]
64  real(RP), private :: betg = 0.0_rp ! of ground [-]
65  real(RP), private :: strgr = 0.0_rp ! rain strage on roof [-]
66  real(RP), private :: strgb = 0.0_rp ! on wall [-]
67  real(RP), private :: strgg = 0.0_rp ! on ground [-]
68  real(RP), private :: capr = 1.2e6_rp ! heat capacity of roof [J m-3 K]
69  real(RP), private :: capb = 1.2e6_rp ! of wall [J m-3 K]
70  real(RP), private :: capg = 1.2e6_rp ! of ground [J m-3 K]
71  real(RP), private :: aksr = 2.28_rp ! thermal conductivity of roof [W m-1 K]
72  real(RP), private :: aksb = 2.28_rp ! of wall [W m-1 K]
73  real(RP), private :: aksg = 2.28_rp ! of ground [W m-1 K]
74  real(RP), private :: albr = 0.2_rp ! surface albedo of roof
75  real(RP), private :: albb = 0.2_rp ! surface albedo of wall
76  real(RP), private :: albg = 0.2_rp ! surface albedo of ground
77  real(RP), private :: epsr = 0.90_rp ! Surface emissivity of roof
78  real(RP), private :: epsb = 0.90_rp ! Surface emissivity of wall
79  real(RP), private :: epsg = 0.90_rp ! Surface emissivity of ground
80  real(RP), private :: z0r = 0.01_rp ! roughness length for momentum of building roof
81  real(RP), private :: z0b = 0.0001_rp ! roughness length for momentum of building wall
82  real(RP), private :: z0g = 0.01_rp ! roughness length for momentum of ground
83  real(RP), private :: trlend = 293.00_rp ! lower boundary condition of roof temperature [K]
84  real(RP), private :: tblend = 293.00_rp ! lower boundary condition of wall temperature [K]
85  real(RP), private :: tglend = 293.00_rp ! lower boundary condition of ground temperature [K]
86  integer , private :: bound = 1 ! Boundary Condition for Roof, Wall, Ground Layer Temp
87  ! [1: Zero-Flux, 2: T = Constant]
88  real(RP), private :: ahdiurnal(1:24) ! AH diurnal profile
89 
90  ! calculated in subroutine urban_param_set
91  real(RP), private :: r ! Normalized roof wight (eq. building coverage ratio)
92  real(RP), private :: rw ! (= 1 - R)
93  real(RP), private :: hgt ! Normalized building height
94  real(RP), private :: z0hr ! roughness length for heat of roof
95  real(RP), private :: z0hb ! roughness length for heat of building wall
96  real(RP), private :: z0hg ! roughness length for heat of ground
97  real(RP), private :: z0c ! Roughness length above canyon for momentum [m]
98  real(RP), private :: z0hc ! Roughness length above canyon for heat [m]
99  real(RP), private :: zdc ! Displacement height [m]
100  real(RP), private :: svf ! Sky view factor [-]
101 
102  real(RP), private :: xxxr = 0.0_rp ! Monin-Obkhov length for roof [-]
103  real(RP), private :: xxxc = 0.0_rp ! Monin-Obkhov length for canopy [-]
104 
105  ! history
106  integer, private :: i_shr, i_shb, i_shg
107  integer, private :: i_lhr, i_lhb, i_lhg
108  integer, private :: i_ghr, i_ghb, i_ghg
109  integer, private :: i_rnr, i_rnb, i_rng, i_rngrd
110 
111  !-----------------------------------------------------------------------------
112 contains
113  !-----------------------------------------------------------------------------
115  subroutine urban_dyn_kusaka01_setup( &
116  UIA, UIS, UIE, UJA, UJS, UJE, &
117  Z0M, Z0H, Z0E )
118  use scale_prc, only: &
119  prc_abort
120  use scale_file_history, only: &
121  file_history_reg
122  implicit none
123 
124  integer, intent(in) :: UIA, UIS, UIE
125  integer, intent(in) :: UJA, UJS, UJE
126  real(RP), intent(out) :: Z0M(uia,uja)
127  real(RP), intent(out) :: Z0H(uia,uja)
128  real(RP), intent(out) :: Z0E(uia,uja)
129 
130  namelist / param_urban_dyn_kusaka01 / &
131  dts_max, &
132  zr, &
133  roof_width, &
134  road_width, &
135  sigma_zed, &
136  ah, &
137  alh, &
138  betr, &
139  betb, &
140  betg, &
141  strgr, &
142  strgb, &
143  strgg, &
144  capr, &
145  capb, &
146  capg, &
147  aksr, &
148  aksb, &
149  aksg, &
150  albr, &
151  albb, &
152  albg, &
153  epsr, &
154  epsb, &
155  epsg, &
156  z0r, &
157  z0b, &
158  z0g, &
159  trlend, &
160  tblend, &
161  tglend, &
162  bound
163 
164  integer :: i, j
165  integer :: ierr
166  !---------------------------------------------------------------------------
167 
168  log_newline
169  log_info("URBAN_DYN_kusaka01_setup",*) 'Setup'
170 
171  !--- read namelist
172  rewind(io_fid_conf)
173  read(io_fid_conf,nml=param_urban_dyn_kusaka01,iostat=ierr)
174  if( ierr < 0 ) then !--- missing
175  log_info("URBAN_DYN_kusaka01_setup",*) 'Not found namelist. Default used.'
176  elseif( ierr > 0 ) then !--- fatal error
177  log_error("URBAN_DYN_kusaka01_setup",*) 'Not appropriate names in namelist PARAM_URBAN_DYN_KUSAKA01. Check!'
178  call prc_abort
179  endif
180  log_nml(param_urban_dyn_kusaka01)
181 
182  ahdiurnal(:) = (/ 0.356, 0.274, 0.232, 0.251, 0.375, 0.647, 0.919, 1.135, 1.249, 1.328, &
183  1.365, 1.363, 1.375, 1.404, 1.457, 1.526, 1.557, 1.521, 1.372, 1.206, &
184  1.017, 0.876, 0.684, 0.512 /)
185 
186  ! set other urban parameters
187  call urban_param_setup
188 
189  ! judge to run slab land model
190  do j = ujs, uje
191  do i = uis, uie
192  z0m(i,j) = z0c
193  z0h(i,j) = z0hc
194  z0e(i,j) = z0hc
195  enddo
196  enddo
197 
198  call file_history_reg( 'URBAN_SHR', 'urban sensible heat flux on roof', 'W/m2', i_shr , ndims=2 )
199  call file_history_reg( 'URBAN_SHB', 'urban sensible heat flux on wall', 'W/m2', i_shb , ndims=2 )
200  call file_history_reg( 'URBAN_SHG', 'urban sensible heat flux on road', 'W/m2', i_shg , ndims=2 )
201  call file_history_reg( 'URBAN_LHR', 'urban latent heat flux on roof', 'W/m2', i_lhr , ndims=2 )
202  call file_history_reg( 'URBAN_LHB', 'urban latent heat flux on wall', 'W/m2', i_lhb , ndims=2 )
203  call file_history_reg( 'URBAN_LHG', 'urban latent heat flux on road', 'W/m2', i_lhg , ndims=2 )
204  call file_history_reg( 'URBAN_GHR', 'urban ground heat flux on roof', 'W/m2', i_ghr , ndims=2 )
205  call file_history_reg( 'URBAN_GHB', 'urban ground heat flux on wall', 'W/m2', i_ghb , ndims=2 )
206  call file_history_reg( 'URBAN_GHG', 'urban ground heat flux on road', 'W/m2', i_ghg , ndims=2 )
207  call file_history_reg( 'URBAN_RNR', 'urban net radiation on roof', 'W/m2', i_rnr , ndims=2 )
208  call file_history_reg( 'URBAN_RNB', 'urban net radiation on wall', 'W/m2', i_rnb , ndims=2 )
209  call file_history_reg( 'URBAN_RNG', 'urban net radiation on road', 'W/m2', i_rng , ndims=2 )
210  call file_history_reg( 'URBAN_RNgrd', 'urban grid average of net radiation', 'W/m2', i_rngrd, ndims=2 )
211 
212  return
213  end subroutine urban_dyn_kusaka01_setup
214 
215  !-----------------------------------------------------------------------------
217 
218  subroutine urban_dyn_kusaka01( &
219  UKA, UKS, UKE, UIA, UIS, UIE, UJA, UJS, UJE, &
220  TMPA, PRSA, &
221  W1, U1, V1, &
222  DENS, QA, LHV, &
223  Z1, PBL, &
224  RHOS, PRSS, &
225  LWD, SWD, &
226  RAIN, SNOW, &
227  CDZ, &
228  fact_urban, &
229  tloc, dsec, dt, &
230  TRL_URB, TBL_URB, TGL_URB, &
231  TR_URB, TB_URB, TG_URB, &
232  TC_URB, QC_URB, UC_URB, &
233  RAINR_URB, RAINB_URB, RAING_URB, &
234  ROFF_URB, &
235  SFC_TEMP, &
236  ALBEDO, &
237  MWFLX, MUFLX, MVFLX, &
238  SHFLX, LHFLX, GHFLX, &
239  Z0M, Z0H, Z0E, &
240  U10, V10, T2, Q2 )
241  use scale_const, only: &
242  rdry => const_rdry, &
243  rvap => const_rvap
244  use scale_atmos_saturation, only: &
245  qsat => atmos_saturation_pres2qsat_all
246  use scale_bulkflux, only: &
247  bulkflux
248  implicit none
249  integer, intent(in) :: UKA, UKS, UKE
250  integer, intent(in) :: UIA, UIS, UIE
251  integer, intent(in) :: UJA, UJS, UJE
252 
253  real(RP), intent(in) :: TMPA(uia,uja)
254  real(RP), intent(in) :: PRSA(uia,uja)
255  real(RP), intent(in) :: W1 (uia,uja)
256  real(RP), intent(in) :: U1 (uia,uja)
257  real(RP), intent(in) :: V1 (uia,uja)
258  real(RP), intent(in) :: DENS(uia,uja)
259  real(RP), intent(in) :: QA (uia,uja)
260  real(RP), intent(in) :: LHV (uia,uja)
261  real(RP), intent(in) :: Z1 (uia,uja)
262  real(RP), intent(in) :: PBL (uia,uja)
263  real(RP), intent(in) :: RHOS(uia,uja) ! density at the surface [kg/m3]
264  real(RP), intent(in) :: PRSS(uia,uja)
265  real(RP), intent(in) :: LWD (uia,uja,2)
266  real(RP), intent(in) :: SWD (uia,uja,2)
267  real(RP), intent(in) :: RAIN(uia,uja)
268  real(RP), intent(in) :: SNOW(uia,uja)
269  real(RP), intent(in) :: CDZ(uka)
270  real(RP), intent(in) :: fact_urban(uia,uja)
271  integer, intent(in) :: tloc
272  real(RP), intent(in) :: dsec
273  real(DP), intent(in) :: dt
274 
275  real(RP), intent(inout) :: TR_URB (uia,uja)
276  real(RP), intent(inout) :: TB_URB (uia,uja)
277  real(RP), intent(inout) :: TG_URB (uia,uja)
278  real(RP), intent(inout) :: TC_URB (uia,uja)
279  real(RP), intent(inout) :: QC_URB (uia,uja)
280  real(RP), intent(inout) :: UC_URB (uia,uja)
281  real(RP), intent(inout) :: TRL_URB (uks:uke,uia,uja)
282  real(RP), intent(inout) :: TBL_URB (uks:uke,uia,uja)
283  real(RP), intent(inout) :: TGL_URB (uks:uke,uia,uja)
284  real(RP), intent(inout) :: RAINR_URB(uia,uja)
285  real(RP), intent(inout) :: RAINB_URB(uia,uja)
286  real(RP), intent(inout) :: RAING_URB(uia,uja)
287  real(RP), intent(inout) :: ROFF_URB (uia,uja)
288 
289  real(RP), intent(out) :: SFC_TEMP(uia,uja)
290  real(RP), intent(out) :: ALBEDO (uia,uja,n_rad_dir,n_rad_rgn)
291  real(RP), intent(out) :: MWFLX (uia,uja)
292  real(RP), intent(out) :: MUFLX (uia,uja)
293  real(RP), intent(out) :: MVFLX (uia,uja)
294  real(RP), intent(out) :: SHFLX (uia,uja)
295  real(RP), intent(out) :: LHFLX (uia,uja)
296  real(RP), intent(out) :: GHFLX (uia,uja)
297  real(RP), intent(out) :: Z0M (uia,uja)
298  real(RP), intent(out) :: Z0H (uia,uja)
299  real(RP), intent(out) :: Z0E (uia,uja)
300  real(RP), intent(out) :: U10 (uia,uja)
301  real(RP), intent(out) :: V10 (uia,uja)
302  real(RP), intent(out) :: T2 (uia,uja)
303  real(RP), intent(out) :: Q2 (uia,uja)
304 
305 
306  ! parameter
307  logical, parameter :: LSOLAR = .false. ! [true=both, false=SSG only]
308 
309  real(RP), parameter :: Uabs_min = 0.1_rp
310 
311  ! work
312  real(RP) :: TR
313  real(RP) :: TB
314  real(RP) :: TG
315  real(RP) :: TC
316  real(RP) :: QC
317  real(RP) :: UC
318  real(RP) :: TRL(uks:uke)
319  real(RP) :: TBL(uks:uke)
320  real(RP) :: TGL(uks:uke)
321  real(RP) :: RAINR
322  real(RP) :: RAINB
323  real(RP) :: RAING
324  real(RP) :: ROFF
325  real(RP) :: ALBD_LW
326  real(RP) :: ALBD_SW
327 
328  real(RP) :: SHR (uia,uja)
329  real(RP) :: SHB (uia,uja)
330  real(RP) :: SHG (uia,uja)
331  real(RP) :: LHR (uia,uja)
332  real(RP) :: LHB (uia,uja)
333  real(RP) :: LHG (uia,uja)
334  real(RP) :: GHR (uia,uja)
335  real(RP) :: GHB (uia,uja)
336  real(RP) :: GHG (uia,uja)
337  real(RP) :: RNR (uia,uja)
338  real(RP) :: RNB (uia,uja)
339  real(RP) :: RNG (uia,uja)
340  real(RP) :: RNgrd(uia,uja)
341 
342  real(RP) :: DZR(uka) ! thickness of each roof layer [m]
343  real(RP) :: DZB(uka) ! thickness of each building layer [m]
344  real(RP) :: DZG(uka) ! thickness of each road layer [m]
345 
346 
347  real(RP) :: Ustar ! friction velocity [m]
348  real(RP) :: Tstar ! friction temperature [K]
349  real(RP) :: Qstar ! friction mixing rate [kg/kg]
350  real(RP) :: Uabs ! modified absolute velocity [m/s]
351  real(RP) :: Ra ! Aerodynamic resistance (=1/Ce) [1/s]
352 
353  real(RP) :: QVsat ! saturation water vapor mixing ratio at surface [kg/kg]
354  real(RP) :: Rtot ! total gas constant
355  real(RP) :: qdry ! dry air mass ratio [kg/kg]
356 
357  real(RP) :: FracU10 ! calculation parameter for U10 [-]
358  real(RP) :: FracT2 ! calculation parameter for T2 [-]
359  real(RP) :: FracQ2 ! calculation parameter for Q2 [-]
360 
361  integer :: k, i, j
362  !---------------------------------------------------------------------------
363 
364  log_progress(*) 'urban / dynamics / Kusaka01'
365 
366  dzr(:) = cdz(:)
367  dzb(:) = cdz(:)
368  dzg(:) = cdz(:)
369 
370  do j = ujs, uje
371  do i = uis, uie
372 
373  if( fact_urban(i,j) > 0.0_rp ) then
374 
375  qdry = 1.0_rp - qa(i,j)
376  rtot = qdry * rdry + qa(i,j) * rvap
377 
378  uabs = max( sqrt( u1(i,j)**2 + v1(i,j)**2 + w1(i,j)**2 ), uabs_min )
379 
380  ! save
381  tr = tr_urb(i,j)
382  tb = tb_urb(i,j)
383  tg = tg_urb(i,j)
384  tc = tc_urb(i,j)
385  qc = qc_urb(i,j)
386  uc = uc_urb(i,j)
387 
388  do k = uks, uke
389  trl(k) = trl_urb(k,i,j)
390  tbl(k) = tbl_urb(k,i,j)
391  tgl(k) = tgl_urb(k,i,j)
392  end do
393 
394  rainr = rainr_urb(i,j)
395  rainb = rainb_urb(i,j)
396  raing = raing_urb(i,j)
397  roff = roff_urb(i,j)
398 
399  call slc_main( uka, uks, uke, uia, uis, uie, uja, ujs, uje, &
400  trl(:), & ! [INOUT]
401  tbl(:), & ! [INOUT]
402  tgl(:), & ! [INOUT]
403  tr, & ! [INOUT]
404  tb, & ! [INOUT]
405  tg, & ! [INOUT]
406  tc, & ! [INOUT]
407  qc, & ! [INOUT]
408  uc, & ! [INOUT]
409  rainr, & ! [INOUT]
410  rainb, & ! [INOUT]
411  raing, & ! [INOUT]
412  roff, & ! [INOUT]
413  albd_lw, & ! [OUT]
414  albd_sw, & ! [OUT]
415  shr(i,j), & ! [OUT]
416  shb(i,j), & ! [OUT]
417  shg(i,j), & ! [OUT]
418  lhr(i,j), & ! [OUT]
419  lhb(i,j), & ! [OUT]
420  lhg(i,j), & ! [OUT]
421  ghr(i,j), & ! [OUT]
422  ghb(i,j), & ! [OUT]
423  ghg(i,j), & ! [OUT]
424  rnr(i,j), & ! [OUT]
425  rnb(i,j), & ! [OUT]
426  rng(i,j), & ! [OUT]
427  sfc_temp(i,j), & ! [OUT]
428  rngrd(i,j), & ! [OUT]
429  shflx(i,j), & ! [OUT]
430  lhflx(i,j), & ! [OUT]
431  ghflx(i,j), & ! [OUT]
432  u10(i,j), & ! [OUT]
433  v10(i,j), & ! [OUT]
434  t2(i,j), & ! [OUT]
435  q2(i,j), & ! [OUT]
436  lsolar, & ! [IN]
437  prsa(i,j), & ! [IN]
438  prss(i,j), & ! [IN]
439  tmpa(i,j), & ! [IN]
440  qa(i,j), & ! [IN]
441  uabs, & ! [IN]
442  u1(i,j), & ! [IN]
443  v1(i,j), & ! [IN]
444  lhv(i,j), & ! [IN]
445  z1(i,j), & ! [IN]
446  swd(i,j,:), & ! [IN]
447  lwd(i,j,:), & ! [IN]
448  rain(i,j), & ! [IN]
449  snow(i,j), & ! [IN]
450  dens(i,j), & ! [IN]
451  dzr(:), dzg(:), dzb(:), & ! [IN]
452  tloc, dsec, dt, & ! (in)
453  i, j ) ! [IN]
454 
455 
456  ! update
457  tr_urb(i,j) = tr
458  tb_urb(i,j) = tb
459  tg_urb(i,j) = tg
460  tc_urb(i,j) = tc
461  qc_urb(i,j) = qc
462  uc_urb(i,j) = uc
463  do k = uks, uke
464  trl_urb(k,i,j) = trl(k)
465  tbl_urb(k,i,j) = tbl(k)
466  tgl_urb(k,i,j) = tgl(k)
467  end do
468  rainr_urb(i,j) = rainr
469  rainb_urb(i,j) = rainb
470  raing_urb(i,j) = raing
471  roff_urb(i,j) = roff
472 
473  albedo(i,j,i_r_direct ,i_r_ir ) = albd_lw
474  albedo(i,j,i_r_diffuse,i_r_ir ) = albd_lw
475  albedo(i,j,i_r_direct ,i_r_nir) = albd_sw
476  albedo(i,j,i_r_diffuse,i_r_nir) = albd_sw
477  albedo(i,j,i_r_direct ,i_r_vis) = albd_sw
478  albedo(i,j,i_r_diffuse,i_r_vis) = albd_sw
479 
480  ! saturation at the surface
481  call qsat( sfc_temp(i,j), prss(i,j), qdry, & ! [IN]
482  qvsat ) ! [OUT]
483 
484  call bulkflux( ustar, & ! [OUT]
485  tstar, & ! [OUT]
486  qstar, & ! [OUT]
487  uabs, & ! [OUT]
488  ra, & ! [OUT]
489  fracu10, & ! [OUT]
490  fract2, & ! [OUT]
491  fracq2, & ! [OUT]
492  tmpa(i,j), & ! [IN]
493  sfc_temp(i,j), & ! [IN]
494  prsa(i,j), & ! [IN]
495  prss(i,j), & ! [IN]
496  qa(i,j), & ! [IN]
497  qvsat, & ! [IN]
498  u1(i,j), & ! [IN]
499  v1(i,j), & ! [IN]
500  z1(i,j), & ! [IN]
501  pbl(i,j), & ! [IN]
502  z0c, & ! [IN]
503  z0hc, & ! [IN]
504  z0hc ) ! [IN]
505 
506  mwflx(i,j) = -rhos(i,j) * ustar**2 / uabs * w1(i,j)
507  muflx(i,j) = -rhos(i,j) * ustar**2 / uabs * u1(i,j)
508  mvflx(i,j) = -rhos(i,j) * ustar**2 / uabs * v1(i,j)
509 
510  z0m(i,j) = z0c
511  z0h(i,j) = z0hc
512  z0e(i,j) = z0hc
513 
514  else
515  sfc_temp(i,j) = 300.0_rp ! constant value
516  albedo(i,j,:,:) = 0.0_rp
517  mwflx(i,j) = 0.0_rp
518  muflx(i,j) = 0.0_rp
519  mvflx(i,j) = 0.0_rp
520  shflx(i,j) = 0.0_rp
521  lhflx(i,j) = 0.0_rp
522  ghflx(i,j) = 0.0_rp
523  z0m(i,j) = 0.0_rp
524  z0h(i,j) = 0.0_rp
525  z0e(i,j) = 0.0_rp
526  u10(i,j) = 0.0_rp
527  v10(i,j) = 0.0_rp
528  t2(i,j) = 0.0_rp
529  q2(i,j) = 0.0_rp
530  shr(i,j) = 0.0_rp
531  shb(i,j) = 0.0_rp
532  shg(i,j) = 0.0_rp
533  lhr(i,j) = 0.0_rp
534  lhb(i,j) = 0.0_rp
535  lhg(i,j) = 0.0_rp
536  ghr(i,j) = 0.0_rp
537  ghb(i,j) = 0.0_rp
538  ghg(i,j) = 0.0_rp
539  rnr(i,j) = 0.0_rp
540  rnb(i,j) = 0.0_rp
541  rng(i,j) = 0.0_rp
542  rngrd(i,j) = 0.0_rp
543  endif
544 
545  end do
546  end do
547 
548  call put_history( uia, uja, &
549  shr(:,:), shb(:,:), shg(:,:), &
550  lhr(:,:), lhb(:,:), lhg(:,:), &
551  ghr(:,:), ghb(:,:), ghg(:,:), &
552  rnr(:,:), rnb(:,:), rng(:,:), &
553  rngrd(:,:) )
554 
555  return
556  end subroutine urban_dyn_kusaka01
557 
558  !-----------------------------------------------------------------------------
559  subroutine slc_main( &
560  UKA, UKS, UKE, UIA, UIS, UIE, UJA, UJS, UJE, &
561  TRL, & ! (inout)
562  tbl, & ! (inout)
563  tgl, & ! (inout)
564  tr, & ! (inout)
565  tb, & ! (inout)
566  tg, & ! (inout)
567  tc, & ! (inout)
568  qc, & ! (inout)
569  uc, & ! (inout)
570  rainr, & ! (inout)
571  rainb, & ! (inout)
572  raing, & ! (inout)
573  roff, & ! (inout)
574  albd_lw_grid, & ! (out)
575  albd_sw_grid, & ! (out)
576  shr, & ! (out)
577  shb, & ! (out)
578  shg, & ! (out)
579  lhr, & ! (out)
580  lhb, & ! (out)
581  lhg, & ! (out)
582  ghr, & ! (out)
583  ghb, & ! (out)
584  ghg, & ! (out)
585  rnr, & ! (out)
586  rnb, & ! (out)
587  rng, & ! (out)
588  rts, & ! (out)
589  rn, & ! (out)
590  sh, & ! (out)
591  lh, & ! (out)
592  ghflx, & ! (out)
593  u10, & ! (out)
594  v10, & ! (out)
595  t2, & ! (out)
596  q2, & ! (out)
597  lsolar, & ! (in)
598  prsa, & ! (in)
599  prss, & ! (in)
600  ta, & ! (in)
601  qa, & ! (in)
602  ua, & ! (in)
603  u1, & ! (in)
604  v1, & ! (in)
605  lhv, & ! (in)
606  za, & ! (in)
607  ssg, & ! (in)
608  llg, & ! (in)
609  rain, & ! (in)
610  snow, & ! (in)
611  rhoo, & ! (in)
612  dzr, dzb, dzg, & ! (in)
613  tloc, dsec, & ! (in)
614  dt, & ! (in)
615  i, j ) ! (in)
616  use scale_prc, only: &
617  prc_myrank, &
618  prc_abort
619  use scale_const, only: &
620  eps => const_eps, & ! small number (machine epsilon)
621  karman => const_karman, & ! Kalman constant [-]
622  cpdry => const_cpdry, & ! Heat capacity of dry air [J/K/kg]
623  grav => const_grav, & ! gravitational constant [m/s2]
624  rdry => const_rdry, & ! specific gas constant (dry) [J/kg/K]
625  rvap => const_rvap, & ! gas constant (water vapor) [J/kg/K]
626  stb => const_stb, & ! stefan-Boltzman constant [MKS unit]
627  tem00 => const_tem00, & ! temperature reference (0 degree C) [K]
628  pre00 => const_pre00 ! pressure reference [Pa]
629  use scale_atmos_hydrometeor, only: &
630  hydrometeor_lhv => atmos_hydrometeor_lhv
631  use scale_atmos_saturation, only: &
632  qsat => atmos_saturation_pres2qsat_all
633  implicit none
634 
635  integer, intent(in) :: UKA, UKS, UKE
636  integer, intent(in) :: UIA, UIS, UIE
637  integer, intent(in) :: UJA, UJS, UJE
638 
639  !-- configuration variables
640  logical , intent(in) :: LSOLAR ! logical [true=both, false=SSG only]
641 
642  !-- Input variables from Coupler to Urban
643  real(RP), intent(in) :: PRSA ! Pressure at 1st atmospheric layer [Pa]
644  real(RP), intent(in) :: PRSS ! Surface Pressure [Pa]
645  real(RP), intent(in) :: TA ! temp at 1st atmospheric level [K]
646  real(RP), intent(in) :: QA ! specific humidity at 1st atmospheric level [kg/kg]
647  real(RP), intent(in) :: UA ! wind speed at 1st atmospheric level [m/s]
648  real(RP), intent(in) :: U1 ! u at 1st atmospheric level [m/s]
649  real(RP), intent(in) :: V1 ! v at 1st atmospheric level [m/s]
650  real(RP), intent(in) :: LHV ! latent heat of vaporization [J/kg]
651  real(RP), intent(in) :: ZA ! height of 1st atmospheric level [m]
652  real(RP), intent(in) :: SSG(2) ! downward total short wave radiation [W/m/m]
653  real(RP), intent(in) :: LLG(2) ! downward long wave radiation [W/m/m]
654  real(RP), intent(in) :: RAIN ! liquid water flux [kg/m2/s]
655  real(RP), intent(in) :: SNOW ! ice water flux [kg/m2/s]
656  real(RP), intent(in) :: RHOO ! air density [kg/m^3]
657  real(RP), intent(in) :: DZR(uka)
658  real(RP), intent(in) :: DZB(uka)
659  real(RP), intent(in) :: DZG(uka)
660  integer, intent(in) :: tloc
661  real(RP), intent(in) :: dsec
662  real(DP), intent(in) :: dt
663 
664  !-- In/Out variables from/to Coupler to/from Urban
665  real(RP), intent(inout) :: TRL(uks:uke) ! layer temperature [K]
666  real(RP), intent(inout) :: TBL(uks:uke) ! layer temperature [K]
667  real(RP), intent(inout) :: TGL(uks:uke) ! layer temperature [K]
668  real(RP), intent(inout) :: TR ! roof temperature [K]
669  real(RP), intent(inout) :: TB ! building wall temperature [K]
670  real(RP), intent(inout) :: TG ! road temperature [K]
671  real(RP), intent(inout) :: TC ! urban-canopy air temperature [K]
672  real(RP), intent(inout) :: QC ! urban-canopy air specific humidity [kg/kg]
673  real(RP), intent(inout) :: UC ! diagnostic canopy wind [m/s]
674  real(RP), intent(inout) :: RAINR ! rain amount in storage on roof [kg/m2]
675  real(RP), intent(inout) :: RAINB ! rain amount in storage on building [kg/m2]
676  real(RP), intent(inout) :: RAING ! rain amount in storage on road [kg/m2]
677  real(RP), intent(inout) :: ROFF ! runoff from urban [kg/m2]
678 
679  !-- Output variables from Urban to Coupler
680  real(RP), intent(out) :: ALBD_SW_grid ! grid mean of surface albedo for SW
681  real(RP), intent(out) :: ALBD_LW_grid ! grid mean of surface albedo for LW ( 1-emiss )
682  real(RP), intent(out) :: RTS ! radiative surface temperature [K]
683  real(RP), intent(out) :: SH ! sensible heat flux [W/m/m]
684  real(RP), intent(out) :: LH ! latent heat flux [W/m/m]
685  real(RP), intent(out) :: GHFLX ! heat flux into the ground [W/m/m]
686  real(RP), intent(out) :: RN ! net radition [W/m/m]
687  real(RP), intent(out) :: U10 ! U wind at 10m [m/s]
688  real(RP), intent(out) :: V10 ! V wind at 10m [m/s]
689  real(RP), intent(out) :: T2 ! air temperature at 2m [K]
690  real(RP), intent(out) :: Q2 ! specific humidity at 2m [kg/kg]
691  real(RP), intent(out) :: RNR, RNB, RNG
692  real(RP), intent(out) :: SHR, SHB, SHG
693  real(RP), intent(out) :: LHR, LHB, LHG
694  real(RP), intent(out) :: GHR, GHB, GHG
695  integer , intent(in) :: i, j
696 
697  !-- parameters
698 ! real(RP), parameter :: SRATIO = 0.75_RP ! ratio between direct/total solar [-]
699 ! real(RP), parameter :: TFa = 0.5_RP ! factor a in Tomita (2009)
700 ! real(RP), parameter :: TFb = 1.1_RP ! factor b in Tomita (2009)
701  real(RP), parameter :: redf_min = 1.0e-2_rp ! minimum reduced factor
702  real(RP), parameter :: redf_max = 1.0_rp ! maximum reduced factor
703 ! real(RP), parameter :: CAP_water = 4.185E6_RP ! Heat capacity of water (15 deg) [J m-3 K]
704 ! real(RP), parameter :: AKS_water = 0.59_RP ! Thermal conductivity of water [W m-1 K]
705 
706  !-- Local variables
707 ! logical :: SHADOW = .false.
708  ! true = consider svf and shadow effects,
709  ! false = consider svf effect only
710 
711  real(RP) :: AH_t ! Sensible Anthropogenic heat [W/m^2]
712  real(RP) :: ALH_t ! Latent Anthropogenic heat [W/m^2]
713 
714  real(RP) :: SSGD ! downward direct short wave radiation [W/m/m]
715  real(RP) :: SSGQ ! downward diffuse short wave radiation [W/m/m]
716 
717  real(RP) :: W, VFGS, VFGW, VFWG, VFWS, VFWW
718  real(RP) :: SX, RX
719 
720  real(RP) :: TRP ! TRP: at previous time step [K]
721  real(RP) :: TBP ! TBP: at previous time step [K]
722  real(RP) :: TGP ! TGP: at previous time step [K]
723  real(RP) :: TCP ! TCP: at previous time step [K]
724  real(RP) :: QCP ! QCP: at previous time step [kg/kg]
725  real(RP) :: TRLP(uks:uke) ! Layer temperature at previous step [K]
726  real(RP) :: TBLP(uks:uke) ! Layer temperature at previous step [K]
727  real(RP) :: TGLP(uks:uke) ! Layer temperature at previous step [K]
728  !
729  real(RP) :: RAINRP ! at previous step, rain amount in storage on roof [kg/m2]
730  real(RP) :: RAINBP ! at previous step, rain amount in storage on building [kg/m2]
731  real(RP) :: RAINGP ! at previous step, rain amount in storage on road [kg/m2]
732 
733  !real(RP) :: UST, TST, QST
734 
735  real(RP) :: RAINT
736  real(RP) :: ROFFR, ROFFB, ROFFG ! runoff [kg/m2]
737 
738  real(RP) :: LUP, LDN, RUP
739  real(RP) :: SUP, SDN
740 
741  real(RP) :: LNET, SNET, FLXUV
742  real(RP) :: SW ! shortwave radition [W/m/m]
743  real(RP) :: LW ! longwave radition [W/m/m]
744  real(RP) :: psim,psim2,psim10 ! similality stability shear function for momentum
745  real(RP) :: psih,psih2,psih10 ! similality stability shear function for heat
746 
747  ! for shadow effect model
748  ! real(RP) :: HOUI1, HOUI2, HOUI3, HOUI4, HOUI5, HOUI6, HOUI7, HOUI8
749  ! real(RP) :: SLX, SLX1, SLX2, SLX3, SLX4, SLX5, SLX6, SLX7, SLX8
750  ! real(RP) :: THEATAZ ! Solar Zenith Angle [rad]
751  ! real(RP) :: THEATAS ! = PI/2. - THETAZ
752  ! real(RP) :: FAI ! Latitude [rad]
753 
754  real(RP) :: SR, SB, SG, RR, RB, RG
755  real(RP) :: SR1, SB1, SB2, SG1, SG2
756  real(RP) :: RB1, RB2, RG1, RG2
757  real(RP) :: HR, ELER, G0R
758  real(RP) :: HB, ELEB, G0B
759  real(RP) :: HG, ELEG, G0G
760 
761  real(RP) :: Z
762  real(RP) :: QS0R, QS0B, QS0G
763 
764  real(RP) :: RIBR, BHR, CDR
765  real(RP) :: RIBC, BHC, CDC
766  real(RP) :: ALPHAB, ALPHAG, ALPHAC
767  real(RP) :: CHR, CHB, CHG, CHC
768  real(RP) :: TC1, TC2, QC1, QC2
769 ! real(RP) :: CAPL1, AKSL1
770 
771  real(RP) :: DTS_MAX_onestep = 0.0_rp ! DTS_MAX * dt
772  real(RP) :: resi1,resi2 ! residual
773  real(RP) :: resi1p,resi2p ! residual
774  real(RP) :: G0RP,G0BP,G0GP
775 
776  real(RP) :: XXX, XXX2, XXX10
777  real(RP) :: THA,THC,THS,THS1,THS2
778  real(RP) :: RovCP
779  real(RP) :: EXN ! exner function at the surface
780  real(RP) :: qdry
781 
782  real(RP) :: tahdiurnal
783 
784  integer :: iteration
785 
786  !-----------------------------------------------------------
787  ! Set parameters
788  !-----------------------------------------------------------
789 
790  rovcp = rdry / cpdry
791  tha = ta * ( pre00 / prsa )**rovcp
792 
793  !--- Renew surface and layer temperatures
794 
795  trp = tr
796  tbp = tb
797  tgp = tg
798  tcp = tc
799  qcp = qc
800  !
801  trlp = trl
802  tblp = tbl
803  tglp = tgl
804  !
805  rainrp = rainr
806  rainbp = rainb
807  raingp = raing
808 
809  !--- Calculate AH data at LST
810  if ( tloc == 24 ) then
811  tahdiurnal = ( 1.0_rp-dsec ) * ahdiurnal(tloc ) &
812  + ( dsec ) * ahdiurnal(1 )
813  else
814  tahdiurnal = ( 1.0_rp-dsec ) * ahdiurnal(tloc ) &
815  + ( dsec ) * ahdiurnal(tloc+1)
816  endif
817 
818  ah_t = ah * tahdiurnal
819  alh_t = alh * tahdiurnal
820 
821  !--- limiter for surface temp change
822  dts_max_onestep = dts_max * dt
823 
824  if ( zdc + z0c + 2.0_rp >= za ) then
825  log_error("URBAN_DYN_kusaka01_SLC_main",*) 'ZDC + Z0C + 2m must be less than the 1st level! STOP.'
826  call prc_abort
827  ! "2.0m" has no special meaning, but it is related with BB formulation from Inoue (1963). Please see subroutine "canopy_wind".
828  ! The canopy model is modeled under an assumption that urban canopy lies below the lowest level of atmospheric model.
829  endif
830 
831  w = 2.0_rp * 1.0_rp * hgt
832  vfgs = svf
833  vfgw = 1.0_rp - svf
834  vfwg = ( 1.0_rp - svf ) * ( 1.0_rp - r ) / w
835  vfws = vfwg
836  vfww = 1.0_rp - 2.0_rp * vfwg
837 
838  sx = ssg(1) + ssg(2) ! downward shortwave radiation [W/m2]
839  rx = llg(1) + llg(2) ! downward longwave radiation [W/m2]
840 
841  ssgd = ssg(1) ! downward direct shortwave radiation [W/m2]
842  ssgq = ssg(2) ! downward diffuse shortwave radiation [W/m2]
843 
844  !--- calculate canopy wind
845 
846  call canopy_wind(za, ua, uc)
847 
848  !-----------------------------------------------------------
849  ! Set evaporation efficiency on roof/wall/road
850  !-----------------------------------------------------------
851 
852  !!--- calculate evaporation efficiency
853  raint = 1.0_rp * ( ( rain + snow ) * dt ) ! [kg/m2/s -> kg/m2]
854  call cal_beta(betr, raint, rainr, strgr, roffr)
855 
856  raint = 0.1_rp * ( ( rain + snow ) * dt )
857  call cal_beta(betb, raint, rainb, strgb, roffb)
858 
859  raint = 0.9_rp * ( ( rain + snow ) * dt )
860  call cal_beta(betg, raint, raing, strgg, roffg)
861 
862  roff = roff + r * roffr + rw * ( roffb + roffg )
863 
864  !-----------------------------------------------------------
865  ! Radiation : Net Short Wave Radiation at roof/wall/road
866  !-----------------------------------------------------------
867 
868  if( sx > 0.0_rp ) then ! SSG is downward short
869 
870  ! currently we use no shadow effect model
871  !! IF(.NOT.SHADOW) THEN ! no shadow effects model
872 
873  sr1 = sx * ( 1.0_rp - albr )
874  sg1 = sx * vfgs * ( 1.0_rp - albg )
875  sb1 = sx * vfws * ( 1.0_rp - albb )
876  sg2 = sb1 * albb / ( 1.0_rp - albb ) * vfgw * ( 1.0_rp - albg )
877  sb2 = sg1 * albg / ( 1.0_rp - albg ) * vfwg * ( 1.0_rp - albb )
878 
879  sr = sr1
880  sg = sg1 + sg2
881  sb = sb1 + sb2
882  snet = r * sr + w * sb + rw * sg
883 
884  else
885 
886  sr = 0.0_rp
887  sg = 0.0_rp
888  sb = 0.0_rp
889  snet = 0.0_rp
890 
891  end if
892 
893 
894  exn = ( prss / pre00 )**rovcp ! exner function
895 
896  qdry = 1.0_rp - qa
897 
898  !-----------------------------------------------------------
899  ! Energy balance on roof/wall/road surface
900  !-----------------------------------------------------------
901 
902  !--------------------------------------------------
903  ! Roof
904  !--------------------------------------------------
905 
906  ! new scheme
907 
908  g0rp = 0.0_rp
909  xxxr = 0.0_rp
910  resi1p = 0.0_rp
911  do iteration = 1, 100
912 
913  ths = tr / exn ! potential temp
914 
915  z = za - zdc
916  bhr = log(z0r/z0hr) / 0.4_rp
917  ribr = ( grav * 2.0_rp / (tha+ths) ) * (tha-ths) * (z+z0r) / (ua*ua)
918  call mos(xxxr,chr,cdr,bhr,ribr,z,z0r,ua,tha,ths,rhoo)
919 
920  call qsat( tr, prss, qdry, & ! [IN]
921  qs0r ) ! [OUT]
922 
923  rr = epsr * ( rx - stb * (tr**4) )
924  !HR = RHOO * CPdry * CHR * UA * (TR-TA)
925  hr = rhoo * cpdry * chr * ua * (ths-tha) * exn
926  eler = min( rhoo * chr * ua * betr * (qs0r-qa), real(RAINR/dt,RP) ) * lhv
927  g0r = sr + rr - hr - eler
928 
929  !--- calculate temperature in roof
930  ! if ( STRGR /= 0.0_RP ) then
931  ! CAPL1 = CAP_water * (RAINR / (DZR(1) + RAINR)) + CAPR * (DZR(1) / (DZR(1) + RAINR))
932  ! AKSL1 = AKS_water * (RAINR / (DZR(1) + RAINR)) + AKSR * (DZR(1) / (DZR(1) + RAINR))
933  ! else
934  ! CAPL1 = CAPR
935  ! AKSL1 = AKSR
936  ! endif
937  !! 1st layer's cap, aks are replaced.
938  !! call multi_layer2(UKE,BOUND,G0R,CAPR,AKSR,TRL,DZR,dt,TRLEND,CAPL1,AKSL1)
939 
940  trl = trlp
941  call multi_layer(uke,bound,g0r,capr,aksr,trl,dzr,dt,trlend)
942  resi1 = trl(1) - tr
943 
944  ! LOG_INFO("URBAN_DYN_kusaka01_SLC_main",'(a3,i5,f8.3,6f15.5)') "TR,",iteration,TR,G0R,SR,RR,HR,ELER,resi1
945 
946  if( abs(resi1) < sqrt(eps) ) then
947  tr = trl(1)
948  tr = max( trp - dts_max_onestep, min( trp + dts_max_onestep, tr ) )
949  exit
950  endif
951 
952  if ( resi1*resi1p < 0.0_rp ) then
953  tr = (tr + trl(1)) * 0.5_rp
954  else
955  tr = trl(1)
956  endif
957  tr = max( trp - dts_max_onestep, min( trp + dts_max_onestep, tr ) )
958 
959  resi1p = resi1
960 
961  enddo
962 
963 ! if( .NOT. (resi1 < sqrt(EPS)) ) then
964 ! LOG_WARN("URBAN_DYN_kusaka01_SLC_main",*) 'Warning not converged for TR in URBAN SLC', &
965 ! PRC_myrank, i,j, &
966 ! resi1, G0R
967 ! end if
968 
969  ! output for debug
970  if ( iteration > 100 ) then
971  log_warn("URBAN_DYN_kusaka01_SLC_main",*) 'iteration for TR was not converged',prc_myrank,i,j
972  log_info_cont(*) '---------------------------------------------------------------------------------'
973  log_info_cont(*) 'DEBUG Message --- Residual [K] :', resi1
974  log_newline
975  log_info_cont(*) 'DEBUG Message --- TRP : Initial TR [K] :', trp
976  log_info_cont(*) 'DEBUG Message --- TRLP: Initial TRL [K] :', trlp
977  log_newline
978  log_info_cont(*) 'DEBUG Message --- SX : Shortwave radiation [W/m2] :', sx
979  log_info_cont(*) 'DEBUG Message --- RX : Longwave radiation [W/m2] :', rx
980  log_info_cont(*) 'DEBUG Message --- PRSS: Surface pressure [Pa] :', prss
981  log_info_cont(*) 'DEBUG Message --- PRSA: Pressure at 1st atmos layer [m] :', prsa
982  log_info_cont(*) 'DEBUG Message --- RHOO: Air density [kg/m3] :', rhoo
983  log_info_cont(*) 'DEBUG Message --- ZA : Height at 1st atmos layer [m] :', za
984  log_info_cont(*) 'DEBUG Message --- TA : Temperature at 1st atmos layer [K] :', ta
985  log_info_cont(*) 'DEBUG Message --- UA : Wind speed at 1st atmos layer [m/s] :', ua
986  log_info_cont(*) 'DEBUG Message --- QA : Specific humidity at 1st atmos layer [kg/kg] :', qa
987  log_info_cont(*) 'DEBUG Message --- DZR : Depth of surface layer [m] :', dzr
988  log_newline
989  log_info_cont(*) 'DEBUG Message --- R, W, RW : Normalized height and road width [-] :', r, w,rw
990  log_info_cont(*) 'DEBUG Message --- SVF : Sky View Factors [-] :', svf
991  log_info_cont(*) 'DEBUG Message --- BETR: Evaporation efficiency [-] :', betr
992  log_info_cont(*) 'DEBUG Message --- EPSR: Surface emissivity of roof [-] :', epsr
993  log_info_cont(*) 'DEBUG Message --- CAPR: Heat capacity of roof [J m-3 K] :', capr
994  log_info_cont(*) 'DEBUG Message --- AKSR: Thermal conductivity of roof [W m-1 K] :', aksr
995  log_info_cont(*) 'DEBUG Message --- QS0R: Surface specific humidity [kg/kg] :', qs0r
996  log_info_cont(*) 'DEBUG Message --- ZDC : Desplacement height of canopy [m] :', zdc
997  log_info_cont(*) 'DEBUG Message --- Z0R : Momentum roughness length of roof [m] :', z0r
998  log_info_cont(*) 'DEBUG Message --- Z0HR: Thermal roughness length of roof [m] :', z0hr
999  log_info_cont(*) '---------------------------------------------------------------------------------'
1000  endif
1001 
1002  !--- update only fluxes ----
1003  ths = tr / exn
1004  ribr = ( grav * 2.0_rp / (tha+ths) ) * (tha-ths) * (z+z0r) / (ua*ua)
1005  call mos(xxxr,chr,cdr,bhr,ribr,z,z0r,ua,tha,ths,rhoo)
1006 
1007  call qsat( tr, prss, qdry, & ! [IN]
1008  qs0r ) ! [OUT]
1009 
1010  rr = epsr * ( rx - stb * (tr**4) )
1011  hr = rhoo * cpdry * chr * ua * (ths-tha) * exn
1012  eler = min( rhoo * chr * ua * betr * (qs0r-qa), real(RAINR/dt,RP) ) * lhv
1013  g0r = sr + rr - hr - eler
1014 
1015  trl = trlp
1016  call multi_layer(uke,bound,g0r,capr,aksr,trl,dzr,dt,trlend)
1017  resi1 = trl(1) - tr
1018  tr = trl(1)
1019 
1020  if ( abs(resi1) > dts_max_onestep ) then
1021  if ( abs(resi1) > dts_max_onestep*10.0_rp ) then
1022  log_error("URBAN_DYN_Kusaka01_main",*) 'tendency of TR exceeded a limit! STOP.'
1023  log_error_cont(*) 'previous TR and updated TR(TRL(1)) is ',tr-resi1, tr
1024  call prc_abort
1025  endif
1026  log_warn("URBAN_DYN_Kusaka01_main",*) 'tendency of TR exceeded a limit'
1027  log_info_cont(*) 'previous TR and updated TR(TRL(1)) is ', tr-resi1, tr
1028  endif
1029 
1030  !--------------------------------------------------
1031  ! Wall and Road
1032  !--------------------------------------------------
1033 
1034  ! new scheme
1035 
1036  ! empirical form
1037  alphab = 6.15_rp + 4.18_rp * uc
1038  if( uc > 5.0_rp ) alphab = 7.51_rp * (uc**0.78_rp )
1039  alphag = 6.15_rp + 4.18_rp * uc
1040  if( uc > 5.0_rp ) alphag = 7.51_rp * (uc**0.78_rp )
1041  chb = alphab / rhoo / cpdry / uc
1042  chg = alphag / rhoo / cpdry / uc
1043 
1044  g0bp = 0.0_rp
1045  g0gp = 0.0_rp
1046  xxxc = 0.0_rp
1047  resi1p = 0.0_rp
1048  resi2p = 0.0_rp
1049 
1050  do iteration = 1, 200
1051 
1052  ths1 = tb / exn
1053  ths2 = tg / exn
1054  thc = tc / exn
1055 
1056  z = za - zdc
1057  bhc = log(z0c/z0hc) / 0.4_rp
1058  ribc = ( grav * 2.0_rp / (tha+thc) ) * (tha-thc) * (z+z0c) / (ua*ua)
1059  call mos(xxxc,chc,cdc,bhc,ribc,z,z0c,ua,tha,thc,rhoo)
1060  alphac = chc * rhoo * cpdry * ua
1061 
1062  call qsat( tb, prss, qdry, qs0b )
1063  call qsat( tg, prss, qdry, qs0g )
1064 
1065  tc1 = rw*alphac + rw*alphag + w*alphab
1066  !TC2 = RW*ALPHAC*TA + RW*ALPHAG*TG + W*ALPHAB*TB
1067  tc2 = rw*alphac*tha + w*alphab*ths1 + rw*alphag*ths2
1068  thc = tc2 / tc1
1069  qc1 = rw*(chc*ua) + rw*(chg*betg*uc) + w*(chb*betb*uc)
1070  qc2 = rw*(chc*ua)*qa + rw*(chg*betg*uc)*qs0g + w*(chb*betb*uc)*qs0b
1071  qc = qc2 / qc1
1072 
1073  rg1 = epsg * ( rx * vfgs &
1074  + epsb * vfgw * stb * tb**4 &
1075  - stb * tg**4 )
1076 
1077  rb1 = epsb * ( rx * vfws &
1078  + epsg * vfwg * stb * tg**4 &
1079  + epsb * vfww * stb * tb**4 &
1080  - stb * tb**4 )
1081 
1082  rg2 = epsg * ( (1.0_rp-epsb) * vfgw * vfws * rx &
1083  + (1.0_rp-epsb) * vfgw * vfwg * epsg * stb * tg**4 &
1084  + epsb * (1.0_rp-epsb) * vfgw * vfww * stb * tb**4 )
1085 
1086  rb2 = epsb * ( (1.0_rp-epsg) * vfwg * vfgs * rx &
1087  + (1.0_rp-epsg) * epsb * vfgw * vfwg * stb * tb**4 &
1088  + (1.0_rp-epsb) * vfws * vfww * rx &
1089  + (1.0_rp-epsb) * vfwg * vfww * stb * epsg * tg**4 &
1090  + epsb * (1.0_rp-epsb) * vfww * (1.0_rp-2.0_rp*vfws) * stb * tb**4 )
1091 
1092  rg = rg1 + rg2
1093  rb = rb1 + rb2
1094 
1095  hb = rhoo * cpdry * chb * uc * (ths1-thc) * exn
1096  eleb = min( rhoo * chb * uc * betb * (qs0b-qc), real(RAINB/dt,RP) ) * lhv
1097  g0b = sb + rb - hb - eleb
1098 
1099  hg = rhoo * cpdry * chg * uc * (ths2-thc) * exn
1100  eleg = min( rhoo * chg * uc * betg * (qs0g-qc), real(RAING/dt,RP) ) * lhv
1101  g0g = sg + rg - hg - eleg
1102 
1103  tbl = tblp
1104  call multi_layer(uke,bound,g0b,capb,aksb,tbl,dzb,dt,tblend)
1105  resi1 = tbl(1) - tb
1106 
1107  tgl = tglp
1108  call multi_layer(uke,bound,g0g,capg,aksg,tgl,dzg,dt,tglend)
1109  resi2 = tgl(1) - tg
1110 
1111  !-----------
1112  !print *,HB, RHOO , CPdry , CHB , UC , THS1,THC
1113  !print *,HG, RHOO , CPdry , CHG , UC , THS2,THC
1114  !print *,ELEB ,RHOO , LHV , CHB , UC , BETB , QS0B , QC
1115  !print *,ELEG ,RHOO , LHV , CHG , UC , BETG , QS0G , QC
1116 
1117  !LOG_INFO("SLC_main",'(a3,i5,f8.3,6f15.5)') "TB,",iteration,TB,G0B,SB,RB,HB,ELEB,resi1
1118  !LOG_INFO("SLC_main",'(a3,i5,f8.3,6f15.5)') "TG,",iteration,TG,G0G,SG,RG,HG,ELEG,resi2
1119  !LOG_INFO("SLC_main",'(a3,i5,f8.3,3f15.5)') "TC,",iteration,TC,QC,QS0B,QS0G
1120  !--------
1121  !resi1 = abs(G0B - G0BP)
1122  !resi2 = abs(G0G - G0GP)
1123  !G0BP = G0B
1124  !G0GP = G0G
1125 
1126  if ( abs(resi1) < sqrt(eps) .AND. abs(resi2) < sqrt(eps) ) then
1127  tb = tbl(1)
1128  tg = tgl(1)
1129  tb = max( tbp - dts_max_onestep, min( tbp + dts_max_onestep, tb ) )
1130  tg = max( tgp - dts_max_onestep, min( tgp + dts_max_onestep, tg ) )
1131  exit
1132  endif
1133 
1134  if ( resi1*resi1p < 0.0_rp ) then
1135  tb = (tb + tbl(1)) * 0.5_rp
1136  else
1137  tb = tbl(1)
1138  endif
1139  if ( resi2*resi2p < 0.0_rp ) then
1140  tg = (tg + tgl(1)) * 0.5_rp
1141  else
1142  tg = tgl(1)
1143  endif
1144  tb = max( tbp - dts_max_onestep, min( tbp + dts_max_onestep, tb ) )
1145  tg = max( tgp - dts_max_onestep, min( tgp + dts_max_onestep, tg ) )
1146 
1147  resi1p = resi1
1148  resi2p = resi2
1149 
1150  ! this is for TC, QC
1151  call qsat( tb, prss, qdry, qs0b )
1152  call qsat( tg, prss, qdry, qs0g )
1153 
1154  ths1 = tb / exn
1155  ths2 = tg / exn
1156 
1157  tc1 = rw*alphac + rw*alphag + w*alphab
1158  !TC2 = RW*ALPHAC*THA + RW*ALPHAG*TG + W*ALPHAB*TB
1159  tc2 = rw*alphac*tha + w*alphab*ths1 + rw*alphag*ths2
1160  thc = tc2 / tc1
1161  tc = thc * exn
1162 
1163  qc1 = rw*(chc*ua) + rw*(chg*betg*uc) + w*(chb*betb*uc)
1164  qc2 = rw*(chc*ua)*qa + rw*(chg*betg*uc)*qs0g + w*(chb*betb*uc)*qs0b
1165  qc = qc2 / qc1
1166 
1167  enddo
1168 
1169 ! if( .NOT. (resi1 < sqrt(EPS) .AND. resi2 < sqrt(EPS) ) ) then
1170 ! LOG_INFO("SLC_main",*) 'Warning not converged for TG, TB in URBAN SLC', &
1171 ! PRC_myrank, i,j, &
1172 ! resi1, resi2, TB, TG, TC, G0BP, G0GP, RB, HB, RG, HG, QC, &
1173 ! CHC, CDC, BHC, ALPHAC
1174 ! LOG_INFO_CONT(*) TBP, TGP, TCP, &
1175 ! PRSS, THA, UA, QA, RHOO, UC, QCP, &
1176 ! ZA, ZDC, Z0C, Z0HC, &
1177 ! CHG, CHB, &
1178 ! W, RW, ALPHAG, ALPHAB, BETG, BETB, &
1179 ! RX, VFGS, VFGW, VFWG, VFWW, VFWS, STB, &
1180 ! SB, SG, LHV, TBLP, TGLP
1181 ! LOG_INFO_CONT(*) "6",VFGS, VFGW, VFWG, VFWW, VFWS
1182 ! end if
1183 
1184  ! output for debug
1185  if ( iteration > 200 ) then
1186  log_warn("URBAN_DYN_Kusaka01_main",*) 'iteration for TB/TG was not converged',prc_myrank,i,j
1187  log_info_cont(*) '---------------------------------------------------------------------------------'
1188  log_info_cont(*) 'DEBUG Message --- Residual [K] :', resi1,resi2
1189  log_newline
1190  log_info_cont(*) 'DEBUG Message --- TBP : Initial TB [K] :', tbp
1191  log_info_cont(*) 'DEBUG Message --- TBLP: Initial TBL [K] :', tblp
1192  log_info_cont(*) 'DEBUG Message --- TGP : Initial TG [K] :', tgp
1193  log_info_cont(*) 'DEBUG Message --- TGLP: Initial TGL [K] :', tglp
1194  log_info_cont(*) 'DEBUG Message --- TCP : Initial TC [K] :', tcp
1195  log_info_cont(*) 'DEBUG Message --- QCP : Initial QC [K] :', qcp
1196  log_newline
1197  log_info_cont(*) 'DEBUG Message --- UC : Canopy wind [m/s] :', uc
1198  log_info_cont(*) 'DEBUG Message --- SX : Shortwave radiation [W/m2] :', sx
1199  log_info_cont(*) 'DEBUG Message --- RX : Longwave radiation [W/m2] :', rx
1200  log_info_cont(*) 'DEBUG Message --- PRSS: Surface pressure [Pa] :', prss
1201  log_info_cont(*) 'DEBUG Message --- PRSA: Pressure at 1st atmos layer [m] :', prsa
1202  log_info_cont(*) 'DEBUG Message --- RHOO: Air density [kg/m3] :', rhoo
1203  log_info_cont(*) 'DEBUG Message --- ZA : Height at 1st atmos layer [m] :', za
1204  log_info_cont(*) 'DEBUG Message --- TA : Temperature at 1st atmos layer [K] :', ta
1205  log_info_cont(*) 'DEBUG Message --- UA : Wind speed at 1st atmos layer [m/s] :', ua
1206  log_info_cont(*) 'DEBUG Message --- QA : Specific humidity at 1st atmos layer [kg/kg] :', qa
1207  log_info_cont(*) 'DEBUG Message --- DZB : Depth of surface layer [m] :', dzb
1208  log_info_cont(*) 'DEBUG Message --- DZG : Depth of surface layer [m] :', dzg
1209  log_newline
1210  log_info_cont(*) 'DEBUG Message --- R, W, RW : Normalized height and road width [-] :', r, w,rw
1211  log_info_cont(*) 'DEBUG Message --- SVF : Sky View Factors [-] :', svf
1212  log_info_cont(*) 'DEBUG Message --- BETB,BETG : Evaporation efficiency [-] :', betb,betg
1213  log_info_cont(*) 'DEBUG Message --- EPSB,EPSG : Surface emissivity [-] :', epsb,epsg
1214  log_info_cont(*) 'DEBUG Message --- CAPB,CAPG : Heat capacity [J m-3 K] :', capb,capg
1215  log_info_cont(*) 'DEBUG Message --- AKSB,AKSG : Thermal conductivity [W m-1 K] :', aksb,aksb
1216  log_info_cont(*) 'DEBUG Message --- QS0B,QS0G : Surface specific humidity [kg/kg] :', qs0b,qs0g
1217  log_info_cont(*) 'DEBUG Message --- ZDC : Desplacement height of canopy [m] :', zdc
1218  log_info_cont(*) 'DEBUG Message --- Z0C : Momentum roughness length of canopy [m] :', z0c
1219  log_info_cont(*) 'DEBUG Message --- Z0HC : Thermal roughness length of canopy [m] :', z0hc
1220  log_info_cont(*) '---------------------------------------------------------------------------------'
1221  endif
1222 
1223 
1224  !--- update only fluxes ----
1225  rg1 = epsg * ( rx * vfgs &
1226  + epsb * vfgw * stb * tb**4 &
1227  - stb * tg**4 )
1228  rb1 = epsb * ( rx * vfws &
1229  + epsg * vfwg * stb * tg**4 &
1230  + epsb * vfww * stb * tb**4 &
1231  - stb * tb**4 )
1232 
1233  rg2 = epsg * ( (1.0_rp-epsb) * vfgw * vfws * rx &
1234  + (1.0_rp-epsb) * vfgw * vfwg * epsg * stb * tg**4 &
1235  + epsb * (1.0_rp-epsb) * vfgw * vfww * stb * tb**4 )
1236  rb2 = epsb * ( (1.0_rp-epsg) * vfwg * vfgs * rx &
1237  + (1.0_rp-epsg) * epsb * vfgw * vfwg * stb * tb**4 &
1238  + (1.0_rp-epsb) * vfws * vfww * rx &
1239  + (1.0_rp-epsb) * vfwg * vfww * stb * epsg * tg**4 &
1240  + epsb * (1.0_rp-epsb) * vfww * (1.0_rp-2.0_rp*vfws) * stb * tb**4 )
1241 
1242  rg = rg1 + rg2
1243  rb = rb1 + rb2
1244 
1245  ths1 = tb / exn
1246  ths2 = tg / exn
1247  thc = tc / exn
1248 
1249  hb = rhoo * cpdry * chb * uc * (ths1-thc) * exn
1250  eleb = min( rhoo * chb * uc * betb * (qs0b-qc), real(RAINB/dt,RP) ) * lhv
1251  g0b = sb + rb - hb - eleb
1252 
1253  hg = rhoo * cpdry * chg * uc * (ths2-thc) * exn
1254  eleg = min( rhoo * chg * uc * betg * (qs0g-qc), real(RAING/dt,RP) ) * lhv
1255  g0g = sg + rg - hg - eleg
1256 
1257  tbl = tblp
1258  call multi_layer(uke,bound,g0b,capb,aksb,tbl,dzb,dt,tblend)
1259  resi1 = tbl(1) - tb
1260  tb = tbl(1)
1261 
1262  tgl = tglp
1263  call multi_layer(uke,bound,g0g,capg,aksg,tgl,dzg,dt,tglend)
1264  resi2 = tgl(1) - tg
1265  tg = tgl(1)
1266 
1267  if ( abs(resi1) > dts_max_onestep ) then
1268  if ( abs(resi1) > dts_max_onestep*10.0_rp ) then
1269  log_error("URBAN_DYN_Kusaka01_main",*) 'tendency of TB exceeded a limit! STOP.'
1270  log_error_cont(*) 'previous TB and updated TB(TBL(1)) is ', tb-resi1,tb
1271  call prc_abort
1272  endif
1273  log_warn("URBAN_DYN_Kusaka01_main",*) 'tendency of TB exceeded a limit'
1274  log_info_cont(*) 'previous TB and updated TB(TBL(1)) is ', tb-resi1, tb
1275  endif
1276 
1277  if ( abs(resi2) > dts_max_onestep ) then
1278  if ( abs(resi2) > dts_max_onestep*10.0_rp ) then
1279  log_error("URBAN_DYN_Kusaka01_main",*) 'tendency of TG exceeded a limit! STOP.'
1280  log_error_cont(*) 'previous TG and updated TG(TGL(1)) is ', tg-resi2, tg, resi2
1281  call prc_abort
1282  endif
1283  log_warn("URBAN_DYN_Kusaka01_main",*) 'tendency of TG exceeded a limit'
1284  log_info_cont(*) 'previous TG and updated TG(TGL(1)) is ', tg-resi2, tg
1285  endif
1286 
1287  !-----------------------------------------------------------
1288  ! Total Fluxes from Urban Canopy
1289  !-----------------------------------------------------------
1290 
1291  flxuv = ( r*cdr + rw*cdc ) * ua * ua
1292  sh = ( r*hr + w*hb + rw*hg ) ! Sensible heat flux [W/m/m]
1293  lh = ( r*eler + w*eleb + rw*eleg ) ! Latent heat flux [W/m/m]
1294  ghflx = r*g0r + w*g0b + rw*g0g
1295  lnet = r*rr + w*rb + rw*rg
1296 
1297  !-----------------------------------------------------------
1298  ! Grid average
1299  !-----------------------------------------------------------
1300 
1301  lw = rx - lnet ! Upward longwave radiation [W/m/m]
1302  sw = sx - snet ! Upward shortwave radiation [W/m/m]
1303  rn = (snet+lnet) ! Net radiation [W/m/m]
1304 
1305  !--- shortwave radiation
1306  sdn = r + w * (vfws + vfgs * albg * vfwg) + rw * (vfgs + vfws * albb * vfgw)
1307  sup = r * albr &
1308  + w * ( vfws * albb + vfgs * albg * vfwg *albb ) &
1309  + rw * ( vfgs * albg + vfws * albb * vfgw *albg )
1310 
1311  albd_sw_grid = sup / sdn
1312 
1313  !--- longwave radiation
1314  ldn = r + w*vfws + rw*vfgs
1315  lup = r * (1.0_rp-epsr) &
1316  + w*( (1.0_rp-epsb*vfww)*(1.0_rp-epsb)*vfws - epsb*vfwg*(1.0_rp-epsg)*vfgs ) &
1317  + rw*( (1.0_rp-epsg)*vfgs - epsg*(1.0_rp-vfgs)*(1.0_rp-epsb)*vfws )
1318 
1319  rup = (ldn - lup) * rx - lnet
1320  albd_lw_grid = lup / ldn
1321 
1322 
1323  ! RUP = R * (EPSR * STB * TR**4 ) &
1324  ! + W * (EPSB*STB*(TB**4) - EPSB*EPSG*VFWG*STB*(TG**4) - EPSB*EPSB*VFWW*STB*(TB**4) &
1325  ! - EPSB *(1.0_RP-EPSG) * EPSB * VFGW * VFWG * STB * (TB**4) &
1326  ! - EPSB *(1.0_RP-EPSB) * VFWG * VFWW * STB * EPSG * (TG**4) &
1327  ! - EPSB * EPSB * (1.0_RP-EPSB) * VFWW * VFWW * STB * (TB**4) ) &
1328  ! + RW * (EPSG*STB*(TG**4) - EPSG * EPSB * VFGW * STB * (TB**4) &
1329  ! - EPSG * EPSB * (1.0_RP-EPSB) * (1.0_RP-SVF) * VFWW * STB * TB**4 )
1330 
1331 
1332  shr = hr ! Sensible heat flux on roof [W/m/m]
1333  shb = hb ! Sensible heat flux on wall [W/m/m]
1334  shg = hg ! Sensible heat flux on road [W/m/m]
1335  lhr = eler ! Latent heat flux on road [W/m/m]
1336  lhb = eleb ! Latent heat flux on wall [W/m/m]
1337  lhg = eleg ! Latent heat flux on road [W/m/m]
1338  ghr = g0r ! Ground heat flux on roof [W/m/m]
1339  ghb = g0b ! Ground heat flux on wall [W/m/m]
1340  ghg = g0g ! Ground heat flux on road [W/m/m]
1341  rnr = sr + rr ! Net radiation on roof [W/m/m]
1342  rnb = sb + rb ! Net radiation on building [W/m/m]
1343  rng = sg + rg ! Net radiation on ground [W/m/m]
1344 
1345  ! calculate rain amount remaining on the surface
1346  rainr = max(0.0_rp, rainr-(lhr/lhv)*real(dt,kind=rp)) ! [kg/m/m = mm]
1347  rainb = max(0.0_rp, rainb-(lhb/lhv)*real(dt,kind=rp)) ! [kg/m/m = mm]
1348  raing = max(0.0_rp, raing-(lhg/lhv)*real(dt,kind=rp)) ! [kg/m/m = mm]
1349 
1350  !-----------------------------------------------------------
1351  ! diagnostic GRID AVERAGED TS from upward logwave
1352  !-----------------------------------------------------------
1353 
1354  rts = ( rup / stb / ( 1.0_rp-albd_lw_grid) )**0.25
1355 
1356  !-----------------------------------------------------------
1357  ! diagnostic grid average U10, V10, T2, Q2 from urban
1358  ! Below method would be better to be improved. This is tentative method.
1359  !-----------------------------------------------------------
1360  !UST = sqrt( FLXUV ) ! u* [m/s]
1361  !TST = -SH / RHOO / CPdry / UST ! T* [K]
1362  !QST = -LH / RHOO / LHV / UST ! q* [-]
1363  !Z = ZA - ZDC
1364  !XXX = 0.4*9.81*Z*TST/TA/UST/UST
1365 
1366  xxx = xxxc
1367  call cal_psi(xxx,psim,psih)
1368 
1369  xxx2 = (2.0_rp/z) * xxx
1370  call cal_psi(xxx2,psim2,psih2)
1371 
1372  xxx10 = (10.0_rp/z) * xxx
1373  call cal_psi(xxx10,psim10,psih10)
1374 
1375  !U10 = U1 * ((log(10.0_RP/Z0C)-psim10)/(log(Z/Z0C)-psim)) ! u at 10 m [m/s]
1376  !V10 = V1 * ((log(10.0_RP/Z0C)-psim10)/(log(Z/Z0C)-psim)) ! v at 10 m [m/s]
1377  u10 = u1 * log(10.0_rp/z0c) / log(z/z0c)
1378  v10 = v1 * log(10.0_rp/z0c) / log(z/z0c)
1379 
1380  t2 = rts + (ta-rts)*((log(2.0_rp/z0hc)-psih2)/(log(z/z0hc)-psih))
1381  q2 = qc
1382 
1383  !-----------------------------------------------------------
1384  ! add anthropogenic heat fluxes
1385  !-----------------------------------------------------------
1386 
1387  sh = sh + ah_t ! Sensible heat flux [W/m/m]
1388  lh = lh + alh_t ! Latent heat flux [W/m/m]
1389 
1390  return
1391  end subroutine slc_main
1392 
1393  !-----------------------------------------------------------------------------
1394  subroutine canopy_wind(ZA, UA, UC)
1395  implicit none
1396 
1397  real(RP), intent(in) :: za ! height at 1st atmospheric level [m]
1398  real(RP), intent(in) :: UA ! wind speed at 1st atmospheric level [m/s]
1399  real(RP), intent(out) :: UC ! wind speed at 1st atmospheric level [m/s]
1400 
1401  real(RP) :: UR,ZC,XLB,BB
1402 
1403  if( zr + 2.0_rp < za ) then
1404  ur = ua * log((zr-zdc)/z0c) / log((za-zdc)/z0c)
1405  zc = 0.7_rp * zr
1406  xlb = 0.4_rp * (zr-zdc)
1407  ! BB formulation from Inoue (1963)
1408  bb = 0.4_rp * zr / ( xlb * log((zr-zdc)/z0c) )
1409  uc = ur * exp( -bb * (1.0_rp-zc/zr) )
1410  else
1411  ! PRINT *, 'Warning ZR + 2m is larger than the 1st WRF level'
1412  zc = za / 2.0_rp
1413  uc = ua / 2.0_rp
1414  endif
1415 
1416  uc = max(uc,0.01_rp)
1417 
1418  return
1419  end subroutine canopy_wind
1420 
1421  !-----------------------------------------------------------------------------
1422  subroutine cal_beta(BET, RAIN, WATER, STRG, ROFF)
1423  implicit none
1424 
1425  real(RP), intent(out) :: BET ! evapolation efficiency [-]
1426  real(RP), intent(in) :: RAIN ! precipitation [mm*dt]
1427  real(RP), intent(inout) :: WATER ! rain amount in strage [kg/m2]
1428  real(RP), intent(inout) :: STRG ! rain strage [kg/m2]
1429  real(RP), intent(inout) :: ROFF ! runoff [kg/m2]
1430 
1431  if ( strg == 0.0_rp ) then ! not consider evapolation from urban
1432  bet = 0.0_rp
1433  roff = rain
1434  else
1435  water = water + rain
1436  roff = max(0.0_rp, water-strg)
1437  water = water - max(0.0_rp, water-strg)
1438  bet = min( water / strg, 1.0_rp)
1439  endif
1440 
1441  return
1442  end subroutine cal_beta
1443 
1444  !-----------------------------------------------------------------------------
1445  subroutine cal_psi(zeta,psim,psih)
1446  use scale_const, only: &
1447  pi => const_pi
1448  implicit none
1449 
1450  real(RP), intent(inout) :: zeta ! z/L
1451  real(RP), intent(out) :: psim
1452  real(RP), intent(out) :: psih
1453  real(RP) :: X
1454 
1455  if( zeta >= 1.0_rp ) zeta = 1.0_rp
1456  if( zeta <= -5.0_rp ) zeta = -5.0_rp
1457 
1458  if( zeta > 0.0_rp ) then
1459  psim = -5.0_rp * zeta
1460  psih = -5.0_rp * zeta
1461  else
1462  x = ( 1.0_rp - 16.0_rp * zeta )**0.25_rp
1463  psim = 2.0_rp * log((1.0_rp+x)/2.0_rp) + log((1.0_rp+x*x)/2.0_rp) - 2.0_rp*atan(x) + pi/2.0_rp
1464  psih = 2.0_rp * log((1.0_rp+x*x)/2.0_rp)
1465  end if
1466 
1467  return
1468  end subroutine cal_psi
1469 
1470  !-----------------------------------------------------------------------------
1471  ! XXX: z/L (requires iteration by Newton-Rapson method)
1472  ! B1: Stanton number
1473  ! PSIM: = PSIX of LSM
1474  ! PSIH: = PSIT of LSM
1475  subroutine mos(XXX,CH,CD,B1,RIB,Z,Z0,UA,TA,TSF,RHO)
1476  use scale_const, only: &
1477  cpdry => const_cpdry ! CPP : heat capacity of dry air [J/K/kg]
1478  implicit none
1479 
1480  real(RP), intent(in) :: B1, Z, Z0, UA, TA, TSF, RHO
1481  real(RP), intent(out) :: CD, CH
1482  real(RP), intent(inout) :: XXX, RIB
1483  real(RP) :: XXX0, X, X0, FAIH, DPSIM, DPSIH
1484  real(RP) :: F, DF, XXXP, US, TS, AL, XKB, DD, PSIM, PSIH
1485  integer :: NEWT
1486  integer, parameter :: NEWT_END = 10
1487 
1488  real(RP) :: lnZ
1489 
1490  lnz = log( (z+z0)/z0 )
1491 
1492  if( rib <= -15.0_rp ) rib = -15.0_rp
1493 
1494  if( rib < 0.0_rp ) then
1495 
1496  do newt = 1, newt_end
1497 
1498  if( xxx >= 0.0_rp ) xxx = -1.0e-3_rp
1499 
1500  xxx0 = xxx * z0/(z+z0)
1501 
1502  x = (1.0_rp-16.0_rp*xxx)**0.25
1503  x0 = (1.0_rp-16.0_rp*xxx0)**0.25
1504 
1505  psim = lnz &
1506  - log( (x+1.0_rp)**2 * (x**2+1.0_rp) ) &
1507  + 2.0_rp * atan(x) &
1508  + log( (x+1.0_rp)**2 * (x0**2+1.0_rp) ) &
1509  - 2.0_rp * atan(x0)
1510  faih = 1.0_rp / sqrt( 1.0_rp - 16.0_rp*xxx )
1511  psih = lnz + 0.4_rp*b1 &
1512  - 2.0_rp * log( sqrt( 1.0_rp - 16.0_rp*xxx ) + 1.0_rp ) &
1513  + 2.0_rp * log( sqrt( 1.0_rp - 16.0_rp*xxx0 ) + 1.0_rp )
1514 
1515  dpsim = ( 1.0_rp - 16.0_rp*xxx )**(-0.25) / xxx &
1516  - ( 1.0_rp - 16.0_rp*xxx0 )**(-0.25) / xxx
1517  dpsih = 1.0_rp / sqrt( 1.0_rp - 16.0_rp*xxx ) / xxx &
1518  - 1.0_rp / sqrt( 1.0_rp - 16.0_rp*xxx0 ) / xxx
1519 
1520  f = rib * psim**2 / psih - xxx
1521 
1522  df = rib * ( 2.0_rp*dpsim*psim*psih - dpsih*psim**2 ) &
1523  / psih**2 - 1.0_rp
1524 
1525  xxxp = xxx
1526  xxx = xxxp - f / df
1527 
1528  if( xxx <= -10.0_rp ) xxx = -10.0_rp
1529 
1530  end do
1531 
1532  else if( rib >= 0.142857_rp ) then
1533 
1534  xxx = 0.714_rp
1535  psim = lnz + 7.0_rp * xxx
1536  psih = psim + 0.4_rp * b1
1537 
1538  else
1539 
1540  al = lnz
1541  xkb = 0.4_rp * b1
1542  dd = -4.0_rp * rib * 7.0_rp * xkb * al + (al+xkb)**2
1543  if( dd <= 0.0_rp ) dd = 0.0_rp
1544 
1545  xxx = ( al + xkb - 2.0_rp*rib*7.0_rp*al - sqrt(dd) ) / ( 2.0_rp * ( rib*7.0_rp**2 - 7.0_rp ) )
1546  psim = lnz + 7.0_rp * min( xxx, 0.714_rp )
1547  psih = psim + 0.4_rp * b1
1548 
1549  endif
1550 
1551  us = 0.4_rp * ua / psim ! u*
1552  if( us <= 0.01_rp ) us = 0.01_rp
1553  ts = 0.4_rp * (ta-tsf) / psih ! T*
1554 
1555  cd = us * us / ua**2 ! CD
1556  ch = 0.4_rp * us / psih / ua ! CH
1557  !ALPHA = RHO * CPdry * 0.4_RP * US / PSIH ! RHO*CP*CH*U
1558 
1559  return
1560  end subroutine mos
1561 
1562  !-------------------------------------------------------------------
1563  subroutine multi_layer(KM,BOUND,G0,CAP,AKS,TSL,DZ,DELT,TSLEND)
1564  !
1565  ! calculate temperature in roof/building/road
1566  ! multi-layer heat equation model
1567  ! Solving Heat Equation by Tri Diagonal Matrix Algorithm
1568  !-------------------------------------------------------------------
1569 
1570  implicit none
1571 
1572  real(RP), intent(in) :: G0
1573  real(RP), intent(in) :: CAP
1574  real(RP), intent(in) :: AKS
1575  real(DP), intent(in) :: DELT ! Tim setep [ s ]
1576  real(RP), intent(in) :: TSLEND
1577  integer, intent(in) :: KM
1578  integer, intent(in) :: BOUND
1579  real(RP), intent(in) :: DZ(km)
1580  real(RP), intent(inout) :: TSL(km)
1581  real(RP) :: A(km), B(km), C(km), D(km), P(km), Q(km)
1582  real(RP) :: DZEND
1583  integer :: K
1584 
1585  dzend = dz(km)
1586 
1587  a(1) = 0.0_rp
1588 
1589  b(1) = cap * dz(1) / delt &
1590  + 2.0_rp * aks / (dz(1)+dz(2))
1591  c(1) = -2.0_rp * aks / (dz(1)+dz(2))
1592  d(1) = cap * dz(1) / delt * tsl(1) + g0
1593 
1594  do k = 2, km-1
1595  a(k) = -2.0_rp * aks / (dz(k-1)+dz(k))
1596  b(k) = cap * dz(k) / delt + 2.0_rp * aks / (dz(k-1)+dz(k)) + 2.0_rp * aks / (dz(k)+dz(k+1))
1597  c(k) = -2.0_rp * aks / (dz(k)+dz(k+1))
1598  d(k) = cap * dz(k) / delt * tsl(k)
1599  end do
1600 
1601  if( bound == 1 ) then ! Flux=0
1602  a(km) = -2.0_rp * aks / (dz(km-1)+dz(km))
1603  b(km) = cap * dz(km) / delt + 2.0_rp * aks / (dz(km-1)+dz(km))
1604  c(km) = 0.0_rp
1605  d(km) = cap * dz(km) / delt * tsl(km)
1606  else ! T=constant
1607  a(km) = -2.0_rp * aks / (dz(km-1)+dz(km))
1608  b(km) = cap * dz(km) / delt + 2.0_rp * aks / (dz(km-1)+dz(km)) + 2.0_rp * aks / (dz(km)+dzend)
1609  c(km) = 0.0_rp
1610  d(km) = cap * dz(km) / delt * tsl(km) + 2.0_rp * aks * tslend / (dz(km)+dzend)
1611  end if
1612 
1613  p(1) = -c(1) / b(1)
1614  q(1) = d(1) / b(1)
1615 
1616  do k = 2, km
1617  p(k) = -c(k) / ( a(k) * p(k-1) + b(k) )
1618  q(k) = ( -a(k) * q(k-1) + d(k) ) / ( a(k) * p(k-1) + b(k) )
1619  end do
1620 
1621  tsl(km) = q(km)
1622 
1623  do k = km-1, 1, -1
1624  tsl(k) = p(k) * tsl(k+1) + q(k)
1625  end do
1626 
1627  return
1628  end subroutine multi_layer
1629 
1630  !-------------------------------------------------------------------
1631  subroutine multi_layer2(KM,BOUND,G0,CAP,AKS,TSL,DZ,DELT,TSLEND,CAP1,AKS1)
1632  !
1633  ! calculate temperature in roof/building/road
1634  ! multi-layer heat equation model
1635  ! Solving Heat Equation by Tri Diagonal Matrix Algorithm
1636  !-------------------------------------------------------------------
1637 
1638  implicit none
1639 
1640  real(RP), intent(in) :: G0
1641  real(RP), intent(in) :: CAP
1642  real(RP), intent(in) :: AKS
1643  real(RP), intent(in) :: CAP1 ! for 1st layer
1644  real(RP), intent(in) :: AKS1 ! for 1st layer
1645  real(DP), intent(in) :: DELT ! Time step [ s ]
1646  real(RP), intent(in) :: TSLEND
1647  integer, intent(in) :: KM
1648  integer, intent(in) :: BOUND
1649  real(RP), intent(in) :: DZ(km)
1650  real(RP), intent(inout) :: TSL(km)
1651  real(RP) :: A(km), B(km), C(km), D(km), X(km), P(km), Q(km)
1652  real(RP) :: DZEND
1653  integer :: K
1654 
1655  dzend = dz(km)
1656 
1657  a(1) = 0.0_rp
1658 
1659  b(1) = cap1 * dz(1) / delt &
1660  + 2.0_rp * aks1 / (dz(1)+dz(2))
1661  c(1) = -2.0_rp * aks1 / (dz(1)+dz(2))
1662  d(1) = cap1 * dz(1) / delt * tsl(1) + g0
1663 
1664  do k = 2, km-1
1665  a(k) = -2.0_rp * aks / (dz(k-1)+dz(k))
1666  b(k) = cap * dz(k) / delt + 2.0_rp * aks / (dz(k-1)+dz(k)) + 2.0_rp * aks / (dz(k)+dz(k+1))
1667  c(k) = -2.0_rp * aks / (dz(k)+dz(k+1))
1668  d(k) = cap * dz(k) / delt * tsl(k)
1669  end do
1670 
1671  if( bound == 1 ) then ! Flux=0
1672  a(km) = -2.0_rp * aks / (dz(km-1)+dz(km))
1673  b(km) = cap * dz(km) / delt + 2.0_rp * aks / (dz(km-1)+dz(km))
1674  c(km) = 0.0_rp
1675  d(km) = cap * dz(km) / delt * tsl(km)
1676  else ! T=constant
1677  a(km) = -2.0_rp * aks / (dz(km-1)+dz(km))
1678  b(km) = cap * dz(km) / delt + 2.0_rp * aks / (dz(km-1)+dz(km)) + 2.0_rp * aks / (dz(km)+dzend)
1679  c(km) = 0.0_rp
1680  d(km) = cap * dz(km) / delt * tsl(km) + 2.0_rp * aks * tslend / (dz(km)+dzend)
1681  end if
1682 
1683  p(1) = -c(1) / b(1)
1684  q(1) = d(1) / b(1)
1685 
1686  do k = 2, km
1687  p(k) = -c(k) / ( a(k) * p(k-1) + b(k) )
1688  q(k) = ( -a(k) * q(k-1) + d(k) ) / ( a(k) * p(k-1) + b(k) )
1689  end do
1690 
1691  x(km) = q(km)
1692 
1693  do k = km-1, 1, -1
1694  x(k) = p(k) * x(k+1) + q(k)
1695  end do
1696 
1697  do k = 1, km
1698  tsl(k) = x(k)
1699  enddo
1700 
1701  return
1702  end subroutine multi_layer2
1703 
1704  !-----------------------------------------------------------------------------
1705  subroutine urban_param_setup
1706  implicit none
1707 
1708  real(RP) :: DHGT,THGT,VFWS,VFGS
1709  integer :: k
1710 
1711  ! initialize
1712  r = 0.0_rp
1713  rw = 0.0_rp
1714  hgt = 0.0_rp
1715  z0hr = 0.0_rp
1716  z0hb = 0.0_rp
1717  z0hg = 0.0_rp
1718  z0c = 0.0_rp
1719  z0hc = 0.0_rp
1720  zdc = 0.0_rp
1721  svf = 0.0_rp
1722 
1723  ! set up other urban parameters
1724  z0hr = 0.1_rp * z0r
1725  z0hb = 0.1_rp * z0b
1726  z0hg = 0.1_rp * z0g
1727  zdc = zr * 0.3_rp
1728  z0c = zr * 0.15_rp
1729  z0hc = 0.1_rp * z0c
1730 
1731  ! HGT: Normalized height
1732  hgt = zr / ( road_width + roof_width )
1733 
1734  ! R: Normalized Roof Width (a.k.a. "building coverage ratio")
1735  r = roof_width / ( road_width + roof_width )
1736  rw = 1.0_rp - r
1737 
1738  ! Calculate Sky View Factor:
1739  dhgt = hgt / 100.0_rp
1740  thgt = 0.0_rp
1741  vfws = 0.0_rp
1742  thgt = hgt - dhgt / 2.0_rp
1743  do k = 1, 99
1744  thgt = thgt - dhgt
1745  vfws = vfws + 0.25_rp * ( 1.0_rp - thgt / sqrt( thgt**2 + rw**2 ) )
1746  end do
1747 
1748  vfws = vfws / 99.0_rp
1749  vfws = vfws * 2.0_rp
1750  vfgs = 1.0_rp - 2.0_rp * vfws * hgt / rw
1751  svf = vfgs
1752 
1753  return
1754  end subroutine urban_param_setup
1755 
1756  subroutine put_history( &
1757  UIA, UJA, &
1758  SHR, SHB, SHG, &
1759  LHR, LHB, LHG, &
1760  GHR, GHB, GHG, &
1761  RNR, RNB, RNG, &
1762  RNgrd )
1763  use scale_file_history, only: &
1764  file_history_put
1765  integer, intent(in) :: UIA, UJA
1766  real(RP), intent(in) :: SHR(uia,uja), SHB(uia,uja), SHG(uia,uja)
1767  real(RP), intent(in) :: LHR(uia,uja), LHB(uia,uja), LHG(uia,uja)
1768  real(RP), intent(in) :: GHR(uia,uja), GHB(uia,uja), GHG(uia,uja)
1769  real(RP), intent(in) :: RNR(uia,uja), RNB(uia,uja), RNG(uia,uja)
1770  real(RP), intent(in) :: RNgrd(uia,uja)
1771 
1772  call file_history_put( i_shr, shr(:,:) )
1773  call file_history_put( i_shb, shb(:,:) )
1774  call file_history_put( i_shg, shg(:,:) )
1775  call file_history_put( i_lhr, lhr(:,:) )
1776  call file_history_put( i_lhb, lhb(:,:) )
1777  call file_history_put( i_lhg, lhg(:,:) )
1778  call file_history_put( i_ghr, ghr(:,:) )
1779  call file_history_put( i_ghb, ghb(:,:) )
1780  call file_history_put( i_ghg, ghg(:,:) )
1781  call file_history_put( i_rnr, rnr(:,:) )
1782  call file_history_put( i_rnb, rnb(:,:) )
1783  call file_history_put( i_rng, rng(:,:) )
1784  call file_history_put( i_rngrd, rngrd(:,:) )
1785 
1786  return
1787  end subroutine put_history
1788 
1789 end module scale_urban_dyn_kusaka01
scale_const::const_cpdry
real(rp), public const_cpdry
specific heat (dry air,constant pressure) [J/kg/K]
Definition: scale_const.F90:56
scale_atmos_saturation
module atmosphere / saturation
Definition: scale_atmos_saturation.F90:12
scale_cpl_sfc_index
module coupler / surface-atmospehre
Definition: scale_cpl_sfc_index.F90:11
scale_cpl_sfc_index::i_r_vis
integer, parameter, public i_r_vis
Definition: scale_cpl_sfc_index.F90:31
scale_io::io_fid_conf
integer, public io_fid_conf
Config file ID.
Definition: scale_io.F90:55
scale_file_history::file_history_reg
subroutine, public file_history_reg(name, desc, unit, itemid, standard_name, ndims, dim_type, cell_measures, fill_halo)
Register/Append variable to history file.
Definition: scale_file_history.F90:649
scale_cpl_sfc_index::n_rad_dir
integer, parameter, public n_rad_dir
Definition: scale_cpl_sfc_index.F90:36
scale_cpl_sfc_index::n_rad_rgn
integer, parameter, public n_rad_rgn
Definition: scale_cpl_sfc_index.F90:28
scale_const::const_rdry
real(rp), public const_rdry
specific gas constant (dry air) [J/kg/K]
Definition: scale_const.F90:55
scale_atmos_hydrometeor
module atmosphere / hydrometeor
Definition: scale_atmos_hydrometeor.F90:12
scale_prc
module PROCESS
Definition: scale_prc.F90:11
scale_bulkflux::bulkflux
procedure(bc), pointer, public bulkflux
Definition: scale_bulkflux.F90:80
scale_prc::prc_myrank
integer, public prc_myrank
process num in local communicator
Definition: scale_prc.F90:89
scale_const::const_rvap
real(rp), parameter, public const_rvap
specific gas constant (water vapor) [J/kg/K]
Definition: scale_const.F90:63
scale_prc::prc_abort
subroutine, public prc_abort
Abort Process.
Definition: scale_prc.F90:338
scale_const
module CONSTANT
Definition: scale_const.F90:11
scale_cpl_sfc_index::i_r_direct
integer, parameter, public i_r_direct
Definition: scale_cpl_sfc_index.F90:37
scale_urban_dyn_kusaka01
module urban / dynamics / Kusaka01
Definition: scale_urban_dyn_kusaka01.F90:12
scale_urban_dyn_kusaka01::urban_dyn_kusaka01
subroutine, public urban_dyn_kusaka01(UKA, UKS, UKE, UIA, UIS, UIE, UJA, UJS, UJE, TMPA, PRSA, W1, U1, V1, DENS, QA, LHV, Z1, PBL, RHOS, PRSS, LWD, SWD, RAIN, SNOW, CDZ, fact_urban, tloc, dsec, dt, TRL_URB, TBL_URB, TGL_URB, TR_URB, TB_URB, TG_URB, TC_URB, QC_URB, UC_URB, RAINR_URB, RAINB_URB, RAING_URB, ROFF_URB, SFC_TEMP, ALBEDO, MWFLX, MUFLX, MVFLX, SHFLX, LHFLX, GHFLX, Z0M, Z0H, Z0E, U10, V10, T2, Q2)
Main routine for land submodel.
Definition: scale_urban_dyn_kusaka01.F90:241
scale_prof
module profiler
Definition: scale_prof.F90:11
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:33
scale_cpl_sfc_index::i_r_nir
integer, parameter, public i_r_nir
Definition: scale_cpl_sfc_index.F90:30
scale_urban_dyn_kusaka01::urban_dyn_kusaka01_setup
subroutine, public urban_dyn_kusaka01_setup(UIA, UIS, UIE, UJA, UJS, UJE, Z0M, Z0H, Z0E)
Setup.
Definition: scale_urban_dyn_kusaka01.F90:118
scale_precision
module PRECISION
Definition: scale_precision.F90:14
scale_const::const_pi
real(rp), public const_pi
pi
Definition: scale_const.F90:31
scale_cpl_sfc_index::i_r_ir
integer, parameter, public i_r_ir
Definition: scale_cpl_sfc_index.F90:29
scale_io
module STDIO
Definition: scale_io.F90:10
scale_cpl_sfc_index::i_r_diffuse
integer, parameter, public i_r_diffuse
Definition: scale_cpl_sfc_index.F90:38
scale_precision::rp
integer, parameter, public rp
Definition: scale_precision.F90:41
scale_file_history
module file_history
Definition: scale_file_history.F90:15
Generated by   doxygen 1.8.13