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