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