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