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