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  use scale_mapprojection, only: &
22  base_lon => mapprojection_basepoint_lon
23  !-----------------------------------------------------------------------------
24  implicit none
25  private
26  !-----------------------------------------------------------------------------
27  !
28  !++ Public procedure
29  !
30  public :: urban_dyn_kusaka01_setup
31  public :: urban_dyn_kusaka01_finalize
32  public :: urban_dyn_kusaka01
33 
34  !-----------------------------------------------------------------------------
35  !
36  !++ Public parameters & variables
37  !
38  !-----------------------------------------------------------------------------
39  !
40  !++ Private procedure
41  !
42  private :: slc_main
43  private :: canopy_wind
44  private :: cal_beta
45  private :: cal_psi
46  private :: mos
47  private :: multi_layer
48  private :: urban_param_setup
49  private :: read_urban_param_table
50  private :: read_urban_gridded_data_2d
51  private :: read_urban_gridded_data_3d
52 
53  !-----------------------------------------------------------------------------
54  !
55  !++ Private parameters & variables
56  !
57  !-----------------------------------------------------------------------------
58  !
59  ! from namelist
60  real(RP), private :: DTS_MAX = 0.1_rp ! maximum dT during one step [K/step]
61  ! DTS_MAX * dt
62  integer, private :: BOUND = 1 ! Boundary Condition for Roof, Wall, Ground Layer Temp
63  ! [1: Zero-Flux, 2: T = Constant]
64 #ifdef DEBUG_KUSAKA01
65  logical, private :: debug = .true.
66 #else
67  logical, private :: debug = .false.
68 #endif
69  !$acc declare create(DTS_MAX,BOUND,debug)
70 
71  ! urban parameters
72  real(RP), private :: ZR = 10.0_rp ! roof level (building height) [m]
73  real(RP), private :: roof_width = 9.0_rp ! roof width [m]
74  real(RP), private :: road_width = 11.0_rp ! road width [m]
75  real(RP), private :: SIGMA_ZED = 1.0_rp ! Standard deviation of roof height [m]
76  real(RP), private :: AH_TBL = 17.5_rp ! Sensible Anthropogenic heat from urban subgrid [W/m^2]
77  real(RP), private :: AHL_TBL = 0.0_rp ! Latent Anthropogenic heat from urban subgrid [W/m^2]
78  real(RP), private :: BETR_CONST = -1.0_rp ! Evaporation efficiency of roof [-]
79  real(RP), private :: BETB_CONST = -1.0_rp ! of building [-]
80  real(RP), private :: BETG_CONST = -1.0_rp ! of ground [-]
81  real(RP), private :: STRGR = 0.0_rp ! rain strage on roof [-]
82  real(RP), private :: STRGB = 0.0_rp ! on wall [-]
83  real(RP), private :: STRGG = 0.0_rp ! on ground [-]
84  real(RP), private :: CAPR = 1.2e6_rp ! heat capacity of roof [J m-3 K]
85  real(RP), private :: CAPB = 1.2e6_rp ! of wall [J m-3 K]
86  real(RP), private :: CAPG = 1.2e6_rp ! of ground [J m-3 K]
87  real(RP), private :: AKSR = 2.28_rp ! thermal conductivity of roof [W m-1 K]
88  real(RP), private :: AKSB = 2.28_rp ! of wall [W m-1 K]
89  real(RP), private :: AKSG = 2.28_rp ! of ground [W m-1 K]
90  real(RP), private :: ALBR = 0.2_rp ! surface albedo of roof
91  real(RP), private :: ALBB = 0.2_rp ! surface albedo of wall
92  real(RP), private :: ALBG = 0.2_rp ! surface albedo of ground
93  real(RP), private :: EPSR = 0.90_rp ! Surface emissivity of roof
94  real(RP), private :: EPSB = 0.90_rp ! Surface emissivity of wall
95  real(RP), private :: EPSG = 0.90_rp ! Surface emissivity of ground
96  real(RP), private :: Z0R = 0.01_rp ! roughness length for momentum of building roof
97  real(RP), private :: Z0B = 0.0001_rp ! roughness length for momentum of building wall
98  real(RP), private :: Z0G = 0.01_rp ! roughness length for momentum of ground
99  real(RP), private :: TRLEND = 293.00_rp ! lower boundary condition of roof temperature [K]
100  real(RP), private :: TBLEND = 293.00_rp ! lower boundary condition of wall temperature [K]
101  real(RP), private :: TGLEND = 293.00_rp ! lower boundary condition of ground temperature [K]
102  !$acc declare create(ZR,BETR_CONST,BETB_CONST,BETG_CONST,STRGR,STRGB,STRGG,CAPR,CAPB,CAPG,AKSR,AKSB,AKSG,ALBR,ALBB,ALBG,EPSR,EPSB,EPSG,Z0R,TRLEND,TBLEND,TGLEND)
103 
104  ! calculated in subroutine urban_param_set
105  real(RP), private :: R ! Normalized roof wight (eq. building coverage ratio)
106  real(RP), private :: RW ! (= 1 - R)
107  real(RP), private :: HGT ! Normalized building height
108  real(RP), private :: Z0HR ! roughness length for heat of roof
109  real(RP), private :: Z0HB ! roughness length for heat of building wall
110  real(RP), private :: Z0HG ! roughness length for heat of ground
111  real(RP), private :: Z0C_TBL ! Roughness length above canyon for momentum [m]
112  real(RP), private :: Z0HC_TBL ! Roughness length above canyon for heat [m]
113  real(RP), private :: ZDC_TBL ! Displacement height [m]
114  real(RP), private :: SVF ! Sky view factor [-]
115  !$acc declare create(R,RW,HGT,Z0HR,Z0HB,Z0HG,SVF)
116 
117  ! history
118  integer, private :: I_SHR, I_SHB, I_SHG
119  integer, private :: I_LHR, I_LHB, I_LHG
120  integer, private :: I_GHR, I_GHB, I_GHG
121  integer, private :: I_RNR, I_RNB, I_RNG, I_RNgrd
122 
123 #ifdef DEBUG_KUSAKA01
124  integer, private :: cnt_num1, cnt_num2
125  integer, private :: cnt_itr1, cnt_itr2
126  integer, private :: max_itr1, max_itr2
127  !$acc declare create(cnt_num1,cnt_num2,cnt_itr1,cnt_itr2,max_itr1,max_itr2)
128 #endif
129 
130  !-----------------------------------------------------------------------------
131 contains
132  !-----------------------------------------------------------------------------
134  subroutine urban_dyn_kusaka01_setup( &
135  UIA, UIS, UIE, UJA, UJS, UJE, &
136  fact_urban, &
137  Z0M, Z0H, Z0E, ZD, &
138  AH_URB, AHL_URB, AH_TOFFSET )
139  use scale_prc, only: &
140  prc_myrank, &
141  prc_abort
142  use scale_const, only: &
143  undef => const_undef
144  use scale_file_history, only: &
145  file_history_reg
146  implicit none
147 
148  integer, intent(in) :: uia, uis, uie
149  integer, intent(in) :: uja, ujs, uje
150  real(rp), intent(in) :: fact_urban(uia,uja)
151  real(rp), intent(out) :: z0m(uia,uja)
152  real(rp), intent(out) :: z0h(uia,uja)
153  real(rp), intent(out) :: z0e(uia,uja)
154  real(rp), intent(out) :: zd (uia,uja)
155  real(rp), intent(out) :: ah_urb (uia,uja,1:24)
156  real(rp), intent(out) :: ahl_urb (uia,uja,1:24)
157  real(rp), intent(out) :: ah_toffset
158 
159  character(len=H_LONG) :: urban_dyn_kusaka01_param_in_filename = ''
160  character(len=H_LONG) :: urban_dyn_kusaka01_gridded_z0m_in_filename = ''
161  character(len=H_LONG) :: urban_dyn_kusaka01_gridded_z0m_in_varname = 'URBAN_Z0M'
162  character(len=H_LONG) :: urban_dyn_kusaka01_gridded_z0h_in_filename = ''
163  character(len=H_LONG) :: urban_dyn_kusaka01_gridded_z0h_in_varname = 'URBAN_Z0H'
164  !character(len=H_LONG) :: URBAN_DYN_KUSAKA01_GRIDDED_ZD_IN_FILENAME = '' !< gridded data of ZD
165  !character(len=H_LONG) :: URBAN_DYN_KUSAKA01_GRIDDED_ZD_IN_VARNAME = 'URBAN_ZD' !< var name of gridded data for Zd
166  character(len=H_LONG) :: urban_dyn_kusaka01_gridded_ah_in_filename = ''
167  character(len=H_LONG) :: urban_dyn_kusaka01_gridded_ah_in_varname = 'URBAN_AH'
168  character(len=H_LONG) :: urban_dyn_kusaka01_gridded_ahl_in_filename = ''
169  character(len=H_LONG) :: urban_dyn_kusaka01_gridded_ahl_in_varname = 'URBAN_AHL'
170 
171  namelist / param_urban_dyn_kusaka01 / &
172  dts_max, &
173  bound, &
174  debug, &
175  urban_dyn_kusaka01_param_in_filename, &
176  urban_dyn_kusaka01_gridded_z0m_in_filename, &
177  urban_dyn_kusaka01_gridded_z0h_in_filename, &
178  !URBAN_DYN_KUSAKA01_GRIDDED_ZD_IN_FILENAME, &
179  urban_dyn_kusaka01_gridded_ah_in_filename, &
180  urban_dyn_kusaka01_gridded_ahl_in_filename
181 
182  real(rp) :: udata(uia,uja)
183  real(rp) :: udata2(uia,uja,24)
184  real(rp) :: rtime
185  integer :: itime
186 
187  real(rp) :: ahdiurnal(1:24) ! AH diurnal profile
188 
189  integer :: i, j, k
190  integer :: ierr
191  !---------------------------------------------------------------------------
192 
193  log_newline
194  log_info("URBAN_DYN_kusaka01_setup",*) 'Setup'
195 
196  !$acc data copyin(fact_urban) &
197  !$acc copyout(Z0M, Z0H, Z0E, ZD, AH_URB, AHL_URB)
198 
199  !--- read namelist
200  rewind(io_fid_conf)
201  read(io_fid_conf,nml=param_urban_dyn_kusaka01,iostat=ierr)
202  if( ierr < 0 ) then !--- missing
203  log_info("URBAN_DYN_kusaka01_setup",*) 'Not found namelist. Default used.'
204  elseif( ierr > 0 ) then !--- fatal error
205  log_error("URBAN_DYN_kusaka01_setup",*) 'Not appropriate names in namelist PARAM_URBAN_DYN_KUSAKA01. Check!'
206  call prc_abort
207  endif
208  log_nml(param_urban_dyn_kusaka01)
209 
210  !-- read urban parameter from file
211  if( urban_dyn_kusaka01_param_in_filename /= '' ) then
212  call read_urban_param_table( trim(urban_dyn_kusaka01_param_in_filename) )
213  endif
214 
215  ! set other urban parameters
216  call urban_param_setup
217 
218  ! Local time: 1 to 24
219  ahdiurnal(:) = (/ 0.356, 0.274, 0.232, 0.251, 0.375, 0.647, 0.919, 1.135, 1.249, 1.328, &
220  1.365, 1.363, 1.375, 1.404, 1.457, 1.526, 1.557, 1.521, 1.372, 1.206, &
221  1.017, 0.876, 0.684, 0.512 /)
222 
223  ! Shift to UTC based on local solar timezone of domain center
224  rtime = modulo(base_lon, 360.0_rp) / 15.0_rp
225  itime = nint(rtime)
226  ah_toffset = rtime - itime ! currently not used: difference of time from domain center
227  ahdiurnal(:) = cshift(ahdiurnal(:), -1*itime) ! convert from LT to UTC
228 
229  do j = ujs, uje
230  do i = uis, uie
231  z0m(i,j) = z0c_tbl
232  z0h(i,j) = z0hc_tbl
233  z0e(i,j) = z0hc_tbl
234  zd(i,j) = zdc_tbl
235  do k = 1, 24 ! UTC
236  ah_urb(i,j,k) = ah_tbl * ahdiurnal(k) * fact_urban(i,j)
237  ahl_urb(i,j,k) = ahl_tbl * ahdiurnal(k) * fact_urban(i,j)
238  enddo
239  enddo
240  enddo
241 
242  !-- replace gridded Z0M data if there is a file
243  if( urban_dyn_kusaka01_gridded_z0m_in_filename /= '' ) then
244  udata(:,:) = 0.0_rp
245  call read_urban_gridded_data_2d( &
246  uia, uja, &
247  urban_dyn_kusaka01_gridded_z0m_in_filename, &
248  urban_dyn_kusaka01_gridded_z0m_in_varname, &
249  udata(:,:) )
250 
251  ! replace to gridded data
252  do j = ujs, uje
253  do i = uis, uie
254  if( udata(i,j) /= undef )then
255  if ( udata(i,j) > 0.0_rp ) then
256  z0m(i,j) = udata(i,j)
257  else if ( udata(i,j) < 0.0_rp ) then
258  log_error("URBAN_DYN_kusaka01_setup",*) 'Gridded Z0M data includes data less than 0. Please check data!',i,j
259  call prc_abort
260  else ! Z0M = 0[m]
261  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
262  endif
263  endif
264  enddo
265  enddo
266  endif
267 
268  !-- read gridded Z0H & Z0E data from a file
269  if( urban_dyn_kusaka01_gridded_z0h_in_filename /= '' ) then
270  udata = 0.0_rp
271  call read_urban_gridded_data_2d( &
272  uia, uja, &
273  urban_dyn_kusaka01_gridded_z0h_in_filename, &
274  urban_dyn_kusaka01_gridded_z0h_in_varname, &
275  udata )
276 
277  ! replace to gridded data
278  do j = ujs, uje
279  do i = uis, uie
280  if( udata(i,j) /= undef )then
281  if ( udata(i,j) > 0.0_rp ) then
282  z0h(i,j) = udata(i,j)
283  z0e(i,j) = udata(i,j)
284  else if ( udata(i,j) < 0.0_rp ) then
285  log_error("URBAN_DYN_kusaka01_setup",*) 'Gridded Z0H data includes data less than 0. Please check data!',i,j
286  call prc_abort
287  else ! Z0H = 0[m]
288  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
289  endif
290  endif
291  enddo
292  enddo
293  endif
294 
295  ! currently NOT USED
296  !-- read gridded ZD data from a file
297  !if( URBAN_DYN_KUSAKA01_GRIDDED_ZD_IN_FILENAME /= '' ) then
298  ! udata = 0.0_RP
299  ! call read_urban_gridded_data_2D( &
300  ! UIA, UJA, &
301  ! URBAN_DYN_KUSAKA01_GRIDDED_ZD_IN_FILENAME, &
302  ! URBAN_DYN_KUSAKA01_GRIDDED_ZD_IN_VARNAME, &
303  ! udata )
304  !
305  ! ! replace to gridded data
306  ! do j = UJS, UJE
307  ! do i = UIS, UIE
308  ! if( udata(i,j) /= UNDEF )then
309  ! if ( udata(i,j) >= 0.0_RP ) then
310  ! ZD(i,j) = udata(i,j)
311  ! else
312  ! LOG_ERROR("URBAN_DYN_kusaka01_setup",*) 'Gridded ZD data includes data less than 0. Please check data!',i,j
313  ! call PRC_abort
314  ! endif
315  ! endif
316  ! enddo
317  ! enddo
318  !endif
319 
320  !-- read gridded AH data from a file
321  if( urban_dyn_kusaka01_gridded_ah_in_filename /= '' ) then
322  udata2 = 0.0_rp
323  call read_urban_gridded_data_3d( &
324  uia, uja, &
325  urban_dyn_kusaka01_gridded_ah_in_filename, &
326  urban_dyn_kusaka01_gridded_ah_in_varname, &
327  udata2 )
328 
329  ! replace to gridded data
330  do k = 1, 24
331  do j = ujs, uje
332  do i = uis, uie
333  if( udata2(i,j,k) /= undef )then
334  ah_urb(i,j,k) = udata2(i,j,k)
335  endif
336  enddo
337  enddo
338  enddo
339  endif
340 
341  !-- read gridded AHL data from a file
342  if( urban_dyn_kusaka01_gridded_ahl_in_filename /= '' ) then
343  udata2 = 0.0_rp
344  call read_urban_gridded_data_3d( &
345  uia, uja, &
346  urban_dyn_kusaka01_gridded_ahl_in_filename, &
347  urban_dyn_kusaka01_gridded_ahl_in_varname, &
348  udata2 )
349 
350  ! replace to gridded data
351  do k = 1, 24
352  do j = ujs, uje
353  do i = uis, uie
354  if( udata2(i,j,k) /= undef )then
355  ahl_urb(i,j,k) = udata2(i,j,k)
356  endif
357  enddo
358  enddo
359  enddo
360  endif
361 
362  call file_history_reg( 'URBAN_SHR', 'urban sensible heat flux on roof', 'W/m2', i_shr , ndims=2 )
363  call file_history_reg( 'URBAN_SHB', 'urban sensible heat flux on wall', 'W/m2', i_shb , ndims=2 )
364  call file_history_reg( 'URBAN_SHG', 'urban sensible heat flux on road', 'W/m2', i_shg , ndims=2 )
365  call file_history_reg( 'URBAN_LHR', 'urban latent heat flux on roof', 'W/m2', i_lhr , ndims=2 )
366  call file_history_reg( 'URBAN_LHB', 'urban latent heat flux on wall', 'W/m2', i_lhb , ndims=2 )
367  call file_history_reg( 'URBAN_LHG', 'urban latent heat flux on road', 'W/m2', i_lhg , ndims=2 )
368  call file_history_reg( 'URBAN_GHR', 'urban ground heat flux on roof', 'W/m2', i_ghr , ndims=2 )
369  call file_history_reg( 'URBAN_GHB', 'urban ground heat flux on wall', 'W/m2', i_ghb , ndims=2 )
370  call file_history_reg( 'URBAN_GHG', 'urban ground heat flux on road', 'W/m2', i_ghg , ndims=2 )
371  call file_history_reg( 'URBAN_RNR', 'urban net radiation on roof', 'W/m2', i_rnr , ndims=2 )
372  call file_history_reg( 'URBAN_RNB', 'urban net radiation on wall', 'W/m2', i_rnb , ndims=2 )
373  call file_history_reg( 'URBAN_RNG', 'urban net radiation on road', 'W/m2', i_rng , ndims=2 )
374  call file_history_reg( 'URBAN_RNgrd', 'urban grid average of net radiation', 'W/m2', i_rngrd, ndims=2 )
375 
376  !$acc update device(Z0M, Z0H, Z0E, ZD, AH_URB, AHL_URB)
377  !$acc end data
378 
379  !$acc update device(DTS_MAX,BOUND,debug)
380  !$acc update device(ZR,BETR_CONST,BETB_CONST,BETG_CONST,STRGR,STRGB,STRGG,CAPR,CAPB,CAPG,AKSR,AKSB,AKSG,ALBR,ALBB,ALBG,EPSR,EPSB,EPSG,Z0R,TRLEND,TBLEND,TGLEND)
381  !$acc update device(R,RW,HGT,Z0HR,Z0HB,Z0HG,SVF)
382 
383 #ifdef DEBUG_KUSAKA01
384  cnt_num1 = 0
385  cnt_num2 = 0
386  cnt_itr1 = 0
387  cnt_itr2 = 0
388  max_itr1 = 0
389  max_itr2 = 0
390  !$acc update device(cnt_num1,cnt_num2,cnt_itr1,cnt_itr2,max_itr1,max_itr2)
391 #endif
392 
393  return
394  end subroutine urban_dyn_kusaka01_setup
395 
396  !-----------------------------------------------------------------------------
398  subroutine urban_dyn_kusaka01_finalize
399  use mpi
400  use scale_prc, only: &
401  prc_local_comm_world, &
402  prc_ismaster
403  implicit none
404 
405 #ifdef DEBUG_KUSAKA01
406  integer :: iwork1(4), iwork2(2), ierr
407  !$acc update host(cnt_num1,cnt_num2,cnt_itr1,cnt_itr2,max_itr1,max_itr2)
408  iwork1(1) = cnt_num1
409  iwork1(2) = cnt_num2
410  iwork1(3) = cnt_itr1
411  iwork1(4) = cnt_itr2
412  call mpi_reduce( mpi_in_place, iwork1, 4, mpi_integer, mpi_sum, 0, prc_local_comm_world, ierr)
413  iwork2(1) = max_itr1
414  iwork2(2) = max_itr2
415  call mpi_reduce( mpi_in_place, iwork2, 2, mpi_integer, mpi_max, 0, prc_local_comm_world, ierr)
416 
417  if ( prc_ismaster ) then
418  log_info("URBAN_DYN_kusaka01_finalize",*) "Averaged iteration count"
419  log_info_cont(*) "TR: ", real(iwork1(3),dp) / iwork1(1), ", (max ", iwork2(1), ")"
420  log_info_cont(*) "TB/TG: ", real(iwork1(4),dp) / iwork1(2), ", (max ", iwork2(2), ")"
421  end if
422 #endif
423 
424  end subroutine urban_dyn_kusaka01_finalize
425 
426  !-----------------------------------------------------------------------------
428 
429  subroutine urban_dyn_kusaka01( &
430  UKA, UKS, UKE, UIA, UIS, UIE, UJA, UJS, UJE, &
431  TMPA, PRSA, &
432  U1, V1, &
433  DENS, QA, LHV, &
434  Z1, &
435  RHOS, PRSS, &
436  LWD, SWD, &
437  RAIN, EFLX, &
438  Z0M, Z0H, Z0E, &
439  ZD, &
440  CDZ, &
441  TanSL_X, TanSL_Y, &
442  fact_urban, &
443  dt, &
444  TRL_URB, TBL_URB, TGL_URB, &
445  TR_URB, TB_URB, TG_URB, &
446  TC_URB, QC_URB, UC_URB, &
447  RAINR_URB, RAINB_URB, RAING_URB, &
448  ROFF_URB, &
449  SFC_TEMP, &
450  ALBEDO, &
451  MWFLX, MUFLX, MVFLX, &
452  SHFLX, LHFLX, GHFLX, &
453  Ustar, Tstar, Qstar, Wstar, &
454  RLmo, &
455  U10, V10, T2, Q2 )
456  use scale_prc, only: &
457  prc_myrank, &
458  prc_abort
459  use scale_const, only: &
460  eps => const_eps, &
461  undef => const_undef, &
462  rdry => const_rdry, &
463  rvap => const_rvap
464  use scale_bulkflux, only: &
465  bulkflux
466  implicit none
467  integer, intent(in) :: uka, uks, uke
468  integer, intent(in) :: uia, uis, uie
469  integer, intent(in) :: uja, ujs, uje
470 
471  real(rp), intent(in) :: tmpa(uia,uja)
472  real(rp), intent(in) :: prsa(uia,uja)
473  real(rp), intent(in) :: u1 (uia,uja)
474  real(rp), intent(in) :: v1 (uia,uja)
475  real(rp), intent(in) :: dens(uia,uja)
476  real(rp), intent(in) :: qa (uia,uja)
477  real(rp), intent(in) :: lhv (uia,uja)
478  real(rp), intent(in) :: z1 (uia,uja)
479  real(rp), intent(in) :: rhos(uia,uja) ! density at the surface [kg/m3]
480  real(rp), intent(in) :: prss(uia,uja)
481  real(rp), intent(in) :: lwd (uia,uja,2)
482  real(rp), intent(in) :: swd (uia,uja,2)
483  real(rp), intent(in) :: rain(uia,uja)
484  real(rp), intent(in) :: eflx(uia,uja)
485  real(rp), intent(in) :: z0m (uia,uja)
486  real(rp), intent(in) :: z0h (uia,uja)
487  real(rp), intent(in) :: z0e (uia,uja)
488  real(rp), intent(in) :: zd (uia,uja)
489  real(rp), intent(in) :: cdz(uka)
490  real(rp), intent(in) :: tansl_x(uia,uja)
491  real(rp), intent(in) :: tansl_y(uia,uja)
492  real(rp), intent(in) :: fact_urban(uia,uja)
493  real(dp), intent(in) :: dt
494 
495  real(rp), intent(inout) :: tr_urb (uia,uja)
496  real(rp), intent(inout) :: tb_urb (uia,uja)
497  real(rp), intent(inout) :: tg_urb (uia,uja)
498  real(rp), intent(inout) :: tc_urb (uia,uja)
499  real(rp), intent(inout) :: qc_urb (uia,uja)
500  real(rp), intent(inout) :: uc_urb (uia,uja)
501  real(rp), intent(inout) :: trl_urb (uks:uke,uia,uja)
502  real(rp), intent(inout) :: tbl_urb (uks:uke,uia,uja)
503  real(rp), intent(inout) :: tgl_urb (uks:uke,uia,uja)
504  real(rp), intent(inout) :: rainr_urb(uia,uja)
505  real(rp), intent(inout) :: rainb_urb(uia,uja)
506  real(rp), intent(inout) :: raing_urb(uia,uja)
507  real(rp), intent(out) :: roff_urb (uia,uja)
508 
509  real(rp), intent(out) :: sfc_temp(uia,uja)
510  real(rp), intent(out) :: albedo (uia,uja,n_rad_dir,n_rad_rgn)
511  real(rp), intent(out) :: mwflx (uia,uja)
512  real(rp), intent(out) :: muflx (uia,uja)
513  real(rp), intent(out) :: mvflx (uia,uja)
514  real(rp), intent(out) :: shflx (uia,uja)
515  real(rp), intent(out) :: lhflx (uia,uja)
516  real(rp), intent(out) :: ghflx (uia,uja)
517  real(rp), intent(out) :: ustar (uia,uja)
518  real(rp), intent(out) :: tstar (uia,uja)
519  real(rp), intent(out) :: qstar (uia,uja)
520  real(rp), intent(out) :: wstar (uia,uja)
521  real(rp), intent(out) :: rlmo (uia,uja)
522  real(rp), intent(out) :: u10 (uia,uja)
523  real(rp), intent(out) :: v10 (uia,uja)
524  real(rp), intent(out) :: t2 (uia,uja)
525  real(rp), intent(out) :: q2 (uia,uja)
526 
527  ! parameter
528  logical, parameter :: lsolar = .false. ! [true=both, false=SSG only]
529 
530  real(rp), parameter :: uabs_min = 0.1_rp
531 
532  ! work
533  real(rp) :: tr
534  real(rp) :: tb
535  real(rp) :: tg
536  real(rp) :: tc
537  real(rp) :: qc
538  real(rp) :: uc
539  real(rp) :: trl(uks:uke)
540  real(rp) :: tbl(uks:uke)
541  real(rp) :: tgl(uks:uke)
542  real(rp) :: rainr
543  real(rp) :: rainb
544  real(rp) :: raing
545  real(rp) :: albd_lw
546  real(rp) :: albd_sw
547 
548  real(rp) :: shr (uia,uja)
549  real(rp) :: shb (uia,uja)
550  real(rp) :: shg (uia,uja)
551  real(rp) :: lhr (uia,uja)
552  real(rp) :: lhb (uia,uja)
553  real(rp) :: lhg (uia,uja)
554  real(rp) :: ghr (uia,uja)
555  real(rp) :: ghb (uia,uja)
556  real(rp) :: ghg (uia,uja)
557  real(rp) :: rnr (uia,uja)
558  real(rp) :: rnb (uia,uja)
559  real(rp) :: rng (uia,uja)
560  real(rp) :: rngrd(uia,uja)
561 
562  real(rp) :: dzr(uka) ! thickness of each roof layer [m]
563  real(rp) :: dzb(uka) ! thickness of each building layer [m]
564  real(rp) :: dzg(uka) ! thickness of each road layer [m]
565 
566  real(rp) :: swdt(2)
567  real(rp) :: lwdt(2)
568 
569  real(rp) :: uabs ! modified absolute velocity [m/s]
570 
571  real(rp) :: mflx
572  real(rp) :: w
573 
574  ! work
575  real(rp) :: trlp(uka)
576  real(rp) :: tblp(uka)
577  real(rp) :: tglp(uka)
578  real(rp) :: a(uka)
579  real(rp) :: b(uka)
580  real(rp) :: c(uka)
581  real(rp) :: d(uka)
582  real(rp) :: p(uka)
583  real(rp) :: q(uka)
584 
585  logical :: converged
586 
587 #ifdef DEBUG_KUSAKA01
588  integer :: cn1, cn2, ci1, ci2, mi1, mi2
589 #endif
590 
591  integer :: k, i, j
592  !---------------------------------------------------------------------------
593 
594  log_progress(*) 'urban / dynamics / Kusaka01'
595 
596  !$acc data copyin(TMPA,PRSA,U1,V1,DENS,QA,LHV,Z1,RHOS,PRSS,LWD,SWD,RAIN,EFLX,Z0M,Z0H,Z0E,ZD,CDZ,TanSL_X,TanSL_Y,fact_urban) &
597  !$acc copy(TR_URB,TB_URB,TG_URB,TC_URB,QC_URB,UC_URB,TRL_URB,TBL_URB,TGL_URB,RAINR_URB,RAINB_URB,RAING_URB) &
598  !$acc copyout(ROFF_URB,SFC_TEMP,ALBEDO,MWFLX,MUFLX,MVFLX,SHFLX,LHFLX,GHFLX,Ustar,Tstar,Qstar,Wstar,RLmo,U10,V10,T2,Q2) &
599  !$acc create(SHR,SHB,SHG,LHR,LHB,LHG,GHR,GHB,GHG,RNR,RNB,RNG,RNgrd,DZR,DZB,DZG)
600 
601 
602  !$acc kernels
603  dzr(:) = cdz(:)
604  dzb(:) = cdz(:)
605  dzg(:) = cdz(:)
606  !$acc end kernels
607 
608  converged = .true.
609 
610 #ifdef DEBUG_KUSAKA01
611  cn1 = 0; cn2 = 0; ci1 = 0; ci2 = 0; mi1 = 0; mi2 = 0
612 #endif
613 
614  !$omp parallel do schedule(dynamic) collapse(2) &
615 #ifdef DEBUG_KUSAKA01
616  !$omp reduction(+:cn1,cn2,ci1,ci2) reduction(max:mi1,mi2) &
617 #endif
618  !$omp private(w,Uabs,TR,TB,TG,TC,QC,UC,TRL,TBL,TGL,RAINR,RAINB,RAING,ALBD_LW,ALBD_SW,MFLX,SWDt,LWDt, &
619  !$omp TRLP,TBLP,TGLP,A,B,C,D,P,Q)
620  !$acc parallel
621  !$acc loop gang collapse(2) reduction(.and.: converged) independent &
622 #ifdef DEBUG_KUSAKA01
623  !$acc reduction(+:cn1,cn2,ci1,ci2) reduction(max:mi1,mi2) &
624 #endif
625  !$acc private(w,Uabs,TR,TB,TG,TC,QC,UC,TRL,TBL,TGL,RAINR,RAINB,RAING,ALBD_LW,ALBD_SW,MFLX,SWDt,LWDt, &
626  !$acc TRLP,TBLP,TGLP,A,B,C,D,P,Q)
627  do j = ujs, uje
628  do i = uis, uie
629 
630  if( fact_urban(i,j) > 0.0_rp ) then
631 
632  w = u1(i,j) * tansl_x(i,j) + v1(i,j) * tansl_y(i,j)
633  uabs = sqrt( u1(i,j)**2 + v1(i,j)**2 + w**2 )
634 
635  ! save
636  tr = tr_urb(i,j)
637  tb = tb_urb(i,j)
638  tg = tg_urb(i,j)
639  tc = tc_urb(i,j)
640  qc = qc_urb(i,j)
641  uc = uc_urb(i,j)
642 
643  do k = uks, uke
644  trl(k) = trl_urb(k,i,j)
645  tbl(k) = tbl_urb(k,i,j)
646  tgl(k) = tgl_urb(k,i,j)
647  end do
648 
649  rainr = rainr_urb(i,j)
650  rainb = rainb_urb(i,j)
651  raing = raing_urb(i,j)
652 
653  swdt(:) = swd(i,j,:)
654  lwdt(:) = lwd(i,j,:)
655 
656  call slc_main( uka, uks, uke, &
657  trl(:), & ! [INOUT]
658  tbl(:), & ! [INOUT]
659  tgl(:), & ! [INOUT]
660  tr, & ! [INOUT]
661  tb, & ! [INOUT]
662  tg, & ! [INOUT]
663  tc, & ! [INOUT]
664  qc, & ! [INOUT]
665  uc, & ! [INOUT]
666  rainr, & ! [INOUT]
667  rainb, & ! [INOUT]
668  raing, & ! [INOUT]
669  roff_urb(i,j), & ! [OUT]
670  albd_lw, & ! [OUT]
671  albd_sw, & ! [OUT]
672  shr(i,j), & ! [OUT]
673  shb(i,j), & ! [OUT]
674  shg(i,j), & ! [OUT]
675  lhr(i,j), & ! [OUT]
676  lhb(i,j), & ! [OUT]
677  lhg(i,j), & ! [OUT]
678  ghr(i,j), & ! [OUT]
679  ghb(i,j), & ! [OUT]
680  ghg(i,j), & ! [OUT]
681  rnr(i,j), & ! [OUT]
682  rnb(i,j), & ! [OUT]
683  rng(i,j), & ! [OUT]
684  sfc_temp(i,j), & ! [OUT]
685  rngrd(i,j), & ! [OUT]
686  shflx(i,j), & ! [OUT]
687  lhflx(i,j), & ! [OUT]
688  ghflx(i,j), & ! [OUT]
689  mflx, & ! [OUT]
690  ustar(i,j), & ! [OUT]
691  tstar(i,j), & ! [OUT]
692  qstar(i,j), & ! [OUT]
693  u10(i,j), & ! [OUT]
694  v10(i,j), & ! [OUT]
695  t2(i,j), & ! [OUT]
696  q2(i,j), & ! [OUT]
697  converged, & ! [OUT]
698  lsolar, & ! [IN]
699  prsa(i,j), & ! [IN]
700  prss(i,j), & ! [IN]
701  tmpa(i,j), & ! [IN]
702  rhos(i,j), & ! [IN]
703  qa(i,j), & ! [IN]
704  uabs, & ! [IN]
705  u1(i,j), & ! [IN]
706  v1(i,j), & ! [IN]
707  lhv(i,j), & ! [IN]
708  z1(i,j), & ! [IN]
709  swdt(:), & ! [IN]
710  lwdt(:), & ! [IN]
711  rain(i,j), & ! [IN]
712  eflx(i,j), & ! [IN]
713  dens(i,j), & ! [IN]
714  z0m(i,j), & ! [IN]
715  z0h(i,j), & ! [IN]
716  zd(i,j), & ! [IN]
717  dzr(:), & ! [IN]
718  dzg(:), dzb(:), & ! [IN]
719  dt, & ! [IN]
720  i, j, & ! [IN]
721 #ifdef DEBUG_KUSAKA01
722  cn1, cn2, ci1, ci2, mi1, mi2, & ! [INOUT]
723 #endif
724  trlp(:), tblp(:), tglp(:), & ! (work)
725  a(:), b(:), c(:), d(:), & ! (work)
726  p(:), q(:) ) ! (work)
727 
728  ! update
729  tr_urb(i,j) = tr
730  tb_urb(i,j) = tb
731  tg_urb(i,j) = tg
732  tc_urb(i,j) = tc
733  qc_urb(i,j) = qc
734  uc_urb(i,j) = uc
735  do k = uks, uke
736  trl_urb(k,i,j) = trl(k)
737  tbl_urb(k,i,j) = tbl(k)
738  tgl_urb(k,i,j) = tgl(k)
739  end do
740  rainr_urb(i,j) = rainr
741  rainb_urb(i,j) = rainb
742  raing_urb(i,j) = raing
743 
744  albedo(i,j,i_r_direct ,i_r_ir ) = albd_lw
745  albedo(i,j,i_r_diffuse,i_r_ir ) = albd_lw
746  albedo(i,j,i_r_direct ,i_r_nir) = albd_sw
747  albedo(i,j,i_r_diffuse,i_r_nir) = albd_sw
748  albedo(i,j,i_r_direct ,i_r_vis) = albd_sw
749  albedo(i,j,i_r_diffuse,i_r_vis) = albd_sw
750 
751  if ( uabs < eps ) then
752  mwflx(i,j) = 0.0_rp
753  muflx(i,j) = 0.0_rp
754  mvflx(i,j) = 0.0_rp
755  else
756  mwflx(i,j) = mflx * w / uabs
757  muflx(i,j) = mflx * u1(i,j) / uabs
758  mvflx(i,j) = mflx * v1(i,j) / uabs
759  end if
760 
761  else
762  sfc_temp(i,j) = undef
763  albedo(i,j,:,:) = undef
764  mwflx(i,j) = undef
765  muflx(i,j) = undef
766  mvflx(i,j) = undef
767  shflx(i,j) = undef
768  lhflx(i,j) = undef
769  ghflx(i,j) = undef
770  ustar(i,j) = undef
771  tstar(i,j) = undef
772  qstar(i,j) = undef
773  wstar(i,j) = undef
774  rlmo(i,j) = undef
775  u10(i,j) = undef
776  v10(i,j) = undef
777  t2(i,j) = undef
778  q2(i,j) = undef
779  shr(i,j) = undef
780  shb(i,j) = undef
781  shg(i,j) = undef
782  lhr(i,j) = undef
783  lhb(i,j) = undef
784  lhg(i,j) = undef
785  ghr(i,j) = undef
786  ghb(i,j) = undef
787  ghg(i,j) = undef
788  rnr(i,j) = undef
789  rnb(i,j) = undef
790  rng(i,j) = undef
791  rngrd(i,j) = undef
792  endif
793 
794  end do
795  end do
796  !$acc end parallel
797 
798  if ( .not. converged ) then
799  log_error("URBAN_DYN_kusaka01_SLC_main",*) "not converged"
800  call prc_abort
801  end if
802 
803  call put_history( uia, uja, &
804  shr(:,:), shb(:,:), shg(:,:), &
805  lhr(:,:), lhb(:,:), lhg(:,:), &
806  ghr(:,:), ghb(:,:), ghg(:,:), &
807  rnr(:,:), rnb(:,:), rng(:,:), &
808  rngrd(:,:) )
809 
810  !$acc end data
811 
812 #ifdef DEBUG_KUSAKA01
813  cnt_num1 = cnt_num1 + cn1
814  cnt_num2 = cnt_num2 + cn2
815  cnt_itr1 = cnt_itr1 + ci1
816  cnt_itr2 = cnt_itr2 + ci2
817  max_itr1 = max(max_itr1, mi1)
818  max_itr2 = max(max_itr2, mi2)
819 #endif
820 
821  return
822  end subroutine urban_dyn_kusaka01
823 
824  !-----------------------------------------------------------------------------
825 !OCL SERIAL
826  subroutine slc_main( &
827  UKA, UKS, UKE, &
828  TRL, & ! (inout)
829  tbl, & ! (inout)
830  tgl, & ! (inout)
831  tr, & ! (inout)
832  tb, & ! (inout)
833  tg, & ! (inout)
834  tc, & ! (inout)
835  qc, & ! (inout)
836  uc, & ! (inout)
837  rainr, & ! (inout)
838  rainb, & ! (inout)
839  raing, & ! (inout)
840  roff, & ! (out)
841  albd_lw_grid, & ! (out)
842  albd_sw_grid, & ! (out)
843  shr, & ! (out)
844  shb, & ! (out)
845  shg, & ! (out)
846  lhr, & ! (out)
847  lhb, & ! (out)
848  lhg, & ! (out)
849  ghr, & ! (out)
850  ghb, & ! (out)
851  ghg, & ! (out)
852  rnr, & ! (out)
853  rnb, & ! (out)
854  rng, & ! (out)
855  rts, & ! (out)
856  rn, & ! (out)
857  sh, & ! (out)
858  lh, & ! (out)
859  ghflx, & ! (out)
860  mflx, & ! (out)
861  ustar, & ! (out)
862  tstar, & ! (out)
863  qstar, & ! (out)
864  u10, & ! (out)
865  v10, & ! (out)
866  t2, & ! (out)
867  q2, & ! (out)
868  converged, & ! (out)
869  lsolar, & ! (in)
870  prsa, & ! (in)
871  prss, & ! (in)
872  ta, & ! (in)
873  rhos, & ! (in)
874  qa, & ! (in)
875  ua, & ! (in)
876  u1, & ! (in)
877  v1, & ! (in)
878  lhv, & ! (in)
879  za, & ! (in)
880  ssg, & ! (in)
881  llg, & ! (in)
882  rain, & ! (in)
883  eflx, & ! (in)
884  rhoo, & ! (in)
885  z0c, & ! (in)
886  z0hc, & ! (in)
887  zdc, & ! (in)
888  dzr, dzb, dzg, & ! (in)
889  dt, & ! (in)
890  i, j, & ! (in)
891 #ifdef DEBUG_KUSAKA01
892  cnt_num1, cnt_num2, & ! (inout)
893  cnt_itr1, cnt_itr2, & ! (inout)
894  max_itr1, max_itr2, & ! (inout)
895 #endif
896  trlp, tblp, tglp, &
897  a, b, c, d, p, q )
898  !$acc routine seq
899  use scale_prc, only: &
900  prc_myrank, &
901  prc_abort
902  use scale_const, only: &
903  eps => const_eps, & ! small number (machine epsilon)
904  karman => const_karman, & ! Kalman constant [-]
905  cpdry => const_cpdry, & ! Heat capacity of dry air [J/K/kg]
906  grav => const_grav, & ! gravitational constant [m/s2]
907  rdry => const_rdry, & ! specific gas constant (dry) [J/kg/K]
908  rvap => const_rvap, & ! gas constant (water vapor) [J/kg/K]
909  stb => const_stb, & ! stefan-Boltzman constant [MKS unit]
910  tem00 => const_tem00, & ! temperature reference (0 degree C) [K]
911  pre00 => const_pre00 ! pressure reference [Pa]
912  use scale_atmos_hydrometeor, only: &
913  hydrometeor_lhv => atmos_hydrometeor_lhv
914  use scale_atmos_saturation, only: &
915  qsat => atmos_saturation_dens2qsat_liq, &
916  dqs_dtem => atmos_saturation_dqs_dtem_dens_liq
917  use scale_bulkflux, only: &
918  bulkflux_diagnose_surface
919  implicit none
920 
921  integer, intent(in) :: uka, uks, uke
922 
923  !-- In/Out variables from/to Coupler to/from Urban
924  real(rp), intent(inout) :: trl(uks:uke) ! layer temperature [K]
925  real(rp), intent(inout) :: tbl(uks:uke) ! layer temperature [K]
926  real(rp), intent(inout) :: tgl(uks:uke) ! layer temperature [K]
927  real(rp), intent(inout) :: tr ! roof temperature [K]
928  real(rp), intent(inout) :: tb ! building wall temperature [K]
929  real(rp), intent(inout) :: tg ! road temperature [K]
930  real(rp), intent(inout) :: tc ! urban-canopy air temperature [K]
931  real(rp), intent(inout) :: qc ! urban-canopy air specific humidity [kg/kg]
932  real(rp), intent(inout) :: uc ! diagnostic canopy wind [m/s]
933  real(rp), intent(inout) :: rainr ! rain amount in storage on roof [kg/m2]
934  real(rp), intent(inout) :: rainb ! rain amount in storage on building [kg/m2]
935  real(rp), intent(inout) :: raing ! rain amount in storage on road [kg/m2]
936  real(rp), intent(out) :: roff ! runoff from urban [kg/m2/s]
937 
938  !-- Output variables from Urban to Coupler
939  real(rp), intent(out) :: albd_sw_grid ! grid mean of surface albedo for SW
940  real(rp), intent(out) :: albd_lw_grid ! grid mean of surface albedo for LW ( 1-emiss )
941  real(rp), intent(out) :: shr, shb, shg
942  real(rp), intent(out) :: lhr, lhb, lhg
943  real(rp), intent(out) :: ghr, ghb, ghg
944  real(rp), intent(out) :: rnr, rnb, rng
945  real(rp), intent(out) :: rts ! radiative surface temperature [K]
946  real(rp), intent(out) :: rn ! net radition [W/m/m]
947  real(rp), intent(out) :: sh ! sensible heat flux [W/m/m]
948  real(rp), intent(out) :: lh ! latent heat flux [W/m/m]
949  real(rp), intent(out) :: ghflx ! heat flux into the ground [W/m/m]
950  real(rp), intent(out) :: mflx ! momentum flux [kg/m/s2]
951  real(rp), intent(out) :: ustar ! friction velocity [m/s]
952  real(rp), intent(out) :: tstar ! temperature scale [K]
953  real(rp), intent(out) :: qstar ! humidity scale [kg/kg]
954  real(rp), intent(out) :: u10 ! U wind at 10m [m/s]
955  real(rp), intent(out) :: v10 ! V wind at 10m [m/s]
956  real(rp), intent(out) :: t2 ! air temperature at 2m [K]
957  real(rp), intent(out) :: q2 ! specific humidity at 2m [kg/kg]
958  logical, intent(out) :: converged
959 
960  !-- configuration variables
961  logical , intent(in) :: lsolar ! logical [true=both, false=SSG only]
962 
963  !-- Input variables from Coupler to Urban
964  real(rp), intent(in) :: prsa ! Pressure at 1st atmospheric layer [Pa]
965  real(rp), intent(in) :: prss ! Surface Pressure [Pa]
966  real(rp), intent(in) :: ta ! temp at 1st atmospheric level [K]
967  real(rp), intent(in) :: rhos ! surface density [kg/m^3]
968  real(rp), intent(in) :: qa ! specific humidity at 1st atmospheric level [kg/kg]
969  real(rp), intent(in) :: ua ! wind speed at 1st atmospheric level [m/s]
970  real(rp), intent(in) :: u1 ! u at 1st atmospheric level [m/s]
971  real(rp), intent(in) :: v1 ! v at 1st atmospheric level [m/s]
972  real(rp), intent(in) :: lhv ! latent heat of vaporization [J/kg]
973  real(rp), intent(in) :: za ! height of 1st atmospheric level [m]
974  real(rp), intent(in) :: ssg(2) ! downward total short wave radiation [W/m/m]
975  real(rp), intent(in) :: llg(2) ! downward long wave radiation [W/m/m]
976  real(rp), intent(in) :: rain ! water flux [kg/m2/s]
977  real(rp), intent(in) :: eflx ! internal energy flux [J/m2/s]
978  real(rp), intent(in) :: rhoo ! air density [kg/m^3]
979  real(rp), intent(in) :: z0c ! Roughness length above canyon for momentum [m]
980  real(rp), intent(in) :: z0hc ! Roughness length above canyon for heat [m]
981  real(rp), intent(in) :: zdc ! Displacement height [m]
982  real(rp), intent(in) :: dzr(uka)
983  real(rp), intent(in) :: dzb(uka)
984  real(rp), intent(in) :: dzg(uka)
985  real(dp), intent(in) :: dt
986 
987  integer, intent(in) :: i, j
988 
989 #ifdef DEBUG_KUSAKA01
990  integer, intent(inout) :: cnt_num1, cnt_num2
991  integer, intent(inout) :: cnt_itr1, cnt_itr2
992  integer, intent(inout) :: max_itr1, max_itr2
993 #endif
994 
995  ! work
996  real(rp), intent(out) :: trlp(uka) ! Layer temperature at previous step [K]
997  real(rp), intent(out) :: tblp(uka) ! Layer temperature at previous step [K]
998  real(rp), intent(out) :: tglp(uka) ! Layer temperature at previous step [K]
999  real(rp), intent(out) :: a(uka)
1000  real(rp), intent(out) :: b(uka)
1001  real(rp), intent(out) :: c(uka)
1002  real(rp), intent(out) :: d(uka)
1003  real(rp), intent(out) :: p(uka)
1004  real(rp), intent(out) :: q(uka)
1005 
1006  !-- parameters
1007 ! real(RP), parameter :: SRATIO = 0.75_RP ! ratio between direct/total solar [-]
1008 ! real(RP), parameter :: TFa = 0.5_RP ! factor a in Tomita (2009)
1009 ! real(RP), parameter :: TFb = 1.1_RP ! factor b in Tomita (2009)
1010  real(rp), parameter :: redf_min = 1.0e-2_rp ! minimum reduced factor
1011  real(rp), parameter :: redf_max = 1.0_rp ! maximum reduced factor
1012 ! real(RP), parameter :: CAP_water = 4.185E6_RP ! Heat capacity of water (15 deg) [J m-3 K]
1013 ! real(RP), parameter :: AKS_water = 0.59_RP ! Thermal conductivity of water [W m-1 K]
1014 
1015  real(rp), parameter :: rain_rate_r = 1.0_rp
1016  real(rp), parameter :: rain_rate_b = 0.1_rp
1017  real(rp), parameter :: rain_rate_g = 0.9_rp
1018 
1019  integer, parameter :: itr_max = 100
1020 
1021  !-- Local variables
1022 ! logical :: SHADOW = .false.
1023  ! true = consider svf and shadow effects,
1024  ! false = consider svf effect only
1025 
1026  real(rp) :: w, vfgs, vfgw, vfwg, vfws, vfww
1027  real(rp) :: rflux_sw, rflux_lw
1028 
1029  real(rp) :: trp ! TRP: at previous time step [K]
1030  real(rp) :: tbp ! TBP: at previous time step [K]
1031  real(rp) :: tgp ! TGP: at previous time step [K]
1032  real(rp) :: tcp ! TCP: at previous time step [K]
1033  real(rp) :: qcp ! QCP: at previous time step [kg/kg]
1034  !
1035  real(rp) :: rainrp ! at previous step, rain amount in storage on roof [kg/m2]
1036  real(rp) :: rainbp ! at previous step, rain amount in storage on building [kg/m2]
1037  real(rp) :: raingp ! at previous step, rain amount in storage on road [kg/m2]
1038 
1039  real(rp) :: raint
1040  real(rp) :: roffr, roffb, roffg ! runoff [kg/m2]
1041 
1042  real(rp) :: lup, ldn, rup
1043  real(rp) :: sup, sdn
1044 
1045  real(rp) :: lnet, snet, flxuv
1046  real(rp) :: sw ! shortwave radition [W/m/m]
1047  real(rp) :: lw ! longwave radition [W/m/m]
1048  real(rp) :: psim,psim2,psim10 ! similality stability shear function for momentum
1049  real(rp) :: psih,psih2,psih10 ! similality stability shear function for heat
1050 
1051  ! for shadow effect model
1052  ! real(RP) :: HOUI1, HOUI2, HOUI3, HOUI4, HOUI5, HOUI6, HOUI7, HOUI8
1053  ! real(RP) :: SLX, SLX1, SLX2, SLX3, SLX4, SLX5, SLX6, SLX7, SLX8
1054  ! real(RP) :: THEATAZ ! Solar Zenith Angle [rad]
1055  ! real(RP) :: THEATAS ! = PI/2. - THETAZ
1056  ! real(RP) :: FAI ! Latitude [rad]
1057 
1058  real(rp) :: sr, sb, sg, rr, rb, rg
1059  real(rp) :: sr1, sb1, sb2, sg1, sg2
1060  real(rp) :: rb1, rb2, rg1, rg2
1061  real(rp) :: hr, evpr, eler, g0r, betr
1062  real(rp) :: hb, evpb, eleb, g0b, betb
1063  real(rp) :: hg, evpg, eleg, g0g, betg
1064  real(rp) :: evprp, evpbp, evpgp
1065 
1066  real(rp) :: z
1067  real(rp) :: qs0r, qs0b, qs0g
1068 
1069  real(rp) :: ribr, bhr, cdr
1070  real(rp) :: ribc, bhc, cdc
1071  real(rp) :: alphab, alphag, alphac
1072  real(rp) :: chr, chb, chg, chc
1073  real(rp) :: tc1, tc2, qc1, qc2
1074 ! real(RP) :: CAPL1, AKSL1
1075 
1076  real(rp) :: xxx, xxx2, xxx10
1077  real(rp) :: xxxr ! Monin-Obkhov length for roof [-]
1078  real(rp) :: xxxc ! Monin-Obkhov length for canopy [-]
1079 
1080  real(rp) :: tha,thc,ths,ths1,ths2
1081  real(rp) :: rovcp
1082  real(rp) :: exn ! exner function at the surface
1083 
1084  real(rp) :: fracu10, fract2
1085 
1086  real(rp) :: dt_rp
1087 
1088  integer :: iteration
1089  real(rp) :: dts_max_onestep ! DTS_MAX * dt
1090  real(rp) :: resi1, resi2, resi3 ! residual
1091  real(rp) :: fact1, fact2, fact3
1092  real(rp) :: threshold
1093 
1094  ! for Newton method
1095  real(rp) :: dtr, dtb, dtg, dac
1096  real(rp) :: dtrp, dtbp, dtgp, dacp
1097  real(rp) :: dr1dtr
1098  real(rp) :: dr1dtb, dr1dtg, dr1dac, dr2dtb, dr2dtg, dr2dac, dr3dtb, dr3dtg, dr3dac
1099  real(rp) :: dtrdg0r, dtbdg0b, dtgdg0g
1100  real(rp) :: dg0rdtr, dg0bdtb, dg0bdtg, dg0bdac, dg0gdtb, dg0gdtg, dg0gdac
1101  real(rp) :: dhr, dhbdtb, dhbdtg, dhbdac, dhgdtb, dhgdtg, dhgdac
1102  real(rp) :: deler, delebdtb, delebdtg, delebdac, delegdtb, delegdtg, delegdac
1103  real(rp) :: drr
1104  real(rp) :: dbetr, dbetb, dbetg, betrp, betbp, betgp
1105  real(rp) :: drb1dtb, drb1dtg, drb2dtb, drb2dtg
1106  real(rp) :: drg1dtb, drg1dtg, drg2dtb, drg2dtg
1107  real(rp) :: dqcdtb, dqcdtg, dqcdac
1108  real(rp) :: dthcdtb, dthcdtg, dthcdac
1109  real(rp) :: dalphacdtb, dalphacdtg, dalphacdac
1110  real(rp) :: dchcdtb, dchcdtg, dchcdac
1111  real(rp) :: chc_tb, chc_tg, chc_ac
1112  real(rp) :: dqs0r, dqs0b, dqs0g
1113  real(rp) :: tdiff, adiff
1114  real(rp) :: xxxtmp
1115  real(rp) :: alphacp
1116  real(rp) :: rdet
1117  real(rp) :: b1, b2
1118  real(rp) :: tmp
1119 
1120  integer :: k
1121 
1122  !-----------------------------------------------------------
1123  ! Set parameters
1124  !-----------------------------------------------------------
1125 
1126  dt_rp = dt
1127 
1128  threshold = sqrt(eps)
1129 
1130  rovcp = rdry / cpdry
1131  tha = ta * ( pre00 / prsa )**rovcp
1132 
1133  !--- Renew surface and layer temperatures
1134 
1135  trp = tr
1136  tbp = tb
1137  tgp = tg
1138  tcp = tc
1139  qcp = qc
1140  !
1141  do k = uks, uke
1142  trlp(k) = trl(k)
1143  tblp(k) = tbl(k)
1144  tglp(k) = tgl(k)
1145  end do
1146  !
1147 
1148 
1149  !--- limiter for surface temp change
1150  dts_max_onestep = dts_max * dt_rp
1151 
1152  ! "2.0m" has no special meaning, but it is related with BB formulation from Inoue (1963).
1153  ! Please see subroutine "canopy_wind".
1154  ! The canopy model is modeled under an assumption that urban canopy lies
1155  ! below the lowest level of atmospheric model.
1156  if ( zdc + z0c + 2.0_rp >= za ) then
1157  log_error("URBAN_DYN_kusaka01_SLC_main",*) 'ZDC + Z0C + 2m must be less than the 1st level! STOP.'
1158  converged = .false.
1159 #ifdef _OPENACC
1160  return
1161 #else
1162  call prc_abort
1163 #endif
1164  endif
1165 
1166  w = 2.0_rp * 1.0_rp * hgt
1167 
1168  rflux_sw = ssg(1) + ssg(2) ! downward shortwave radiation [W/m2]
1169  rflux_lw = llg(1) + llg(2) ! downward longwave radiation [W/m2]
1170 
1171  !--- calculate canopy wind
1172 
1173  call canopy_wind(za, ua, z0c, zdc, uc)
1174 
1175  !-----------------------------------------------------------
1176  ! calculate water content (temporally)
1177  !-----------------------------------------------------------
1178 
1179  raint = rain * dt_rp ! [kg/m2/s -> kg/m2]
1180 
1181  rainr = rainr + raint * rain_rate_r
1182  rainb = rainb + raint * rain_rate_b
1183  raing = raing + raint * rain_rate_g
1184 
1185  vfgs = svf
1186  vfgw = 1.0_rp - svf
1187  vfwg = ( 1.0_rp - svf ) * ( 1.0_rp - r ) / w
1188  vfws = vfwg
1189  vfww = 1.0_rp - 2.0_rp * vfwg
1190 
1191  ! save the initial value
1192  rainrp = rainr
1193  rainbp = rainb
1194  raingp = raing
1195 
1196  !-----------------------------------------------------------
1197  ! Radiation : Net Short Wave Radiation at roof/wall/road
1198  !-----------------------------------------------------------
1199 
1200  if( rflux_sw > 0.0_rp ) then ! SSG is downward short
1201 
1202  ! currently we use no shadow effect model
1203  !! IF(.NOT.SHADOW) THEN ! no shadow effects model
1204 
1205  sr1 = rflux_sw * ( 1.0_rp - albr )
1206  sg1 = rflux_sw * vfgs * ( 1.0_rp - albg )
1207  sb1 = rflux_sw * vfws * ( 1.0_rp - albb )
1208  sg2 = sb1 * albb / ( 1.0_rp - albb ) * vfgw * ( 1.0_rp - albg )
1209  sb2 = sg1 * albg / ( 1.0_rp - albg ) * vfwg * ( 1.0_rp - albb )
1210 
1211  sr = sr1
1212  sg = sg1 + sg2
1213  sb = sb1 + sb2
1214  snet = r * sr + w * sb + rw * sg
1215 
1216  else
1217 
1218  sr = 0.0_rp
1219  sg = 0.0_rp
1220  sb = 0.0_rp
1221  snet = 0.0_rp
1222 
1223  end if
1224 
1225 
1226  exn = ( prss / pre00 )**rovcp ! exner function
1227 
1228  !-----------------------------------------------------------
1229  ! Energy balance on roof/wall/road surface
1230  !-----------------------------------------------------------
1231 
1232  !--------------------------------------------------
1233  ! Roof
1234  !--------------------------------------------------
1235 
1236  ! new scheme
1237 
1238  z = za - zdc
1239  bhr = log(z0r/z0hr) / 0.4_rp
1240 
1241  b1 = capr * dzr(1) / dt_rp + 2.0_rp * aksr / ( dzr(1)+dzr(2) )
1242  b2 = capr * dzr(2) / dt_rp + 2.0_rp * aksr / ( dzr(1)+dzr(2) ) + 2.0_rp * aksr / ( dzr(2)+dzr(3) )
1243  dtrdg0r = 1.0_rp / ( b1 - ( 2.0_rp * aksr / ( dzr(1)+dzr(2) ) )**2 / b2 )
1244  ! consider only change at the layers of k<=2
1245  dtrdg0r = dtrdg0r * 0.5_rp
1246 
1247  evpr = 0.0
1248  evprp = 0.0
1249  betrp = 0.0_rp
1250  xxxr = 0.0_rp
1251  dtrp = 1.0e10_rp
1252  fact1 = 1.0_rp
1253  !$acc loop seq
1254  do iteration = 1, itr_max
1255 
1256  ths = tr / exn ! potential temp
1257 
1258  ribr = ( grav * 2.0_rp / (tha+ths) ) * (tha-ths) * (z+z0r) / (ua*ua+eps)
1259  call mos(xxxr,chr,cdr,bhr,ribr,z,z0r,ua,tha,ths,rhoo,i,j)
1260  ! ignore differential of CHR for Newton method
1261 
1262  call dqs_dtem( tr, rhos, & ! [IN]
1263  dqs0r, qs0r ) ! [OUT]
1264 
1265  rainr = max( rainrp - evpr * dt_rp, 0.0_rp )
1266  call cal_beta(betr, betr_const, rainr, strgr)
1267  dbetr = ( betr - betrp ) / sign( max(abs(dtrp), 1e-10_rp), dtrp )
1268  betrp = betr
1269 
1270  rr = epsr * ( rflux_lw - stb * (tr**4) )
1271  drr = - epsr * stb * 4.0_rp * tr**3
1272 
1273  !HR = RHOO * CPdry * CHR * UA * (TR-TA)
1274  hr = rhoo * cpdry * chr * ua * ( ths - tha ) * exn
1275  dhr = rhoo * cpdry * chr * ua
1276 
1277  evpr = min( rhoo * chr * ua * betr * (qs0r-qa), rainr / dt_rp )
1278  eler = evpr * lhv
1279  deler = rhoo * chr * ua * ( betr * dqs0r + dbetr * qs0r ) * lhv
1280 
1281  g0r = sr + rr - hr - eler
1282  dg0rdtr = drr - dhr - deler
1283 
1284  !--- calculate temperature in roof
1285  ! if ( STRGR /= 0.0_RP ) then
1286  ! CAPL1 = CAP_water * (RAINR / (DZR(1) + RAINR)) + CAPR * (DZR(1) / (DZR(1) + RAINR))
1287  ! AKSL1 = AKS_water * (RAINR / (DZR(1) + RAINR)) + AKSR * (DZR(1) / (DZR(1) + RAINR))
1288  ! else
1289  ! CAPL1 = CAPR
1290  ! AKSL1 = AKSR
1291  ! endif
1292  !! 1st layer's cap, aks are replaced.
1293  !! call multi_layer2(UKE,BOUND,G0R,CAPR,AKSR,TRL,DZR,dt_RP,TRLEND,CAPL1,AKSL1)
1294 
1295  do k = uks, uke
1296  trl(k) = trlp(k)
1297  end do
1298  call multi_layer(uke,bound, &
1299  a, b, c, d, p, q, &
1300  g0r,capr,aksr,trl,dzr,dt_rp,trlend)
1301  resi1 = trl(1) - tr
1302 
1303  if( abs(resi1) < threshold ) then
1304  tr = trl(1)
1305  tr = max( trp - dts_max_onestep, min( trp + dts_max_onestep, tr ) )
1306  exit
1307  endif
1308 
1309  ! Newton method
1310  dr1dtr = dtrdg0r * dg0rdtr - 1.0_rp
1311 
1312  dtr = - resi1 / dr1dtr
1313  dtr = sign( min(abs(dtr), 5.0_rp), dtr )
1314 
1315  if ( iteration > 50 ) then
1316  if ( dtr*dtrp < 0.0_rp ) then
1317  fact1 = max( fact1 * 0.5_rp, 0.01_rp )
1318  else
1319  fact1 = min( fact1 * 1.5_rp, 1.1_rp )
1320  end if
1321  end if
1322  tmp = tr
1323  tr = max( tr + dtr * fact1, 100.0_rp )
1324  dtrp = tr - tmp
1325 
1326  evpr = ( evprp + evpr ) * 0.5_rp
1327  evprp = evpr
1328  enddo
1329 
1330 #ifdef DEBUG_KUSAKA01
1331  cnt_num1 = cnt_num1 + 1
1332  cnt_itr1 = cnt_itr1 + iteration - 1
1333  max_itr1 = max(max_itr1, iteration - 1)
1334 #endif
1335 
1336  ! output for debug
1337  if ( iteration > itr_max .and. debug ) then
1338  log_warn("URBAN_DYN_kusaka01_SLC_main",*) 'iteration for TR was not converged',prc_myrank,i,j
1339  log_warn_cont(*) '---------------------------------------------------------------------------------'
1340  log_warn_cont(*) 'DEBUG Message --- Residual [K] :', resi1
1341  log_warn_cont(*) 'DEBUG Message --- TRP : Initial TR [K] :', trp
1342 #ifdef _OPENACC
1343  log_warn_cont(*) 'DEBUG Message --- TRLP: Initial TRL [K] :', trlp(uks)
1344 #else
1345  log_warn_cont(*) 'DEBUG Message --- TRLP: Initial TRL [K] :', trlp(:)
1346 #endif
1347  log_warn_cont(*) 'DEBUG Message --- rflux_SW : Shortwave radiation [W/m2] :', rflux_sw
1348  log_warn_cont(*) 'DEBUG Message --- rflux_LW : Longwave radiation [W/m2] :', rflux_lw
1349  log_warn_cont(*) 'DEBUG Message --- PRSS: Surface pressure [Pa] :', prss
1350  log_warn_cont(*) 'DEBUG Message --- PRSA: Pressure at 1st atmos layer [m] :', prsa
1351  log_warn_cont(*) 'DEBUG Message --- RHOO: Air density [kg/m3] :', rhoo
1352  log_warn_cont(*) 'DEBUG Message --- RHOS: Surface density [kg/m3] :', rhos
1353  log_warn_cont(*) 'DEBUG Message --- RAINRP: Initial RAINR [kg/m2] :', rainrp
1354  log_warn_cont(*) 'DEBUG Message --- ZA : Height at 1st atmos layer [m] :', za
1355  log_warn_cont(*) 'DEBUG Message --- TA : Temperature at 1st atmos layer [K] :', ta
1356  log_warn_cont(*) 'DEBUG Message --- UA : Wind speed at 1st atmos layer [m/s] :', ua
1357  log_warn_cont(*) 'DEBUG Message --- QA : Specific humidity at 1st atmos layer [kg/kg] :', qa
1358 #ifdef _OPENACC
1359  log_warn_cont(*) 'DEBUG Message --- DZR : Depth of surface layer [m] :', dzr(1)
1360 #else
1361  log_warn_cont(*) 'DEBUG Message --- DZR : Depth of surface layer [m] :', dzr(:)
1362 #endif
1363  log_warn_cont(*) 'DEBUG Message --- R, W, RW : Normalized height and road width [-] :', r, w,rw
1364  log_warn_cont(*) 'DEBUG Message --- SVF : Sky View Factors [-] :', svf
1365  log_warn_cont(*) 'DEBUG Message --- BETR: Evaporation efficiency [-] :', betr
1366  log_warn_cont(*) 'DEBUG Message --- EPSR: Surface emissivity of roof [-] :', epsr
1367  log_warn_cont(*) 'DEBUG Message --- CAPR: Heat capacity of roof [J m-3 K] :', capr
1368  log_warn_cont(*) 'DEBUG Message --- AKSR: Thermal conductivity of roof [W m-1 K] :', aksr
1369  log_warn_cont(*) 'DEBUG Message --- QS0R: Surface specific humidity [kg/kg] :', qs0r
1370  log_warn_cont(*) 'DEBUG Message --- ZDC : Desplacement height of canopy [m] :', zdc
1371  log_warn_cont(*) 'DEBUG Message --- Z0R : Momentum roughness length of roof [m] :', z0r
1372  log_warn_cont(*) 'DEBUG Message --- Z0HR: Thermal roughness length of roof [m] :', z0hr
1373  log_warn_cont(*) '---------------------------------------------------------------------------------'
1374  endif
1375 
1376  !--- update only fluxes ----
1377  ths = tr / exn
1378  ribr = ( grav * 2.0_rp / (tha+ths) ) * (tha-ths) * (z+z0r) / (ua*ua+eps)
1379  call mos(xxxr,chr,cdr,bhr,ribr,z,z0r,ua,tha,ths,rhoo,i,j)
1380 
1381  call qsat( tr, rhos, & ! [IN]
1382  qs0r ) ! [OUT]
1383 
1384  call cal_beta(betr, betr_const, rainr, strgr)
1385 
1386  rr = epsr * ( rflux_lw - stb * (tr**4) )
1387  hr = rhoo * cpdry * chr * ua * (ths-tha) * exn
1388  evpr = min( rhoo * chr * ua * betr * (qs0r-qa), rainr / dt_rp )
1389  eler = evpr * lhv
1390 
1391 ! G0R = SR + RR - HR - ELER + EFLX
1392  g0r = sr + rr - hr - eler
1393  rainr = max( rainrp - evpr * dt_rp, 0.0_rp )
1394 
1395  do k = uks, uke
1396  trl(k) = trlp(k)
1397  end do
1398  call multi_layer(uke,bound, &
1399  a, b, c, d, p, q, &
1400  g0r,capr,aksr,trl,dzr,dt_rp,trlend)
1401  resi1 = trl(1) - tr
1402  tr = trl(1)
1403 
1404  if ( abs(resi1) > dts_max_onestep ) then
1405  if ( abs(resi1) > dts_max_onestep*10.0_rp ) then
1406  log_error("URBAN_DYN_Kusaka01_main",*) 'tendency of TR exceeded a limit! STOP.'
1407  log_error_cont(*) 'previous TR and updated TR(TRL(1)) is ',tr-resi1, tr
1408  converged = .false.
1409 #ifdef _OPENACC
1410  return
1411 #else
1412  call prc_abort
1413 #endif
1414  endif
1415  log_warn("URBAN_DYN_Kusaka01_main",*) 'tendency of TR exceeded a limit'
1416  log_warn_cont(*) 'previous TR and updated TR(TRL(1)) is ', tr-resi1, tr
1417  endif
1418 
1419  !--------------------------------------------------
1420  ! Wall and Road
1421  !--------------------------------------------------
1422 
1423  ! new scheme
1424 
1425  ! empirical form
1426  alphab = 6.15_rp + 4.18_rp * uc
1427  if( uc > 5.0_rp ) alphab = 7.51_rp * (uc**0.78_rp )
1428  alphag = 6.15_rp + 4.18_rp * uc
1429  if( uc > 5.0_rp ) alphag = 7.51_rp * (uc**0.78_rp )
1430  chb = alphab / rhoo / cpdry / uc
1431  chg = alphag / rhoo / cpdry / uc
1432 
1433  z = za - zdc
1434  thc = tc / exn
1435  bhc = log(z0c/z0hc) / 0.4_rp
1436 
1437  b1 = capb * dzb(1) / dt_rp + 2.0_rp * aksb / ( dzb(1)+dzb(2) )
1438  b2 = capb * dzb(2) / dt_rp + 2.0_rp * aksb / ( dzb(1)+dzb(2) ) + 2.0_rp * aksb / ( dzb(2)+dzb(3) )
1439  dtbdg0b = 1.0_rp / ( b1 - ( 2.0_rp * aksb / ( dzb(1)+dzb(2) ) )**2 / b2 )
1440  ! consider only change at the layers of k<=2
1441  dtbdg0b = dtbdg0b * 0.5_rp
1442 
1443  b1 = capg * dzg(1) / dt_rp + 2.0_rp * aksg / ( dzg(1)+dzg(2) )
1444  b2 = capg * dzg(2) / dt_rp + 2.0_rp * aksg / ( dzg(1)+dzg(2) ) + 2.0_rp * aksg / ( dzg(2)+dzg(3) )
1445  dtgdg0g = 1.0_rp / ( b1 - ( 2.0_rp * aksg / ( dzg(1)+dzg(2) ) )**2 / b2 )
1446  ! consider only change at the layers of k<=2
1447  dtgdg0g = dtgdg0g * 0.5_rp
1448 
1449  evpb = 0.0_rp
1450  evpbp = 0.0_rp
1451  evpg = 0.0_rp
1452  evpgp = 0.0_rp
1453  betbp = 0.0_rp
1454  betgp = 0.0_rp
1455  xxxc = 0.0_rp
1456  alphacp = ( alphab + alphag ) * 0.5_rp
1457  alphac = alphacp
1458  dtbp = 1.0e10_rp
1459  dtgp = 1.0e10_rp
1460  dacp = 1.0_rp
1461  fact1 = 1.0_rp
1462  fact2 = 1.0_rp
1463  fact3 = 1.0_rp
1464  !$acc loop seq
1465  do iteration = 1, itr_max
1466 
1467  ths1 = tb / exn
1468  ths2 = tg / exn
1469 
1470  if ( iteration > 1 ) then
1471 
1472  tdiff = tb * sqrt(eps) * 2.0_rp
1473  adiff = sign( alphac * sqrt(eps) * 2.0_rp, dac )
1474 
1475  ! TB = TB + Tdiff
1476  tc1 = rw*alphac + rw*alphag + w*alphab
1477  tc2 = rw*alphac*tha + w*alphab*(tb+tdiff)/exn + rw*alphag*ths2
1478  thc = tc2 / tc1
1479  ribc = ( grav * 2.0_rp / (tha+thc) ) * (tha-thc) * (z+z0c) / (ua*ua+eps)
1480  xxxtmp = xxxc
1481  ribc = ( grav * 2.0_rp / (tha+thc) ) * (tha-thc) * (z+z0c) / (ua*ua+eps)
1482  call mos(xxxtmp,chc_tb,cdc,bhc,ribc,z,z0c,ua,tha,thc,rhoo,i,j)
1483 
1484  ! TG = TG + Tdiff
1485  tc1 = rw*alphac + rw*alphag + w*alphab
1486  tc2 = rw*alphac*tha + w*alphab*ths1 + rw*alphag*(tg+tdiff)/exn
1487  thc = tc2 / tc1
1488  ribc = ( grav * 2.0_rp / (tha+thc) ) * (tha-thc) * (z+z0c) / (ua*ua+eps)
1489  xxxtmp = xxxc
1490  call mos(xxxtmp,chc_tg,cdc,bhc,ribc,z,z0c,ua,tha,thc,rhoo,i,j)
1491 
1492  ! ALPHAC = ALPHAC + Adiff
1493  tc1 = rw*(alphac+adiff) + rw*alphag + w*alphab
1494  tc2 = rw*(alphac+adiff)*tha + w*alphab*ths1 + rw*alphag*ths2
1495  thc = tc2 / tc1
1496  ribc = ( grav * 2.0_rp / (tha+thc) ) * (tha-thc) * (z+z0c) / (ua*ua+eps)
1497  xxxtmp = xxxc
1498  call mos(xxxtmp,chc_ac,cdc,bhc,ribc,z,z0c,ua,tha,thc,rhoo,i,j)
1499 
1500  end if
1501 
1502  tc1 = rw*alphac + rw*alphag + w*alphab
1503  !TC2 = RW*ALPHAC*THA + RW*ALPHAG*TG + W*ALPHAB*TB
1504  tc2 = rw*alphac*tha + w*alphab*ths1 + rw*alphag*ths2
1505  thc = tc2 / tc1
1506  ribc = ( grav * 2.0_rp / (tha+thc) ) * (tha-thc) * (z+z0c) / (ua*ua+eps)
1507  call mos(xxxc,chc,cdc,bhc,ribc,z,z0c,ua,tha,thc,rhoo,i,j)
1508 
1509  if ( iteration > 1 ) then
1510  dchcdtb = ( chc_tb - chc ) / tdiff
1511  dchcdtg = ( chc_tg - chc ) / tdiff
1512  dchcdac = ( chc_ac - chc ) / adiff
1513  else
1514  dchcdtb = 0.0_rp
1515  dchcdtg = 0.0_rp
1516  dchcdac = 0.0_rp
1517  end if
1518 
1519  alphac = chc * rhoo * cpdry * ua
1520  dalphacdtb = dchcdtb * rhoo * cpdry * ua
1521  dalphacdtg = dchcdtg * rhoo * cpdry * ua
1522  dalphacdac = dchcdac * rhoo * cpdry * ua
1523 
1524  call dqs_dtem( tb, rhos, dqs0b, qs0b )
1525  call dqs_dtem( tg, rhos, dqs0g, qs0g )
1526 
1527  rainb = max( ( rainbp - evpb * dt_rp ), 0.0_rp )
1528  raing = max( ( raingp - evpg * dt_rp ), 0.0_rp )
1529  call cal_beta(betb, betb_const, rainb, strgb)
1530  call cal_beta(betg, betg_const, raing, strgg)
1531  dbetb = ( betb - betbp ) / sign( max(abs(dtbp),1e-10_rp), dtbp )
1532  dbetg = ( betg - betgp ) / sign( max(abs(dtgp),1e-10_rp), dtgp )
1533  betbp = betb
1534  betgp = betg
1535 
1536 
1537  tc1 = rw*alphac + rw*alphag + w*alphab
1538  !TC2 = RW*ALPHAC*TA + RW*ALPHAG*TG + W*ALPHAB*TB
1539  tc2 = rw*alphac*tha + w*alphab*ths1 + rw*alphag*ths2
1540  thc = tc2 / tc1
1541  dthcdtb = ( ( rw*dalphacdtb*tha + w*alphab/exn ) * tc1 - rw*dalphacdtb * tc2 ) / tc1**2
1542  dthcdtg = ( ( rw*dalphacdtg*tha + rw*alphag/exn ) * tc1 - rw*dalphacdtg * tc2 ) / tc1**2
1543  dthcdac = ( ( rw* tha ) * tc1 - rw * tc2 ) / tc1**2
1544 
1545  qc1 = rw*(chc*ua) + rw*(chg*betg*uc) + w*(chb*betb*uc)
1546  qc2 = rw*(chc*ua)*qa + rw*(chg*betg*uc)*qs0g + w*(chb*betb*uc)*qs0b
1547  qc = max( qc2 / ( qc1 + eps ), 0.0_rp )
1548  dqcdtb = ( ( rw*(dchcdtb*ua)*qa + w*(chb*uc)*(betb*dqs0b+dbetb*qs0b) ) * qc1 &
1549  - ( rw*(dchcdtb*ua) + w*(chb*uc*dbetb) ) * qc2 ) / ( qc1**2 + eps )
1550  dqcdtg = ( ( rw*(dchcdtg*ua)*qa + rw*(chg*uc)*(betg*dqs0g+dbetg*qs0g) ) * qc1 &
1551  - ( rw*(dchcdtg*ua) + rw*(chg*uc*dbetg) ) * qc2 ) / ( qc1**2 + eps )
1552  dqcdac = ( ( rw*(dchcdac*ua)*qa ) * qc1 &
1553  - ( rw*(dchcdac*ua) ) * qc2 ) / ( qc1**2 + eps )
1554 
1555  rg1 = epsg * ( rflux_lw * vfgs &
1556  + epsb * vfgw * stb * tb**4 &
1557  - stb * tg**4 )
1558  drg1dtb = epsg * epsb * vfgw * stb * 4.0_rp * tb**3
1559  drg1dtg = - epsg * stb * 4.0_rp * tg**3
1560 
1561  rb1 = epsb * ( rflux_lw * vfws &
1562  + epsg * vfwg * stb * tg**4 &
1563  + epsb * vfww * stb * tb**4 &
1564  - stb * tb**4 )
1565  drb1dtb = epsb * ( epsb * vfww - 1.0_rp ) * stb * 4.0_rp * tb**3
1566  drb1dtg = epsb * epsg * vfwg * stb * 4.0_rp * tg**3
1567 
1568  rg2 = epsg * ( (1.0_rp-epsb) * vfgw * vfws * rflux_lw &
1569  + (1.0_rp-epsb) * vfgw * vfwg * epsg * stb * tg**4 &
1570  + epsb * (1.0_rp-epsb) * vfgw * vfww * stb * tb**4 )
1571  drg2dtb = epsg * epsb * (1.0_rp-epsb) * vfgw * vfww * stb * 4.0_rp * tb**3
1572  drg2dtg = epsg * (1.0_rp-epsb) * vfgw * vfwg * epsg * stb * 4.0_rp * tg**3
1573 
1574  rb2 = epsb * ( (1.0_rp-epsg) * vfwg * vfgs * rflux_lw &
1575  + (1.0_rp-epsb) * vfws * vfww * rflux_lw &
1576  + (1.0_rp-epsb) * epsg * vfwg * vfww * stb * tg**4 &
1577  + (1.0_rp-epsg) * epsb * vfgw * vfwg * stb * tb**4 &
1578  + (1.0_rp-epsb) * epsb * vfww * (1.0_rp-2.0_rp*vfws) * stb * tb**4 )
1579  drb2dtb = epsb**2 * ( (1.0_rp-epsg) * vfgw * vfwg + (1.0_rp-epsb) * vfww * (1.0_rp-2.0_rp*vfws) ) * stb * 4.0_rp * tb**3
1580  drb2dtg = epsb * (1.0_rp-epsb) * epsg * vfwg * vfww * stb * 4.0_rp * tg**3
1581 
1582  rg = rg1 + rg2
1583  rb = rb1 + rb2
1584 
1585  hb = rhoo * cpdry * chb * uc * ( ths1 - thc ) * exn
1586  dhbdtb = rhoo * cpdry * chb * uc * ( 1.0_rp - dthcdtb * exn )
1587  dhbdtg = rhoo * cpdry * chb * uc * ( - dthcdtg * exn )
1588  dhbdac = rhoo * cpdry * chb * uc * ( - dthcdac * exn )
1589 
1590  evpb = min( rhoo * chb * uc * betb * (qs0b-qc), rainb / dt_rp )
1591  eleb = evpb * lhv
1592  delebdtb = rhoo * chb * uc * ( betb * ( dqs0b - dqcdtb ) + dbetb * (qs0b-qc) ) * lhv
1593  delebdtg = rhoo * chb * uc * betb * ( - dqcdtg ) * lhv
1594  delebdac = rhoo * chb * uc * betb * ( - dqcdac ) * lhv
1595 
1596 ! G0B = SB + RB - HB - ELEB + EFLX
1597  g0b = sb + rb - hb - eleb
1598  dg0bdtb = drb1dtb + drb2dtb - dhbdtb - delebdtb
1599  dg0bdtg = drb1dtg + drb2dtg - dhbdtg - delebdtg
1600  dg0bdac = - dhbdac - delebdac
1601 
1602  hg = rhoo * cpdry * chg * uc * ( ths2 - thc ) * exn
1603  dhgdtb = rhoo * cpdry * chg * uc * ( - dthcdtb * exn )
1604  dhgdtg = rhoo * cpdry * chg * uc * ( 1.0_rp - dthcdtg * exn )
1605  dhgdac = rhoo * cpdry * chg * uc * ( - dthcdac * exn )
1606 
1607  evpg = min( rhoo * chg * uc * betg * (qs0g-qc), raing / dt_rp )
1608  eleg = evpg * lhv
1609  delegdtb = rhoo * chg * uc * betg * ( - dqcdtb ) * lhv
1610  delegdtg = rhoo * chg * uc * ( betg * ( dqs0g - dqcdtg ) + dbetg * (qs0g-qc) ) * lhv
1611  delegdac = rhoo * chg * uc * betg * ( - dqcdac ) * lhv
1612 
1613 ! G0G = SG + RG - HG - ELEG + EFLX
1614  g0g = sg + rg - hg - eleg
1615  dg0gdtb = drg1dtb + drg2dtb - dhgdtb - delegdtb
1616  dg0gdtg = drg1dtg + drg2dtg - dhgdtg - delegdtg
1617  dg0gdac = - dhgdac - delegdac
1618 
1619  do k = uks, uke
1620  tbl(k) = tblp(k)
1621  end do
1622  call multi_layer(uke,bound, &
1623  a, b, c, d, p, q, &
1624  g0b,capb,aksb,tbl,dzb,dt_rp,tblend)
1625  resi1 = tbl(1) - tb
1626 
1627  do k = uks, uke
1628  tgl(k) = tglp(k)
1629  end do
1630  call multi_layer(uke,bound, &
1631  a, b, c, d, p, q, &
1632  g0g,capg,aksg,tgl,dzg,dt_rp,tglend)
1633  resi2 = tgl(1) - tg
1634 
1635 
1636  resi3 = alphac - alphacp
1637 
1638  if ( abs(resi1) < threshold .AND. abs(resi2) < threshold .AND. abs(resi3) < threshold ) then
1639  tb = tbl(1)
1640  tg = tgl(1)
1641  tb = max( tbp - dts_max_onestep, min( tbp + dts_max_onestep, tb ) )
1642  tg = max( tgp - dts_max_onestep, min( tgp + dts_max_onestep, tg ) )
1643  exit
1644  endif
1645 
1646  ! Newton method
1647  dr1dtb = dtbdg0b * dg0bdtb - 1.0_rp
1648  dr1dtg = dtbdg0b * dg0bdtg
1649  dr1dac = dtbdg0b * dg0bdac
1650  dr2dtb = dtgdg0g * dg0gdtb
1651  dr2dtg = dtgdg0g * dg0gdtg - 1.0_rp
1652  dr2dac = dtgdg0g * dg0gdac
1653  dr3dtb = dalphacdtb
1654  dr3dtg = dalphacdtg
1655  dr3dac = dalphacdac - 1.0_rp
1656  if ( iteration > 3 ) dr3dac = min( dr3dac, -0.02_rp )
1657 
1658  rdet = 1.0_rp &
1659  / ( dr1dtb * dr2dtg * dr3dac + dr1dtg * dr2dac * dr3dtb + dr1dac * dr2dtb * dr3dtg &
1660  - dr1dac * dr2dtg * dr3dtb - dr1dtg * dr2dtb * dr3dac - dr1dtb * dr2dac * dr3dtg )
1661  dtb = ( - ( dr2dtg * dr3dac - dr2dac * dr3dtg ) * resi1 &
1662  + ( dr1dtg * dr3dac - dr1dac * dr3dtg ) * resi2 &
1663  - ( dr1dtg * dr2dac - dr1dac * dr2dtg ) * resi3 ) * rdet
1664  dtg = ( ( dr2dtb * dr3dac - dr2dac * dr3dtb ) * resi1 &
1665  - ( dr1dtb * dr3dac - dr1dac * dr3dtb ) * resi2 &
1666  + ( dr1dtb * dr2dac - dr1dac * dr2dtb ) * resi3 ) * rdet
1667  dac = ( - ( dr2dtb * dr3dtg - dr2dtg * dr3dtb ) * resi1 &
1668  + ( dr1dtb * dr3dtg - dr1dtg * dr3dtb ) * resi2 &
1669  - ( dr1dtb * dr2dtg - dr1dtg * dr2dtb ) * resi3 ) * rdet
1670 
1671  dtb = sign( min(abs(dtb), 5.0_rp), dtb )
1672  dtg = sign( min(abs(dtg), 5.0_rp), dtg )
1673  dac = sign( min(abs(dac), 50.0_rp), dac )
1674 
1675  if ( iteration > 50 ) then
1676  if ( dtb*dtbp < 0.0_rp ) then
1677  fact1 = max( fact1 * 0.5_rp, 1e-10_rp )
1678  else
1679  fact1 = min( fact1 * 1.5_rp, 1.1_rp )
1680  end if
1681  if ( dtg*dtgp < 0.0_rp ) then
1682  fact2 = max( fact2 * 0.5_rp, 1e-10_rp )
1683  else
1684  fact2 = min( fact2 * 1.5_rp, 1.1_rp )
1685  end if
1686  if ( dac*dacp < 0.0_rp ) then
1687  fact3 = max( fact3 * 0.5_rp, 1e-10_rp )
1688  else
1689  fact3 = min( fact3 * 1.5_rp, 1.1_rp )
1690  end if
1691  end if
1692 
1693  tmp = tb
1694  tb = max( tb + dtb * fact1, 100.0_rp )
1695  dtbp = tb - tmp
1696  tmp = tg
1697  tg = max( tg + dtg * fact2, 100.0_rp )
1698  dtgp = tg - tmp
1699  alphac = max( alphacp + dac * fact3, eps )
1700  dacp = alphac - alphacp
1701 
1702  if ( abs(dtbp) < threshold .AND. abs(dtgp) < threshold .AND. abs(dacp) < threshold ) then
1703  tb = tbl(1)
1704  tg = tgl(1)
1705  tb = max( tbp - dts_max_onestep, min( tbp + dts_max_onestep, tb ) )
1706  tg = max( tgp - dts_max_onestep, min( tgp + dts_max_onestep, tg ) )
1707  exit
1708  endif
1709 
1710  alphacp = alphac
1711 
1712  evpb = ( evpbp + evpb ) * 0.5_rp
1713  evpbp = evpb
1714 
1715  evpg = ( evpgp + evpg ) * 0.5_rp
1716  evpgp = evpg
1717 
1718  enddo
1719 
1720 #ifdef DEBUG_KUSAKA01
1721  cnt_num2 = cnt_num2 + 1
1722  cnt_itr2 = cnt_itr2 + iteration - 1
1723  max_itr2 = max(max_itr2, iteration - 1)
1724 #endif
1725 
1726  ! output for debug
1727  if ( iteration > itr_max .and. debug ) then
1728  log_warn("URBAN_DYN_Kusaka01_main",*) 'iteration for TB/TG was not converged',prc_myrank,i,j
1729  log_warn_cont(*) '---------------------------------------------------------------------------------'
1730  log_warn_cont(*) 'DEBUG Message --- Residual [K] :', resi1, resi2, resi3
1731  log_warn_cont(*) 'DEBUG Message --- TBP : Initial TB [K] :', tbp
1732 #ifdef _OPENACC
1733  log_warn_cont(*) 'DEBUG Message --- TBLP: Initial TBL [K] :', tblp(uks)
1734 #else
1735  log_warn_cont(*) 'DEBUG Message --- TBLP: Initial TBL [K] :', tblp(:)
1736 #endif
1737  log_warn_cont(*) 'DEBUG Message --- TGP : Initial TG [K] :', tgp
1738 #ifdef _OPENACC
1739  log_warn_cont(*) 'DEBUG Message --- TGLP: Initial TGL [K] :', tglp(uks)
1740 #else
1741  log_warn_cont(*) 'DEBUG Message --- TGLP: Initial TGL [K] :', tglp(:)
1742 #endif
1743  log_warn_cont(*) 'DEBUG Message --- TCP : Initial TC [K] :', tcp
1744  log_warn_cont(*) 'DEBUG Message --- QCP : Initial QC [K] :', qcp
1745  log_warn_cont(*) 'DEBUG Message --- UC : Canopy wind [m/s] :', uc
1746  log_warn_cont(*) 'DEBUG Message --- rflux_SW : Shortwave radiation [W/m2] :', rflux_sw
1747  log_warn_cont(*) 'DEBUG Message --- rflux_LW : Longwave radiation [W/m2] :', rflux_lw
1748  log_warn_cont(*) 'DEBUG Message --- PRSS: Surface pressure [Pa] :', prss
1749  log_warn_cont(*) 'DEBUG Message --- PRSA: Pressure at 1st atmos layer [m] :', prsa
1750  log_warn_cont(*) 'DEBUG Message --- RHOO: Air density [kg/m3] :', rhoo
1751  log_warn_cont(*) 'DEBUG Message --- RHOS: Surface density [kg/m3] :', rhos
1752  log_warn_cont(*) 'DEBUG Message --- RAINBP: Initial RAINB [kg/m2] :', rainbp
1753  log_warn_cont(*) 'DEBUG Message --- RAINGP: Initial RAING [kg/m2] :', raingp
1754  log_warn_cont(*) 'DEBUG Message --- ZA : Height at 1st atmos layer [m] :', za
1755  log_warn_cont(*) 'DEBUG Message --- TA : Temperature at 1st atmos layer [K] :', ta
1756  log_warn_cont(*) 'DEBUG Message --- UA : Wind speed at 1st atmos layer [m/s] :', ua
1757  log_warn_cont(*) 'DEBUG Message --- QA : Specific humidity at 1st atmos layer [kg/kg] :', qa
1758 #ifdef _OPENACC
1759  log_warn_cont(*) 'DEBUG Message --- DZB : Depth of surface layer [m] :', dzb(1)
1760  log_warn_cont(*) 'DEBUG Message --- DZG : Depth of surface layer [m] :', dzg(1)
1761 #else
1762  log_warn_cont(*) 'DEBUG Message --- DZB : Depth of surface layer [m] :', dzb(:)
1763  log_warn_cont(*) 'DEBUG Message --- DZG : Depth of surface layer [m] :', dzg(:)
1764 #endif
1765  log_warn_cont(*) 'DEBUG Message --- R, W, RW : Normalized height and road width [-] :', r, w,rw
1766  log_warn_cont(*) 'DEBUG Message --- SVF : Sky View Factors [-] :', svf
1767  log_warn_cont(*) 'DEBUG Message --- BETB,BETG : Evaporation efficiency [-] :', betb,betg
1768  log_warn_cont(*) 'DEBUG Message --- EPSB,EPSG : Surface emissivity [-] :', epsb,epsg
1769  log_warn_cont(*) 'DEBUG Message --- CAPB,CAPG : Heat capacity [J m-3 K] :', capb,capg
1770  log_warn_cont(*) 'DEBUG Message --- AKSB,AKSG : Thermal conductivity [W m-1 K] :', aksb,aksb
1771  log_warn_cont(*) 'DEBUG Message --- QS0B,QS0G : Surface specific humidity [kg/kg] :', qs0b,qs0g
1772  log_warn_cont(*) 'DEBUG Message --- ZDC : Desplacement height of canopy [m] :', zdc
1773  log_warn_cont(*) 'DEBUG Message --- Z0M : Momentum roughness length of canopy [m] :', z0c
1774  log_warn_cont(*) 'DEBUG Message --- Z0H/Z0E : Thermal roughness length of canopy [m] :', z0hc
1775  log_warn_cont(*) 'DEBUG Message --- LHV : Latent heat of vapor [J kg-1] :', lhv
1776  log_warn_cont(*) 'DEBUG Message --- dt : Time step [s] :', dt_rp
1777  log_warn_cont(*) '---------------------------------------------------------------------------------'
1778  endif
1779 
1780  !--- update only fluxes ----
1781 
1782  ! this is for TC, QC
1783  call qsat( tb, rhos, qs0b )
1784  call qsat( tg, rhos, qs0g )
1785 
1786  call cal_beta(betb, betb_const, rainb, strgb)
1787  call cal_beta(betg, betg_const, raing, strgg)
1788 
1789 
1790  rg1 = epsg * ( rflux_lw * vfgs &
1791  + epsb * vfgw * stb * tb**4 &
1792  - stb * tg**4 )
1793  rb1 = epsb * ( rflux_lw * vfws &
1794  + epsg * vfwg * stb * tg**4 &
1795  + epsb * vfww * stb * tb**4 &
1796  - stb * tb**4 )
1797 
1798  rg2 = epsg * ( (1.0_rp-epsb) * vfgw * vfws * rflux_lw &
1799  + (1.0_rp-epsb) * vfgw * vfwg * epsg * stb * tg**4 &
1800  + epsb * (1.0_rp-epsb) * vfgw * vfww * stb * tb**4 )
1801  rb2 = epsb * ( (1.0_rp-epsg) * vfwg * vfgs * rflux_lw &
1802  + (1.0_rp-epsg) * epsb * vfgw * vfwg * stb * tb**4 &
1803  + (1.0_rp-epsb) * vfws * vfww * rflux_lw &
1804  + (1.0_rp-epsb) * vfwg * vfww * stb * epsg * tg**4 &
1805  + epsb * (1.0_rp-epsb) * vfww * (1.0_rp-2.0_rp*vfws) * stb * tb**4 )
1806 
1807  rg = rg1 + rg2
1808  rb = rb1 + rb2
1809 
1810  ths1 = tb / exn
1811  ths2 = tg / exn
1812  thc = tc / exn
1813 
1814  hb = rhoo * cpdry * chb * uc * (ths1-thc) * exn
1815  evpb = min( rhoo * chb * uc * betb * (qs0b-qc), rainb / dt_rp )
1816  eleb = evpb * lhv
1817 ! G0B = SB + RB - HB - ELEB + EFLX
1818  g0b = sb + rb - hb - eleb
1819  rainb = max( rainbp - evpb * dt_rp, 0.0_rp )
1820 
1821  hg = rhoo * cpdry * chg * uc * (ths2-thc) * exn
1822  evpg = min( rhoo * chg * uc * betg * (qs0g-qc), raing / dt_rp )
1823  eleg = evpg * lhv
1824 ! G0G = SG + RG - HG - ELEG + EFLX
1825  g0g = sg + rg - hg - eleg
1826  raing = max( raingp - evpg * dt_rp, 0.0_rp )
1827 
1828  do k = uks, uke
1829  tbl(k) = tblp(k)
1830  end do
1831  call multi_layer(uke,bound, &
1832  a, b, c, d, p, q, &
1833  g0b,capb,aksb,tbl,dzb,dt_rp,tblend)
1834  resi1 = tbl(1) - tb
1835  tb = tbl(1)
1836 
1837  do k = uks, uke
1838  tgl(k) = tglp(k)
1839  end do
1840  call multi_layer(uke,bound, &
1841  a, b, c, d, p, q, &
1842  g0g,capg,aksg,tgl,dzg,dt_rp,tglend)
1843  resi2 = tgl(1) - tg
1844  tg = tgl(1)
1845 
1846  if ( abs(resi1) > dts_max_onestep ) then
1847  if ( abs(resi1) > dts_max_onestep*10.0_rp ) then
1848  log_error("URBAN_DYN_Kusaka01_main",*) 'tendency of TB exceeded a limit! STOP.'
1849  log_error_cont(*) 'previous TB and updated TB(TBL(1)) is ', tb-resi1,tb
1850  converged = .false.
1851 #ifdef _OPENACC
1852  return
1853 #else
1854  call prc_abort
1855 #endif
1856  endif
1857  log_warn("URBAN_DYN_Kusaka01_main",*) 'tendency of TB exceeded a limit'
1858  log_warn_cont(*) 'previous TB and updated TB(TBL(1)) is ', tb-resi1, tb
1859  endif
1860 
1861  if ( abs(resi2) > dts_max_onestep ) then
1862  if ( abs(resi2) > dts_max_onestep*10.0_rp ) then
1863  log_error("URBAN_DYN_Kusaka01_main",*) 'tendency of TG exceeded a limit! STOP.'
1864  log_error_cont(*) 'previous TG and updated TG(TGL(1)) is ', tg-resi2, tg, resi2
1865  converged = .false.
1866 #ifdef _OPENACC
1867  return
1868 #else
1869  call prc_abort
1870 #endif
1871  endif
1872  log_warn("URBAN_DYN_Kusaka01_main",*) 'tendency of TG exceeded a limit'
1873  log_warn_cont(*) 'previous TG and updated TG(TGL(1)) is ', tg-resi2, tg
1874  endif
1875 
1876  !-----------------------------------------------------------
1877  ! Total Fluxes from Urban Canopy
1878  !-----------------------------------------------------------
1879 
1880  flxuv = ( r*cdr + rw*cdc ) * ua * ua
1881  mflx = - rhoo * flxuv ! Momentum flux [kg/m/s2]
1882  sh = ( r*hr + w*hb + rw*hg ) ! Sensible heat flux [W/m/m]
1883  lh = ( r*eler + w*eleb + rw*eleg ) ! Latent heat flux [W/m/m]
1884  ghflx = r*g0r + w*g0b + rw*g0g
1885  lnet = r*rr + w*rb + rw*rg
1886 
1887  !-----------------------------------------------------------
1888  ! Grid average
1889  !-----------------------------------------------------------
1890 
1891  lw = rflux_lw - lnet ! Upward longwave radiation [W/m/m]
1892  sw = rflux_sw - snet ! Upward shortwave radiation [W/m/m]
1893  rn = (snet+lnet) ! Net radiation [W/m/m]
1894 
1895  !--- shortwave radiation
1896  sdn = r + w * (vfws + vfgs * albg * vfwg) + rw * (vfgs + vfws * albb * vfgw)
1897  sup = r * albr &
1898  + w * ( vfws * albb + vfgs * albg * vfwg *albb ) &
1899  + rw * ( vfgs * albg + vfws * albb * vfgw *albg )
1900 
1901  albd_sw_grid = sup / sdn
1902 
1903  !--- longwave radiation
1904  ldn = r + w*vfws + rw*vfgs
1905  lup = r * (1.0_rp-epsr) &
1906  + w*( (1.0_rp-epsb*vfww)*(1.0_rp-epsb)*vfws - epsb*vfwg*(1.0_rp-epsg)*vfgs ) &
1907  + rw*( (1.0_rp-epsg)*vfgs - epsg*(1.0_rp-vfgs)*(1.0_rp-epsb)*vfws )
1908 
1909  rup = (ldn - lup) * rflux_lw - lnet
1910  albd_lw_grid = lup / ldn
1911 
1912 
1913  ! RUP = R * (EPSR * STB * TR**4 ) &
1914  ! + W * (EPSB*STB*(TB**4) - EPSB*EPSG*VFWG*STB*(TG**4) - EPSB*EPSB*VFWW*STB*(TB**4) &
1915  ! - EPSB *(1.0_RP-EPSG) * EPSB * VFGW * VFWG * STB * (TB**4) &
1916  ! - EPSB *(1.0_RP-EPSB) * VFWG * VFWW * STB * EPSG * (TG**4) &
1917  ! - EPSB * EPSB * (1.0_RP-EPSB) * VFWW * VFWW * STB * (TB**4) ) &
1918  ! + RW * (EPSG*STB*(TG**4) - EPSG * EPSB * VFGW * STB * (TB**4) &
1919  ! - EPSG * EPSB * (1.0_RP-EPSB) * (1.0_RP-SVF) * VFWW * STB * TB**4 )
1920 
1921 
1922  shr = hr ! Sensible heat flux on roof [W/m/m]
1923  shb = hb ! Sensible heat flux on wall [W/m/m]
1924  shg = hg ! Sensible heat flux on road [W/m/m]
1925  lhr = eler ! Latent heat flux on road [W/m/m]
1926  lhb = eleb ! Latent heat flux on wall [W/m/m]
1927  lhg = eleg ! Latent heat flux on road [W/m/m]
1928  ghr = g0r ! Ground heat flux on roof [W/m/m]
1929  ghb = g0b ! Ground heat flux on wall [W/m/m]
1930  ghg = g0g ! Ground heat flux on road [W/m/m]
1931  rnr = sr + rr ! Net radiation on roof [W/m/m]
1932  rnb = sb + rb ! Net radiation on building [W/m/m]
1933  rng = sg + rg ! Net radiation on ground [W/m/m]
1934 
1935  !-----------------------------------------------------------
1936  ! calculate the ranoff
1937  !-----------------------------------------------------------
1938  roffr = max( rainr - strgr, 0.0_rp )
1939  roffb = max( rainb - strgb, 0.0_rp )
1940  roffg = max( raing - strgg, 0.0_rp )
1941 
1942  roff = r * roffr + rw * ( roffb + roffg )
1943 
1944  rainr = max( rainr - roffr, 0.0_rp )
1945  rainb = max( rainb - roffb, 0.0_rp )
1946  raing = max( raing - roffg, 0.0_rp )
1947 
1948  !-----------------------------------------------------------
1949  ! diagnostic GRID AVERAGED TS from upward logwave
1950  !-----------------------------------------------------------
1951 
1952  rts = ( rup / stb / ( 1.0_rp-albd_lw_grid) )**0.25
1953 
1954  !-----------------------------------------------------------
1955  ! diagnostic grid average U10, V10, T2, Q2 from urban
1956  ! Below method would be better to be improved. This is tentative method.
1957  !-----------------------------------------------------------
1958  ustar = sqrt( flxuv ) ! u* [m/s]
1959  tstar = -sh / rhoo / cpdry / ustar ! T* [K]
1960  qstar = -lh / rhoo / lhv / ustar ! q* [-]
1961  !Z = ZA - ZDC
1962  !XXX = 0.4*9.81*Z*TST/TA/UST/UST
1963 
1964  xxx = xxxc
1965  call cal_psi(xxx,psim,psih)
1966 
1967  xxx2 = (2.0_rp/z) * xxx
1968  call cal_psi(xxx2,psim2,psih2)
1969 
1970  xxx10 = (10.0_rp/z) * xxx
1971  call cal_psi(xxx10,psim10,psih10)
1972 
1973 
1974  fracu10 = ( log(10.0_rp/z0c ) - psim10 ) / ( log(z/z0c ) - psim )
1975  fract2 = ( log( 2.0_rp/z0hc) - psih2 ) / ( log(z/z0hc) - psih )
1976 
1977  call bulkflux_diagnose_surface( u1, v1, & ! [IN]
1978  ta, qa, & ! [IN]
1979  rts, 1.0_rp, & ! [IN]
1980  z, & ! [IN]
1981  z0c, z0hc, 1.0_rp, & ! [IN]
1982  u10, v10, t2, q2, & ! [OUT] ! Q2 is dummy
1983  fracu10, fract2, 1.0_rp ) ! [IN]
1984 
1985  q2 = qc
1986 
1987  return
1988  end subroutine slc_main
1989 
1990  !-----------------------------------------------------------------------------
1991  subroutine canopy_wind(ZA, UA, Z0C, ZDC, UC)
1992  !$acc routine
1993  implicit none
1994 
1995  real(rp), intent(in) :: za ! height at 1st atmospheric level [m]
1996  real(rp), intent(in) :: ua ! wind speed at 1st atmospheric level [m/s]
1997  real(rp), intent(in) :: z0c ! Roughness length above canyon for momentum [m]
1998  real(rp), intent(in) :: zdc ! Displacement height [m]
1999  real(rp), intent(out) :: uc ! wind speed at 1st atmospheric level [m/s]
2000 
2001  real(rp) :: ur,zc,xlb,bb
2002 
2003  if( zr + 2.0_rp < za ) then
2004  ur = ua * log((zr-zdc)/z0c) / log((za-zdc)/z0c)
2005  zc = 0.7_rp * zr
2006  xlb = 0.4_rp * (zr-zdc)
2007  ! BB formulation from Inoue (1963)
2008  bb = 0.4_rp * zr / ( xlb * log((zr-zdc)/z0c) )
2009  uc = ur * exp( -bb * (1.0_rp-zc/zr) )
2010  else
2011  ! PRINT *, 'Warning ZR + 2m is larger than the 1st WRF level'
2012  zc = za / 2.0_rp
2013  uc = ua / 2.0_rp
2014  endif
2015 
2016  uc = max(uc,0.01_rp)
2017 
2018  return
2019  end subroutine canopy_wind
2020 
2021  !-----------------------------------------------------------------------------
2022  subroutine cal_beta(BET, BET_CONST, WATER, STRG)
2023  !$acc routine
2024  implicit none
2025 
2026  real(rp), intent(out) :: bet ! evapolation efficiency [-]
2027  real(rp), intent(in) :: bet_const ! prescribed value
2028  real(rp), intent(in) :: water ! rain amount in strage [kg/m2]
2029  real(rp), intent(in) :: strg ! rain strage [kg/m2]
2030 
2031  if ( bet_const > 0.0_rp ) then
2032  bet = bet_const
2033  else if ( strg == 0.0_rp ) then ! not consider evapolation from urban
2034  bet = 0.0_rp
2035  else
2036  bet = max( min( water / strg, 1.0_rp ), &
2037  1.0e-10_rp ) ! When WATER < STRG/1e10, fix the beta value so that tiny amoumts of water do not coninue to remain
2038  endif
2039 
2040  return
2041  end subroutine cal_beta
2042 
2043  !-----------------------------------------------------------------------------
2044  subroutine cal_psi(zeta,psim,psih)
2045  !$acc routine
2046  use scale_const, only: &
2047  pi => const_pi
2048  implicit none
2049 
2050  real(rp), intent(inout) :: zeta ! z/L
2051  real(rp), intent(out) :: psim
2052  real(rp), intent(out) :: psih
2053  real(rp) :: x
2054 
2055  if( zeta >= 1.0_rp ) zeta = 1.0_rp
2056  if( zeta <= -5.0_rp ) zeta = -5.0_rp
2057 
2058  if( zeta > 0.0_rp ) then
2059  psim = -5.0_rp * zeta
2060  psih = -5.0_rp * zeta
2061  else
2062  x = ( 1.0_rp - 16.0_rp * zeta )**0.25_rp
2063  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
2064  psih = 2.0_rp * log((1.0_rp+x*x)/2.0_rp)
2065  end if
2066 
2067  return
2068  end subroutine cal_psi
2069 
2070  !-----------------------------------------------------------------------------
2071  ! XXX: z/L (requires iteration by Newton-Rapson method)
2072  ! B1: Stanton number
2073  ! PSIM: = PSIX of LSM
2074  ! PSIH: = PSIT of LSM
2075 !OCL SERIAL
2076  subroutine mos(XXX,CH,CD,B1,RIB,Z,Z0,UA,TA,TSF,RHO,i,j)
2077  !$acc routine
2078  use scale_const, only: &
2079  eps => const_eps, &
2080  cpdry => const_cpdry ! CPP : heat capacity of dry air [J/K/kg]
2081  implicit none
2082 integer,intent(in)::i,j
2083  real(rp), intent(in) :: b1, z, z0, ua, ta, tsf, rho
2084  real(rp), intent(out) :: cd, ch
2085  real(rp), intent(inout) :: xxx, rib
2086  real(rp) :: xxx0, x, x0, faih, dpsim, dpsih
2087  real(rp) :: f, df, xxxp, us, ts, al, xkb, dd, psim, psih
2088  integer :: newt
2089  integer, parameter :: newt_end = 10
2090 
2091  real(rp) :: lnz
2092  real(rp) :: sqx, sqx0
2093 
2094  lnz = log( (z+z0)/z0 )
2095 
2096  if( rib < 0.0_rp ) then
2097 
2098  rib = max( rib, -15.0_rp )
2099 
2100  do newt = 1, newt_end
2101 
2102  xxx = min( xxx, -1.0e-3_rp )
2103 
2104  xxx0 = xxx * z0/(z+z0)
2105 
2106  sqx = sqrt( 1.0_rp - 16.0_rp * xxx )
2107  sqx0 = sqrt( 1.0_rp - 16.0_rp * xxx0 )
2108 
2109  x = sqrt( sqx )
2110  x0 = sqrt( sqx0 )
2111 
2112  psim = lnz &
2113  - log( (x+1.0_rp)**2 * (sqx + 1.0_rp) ) &
2114  + 2.0_rp * atan(x) &
2115  + log( (x+1.0_rp)**2 * (sqx0 + 1.0_rp) ) &
2116  - 2.0_rp * atan(x0)
2117  faih = 1.0_rp / sqx
2118  psih = lnz + 0.4_rp*b1 &
2119  - 2.0_rp * log( sqx + 1.0_rp ) &
2120  + 2.0_rp * log( sqx0 + 1.0_rp )
2121 
2122  dpsim = 1.0_rp / ( x * xxx ) &
2123  - 1.0_rp / ( x0 * xxx )
2124  dpsih = 1.0_rp / ( sqx * xxx ) &
2125  - 1.0_rp / ( sqx0 * xxx )
2126 
2127  f = rib * psim**2 / psih - xxx
2128 
2129  df = rib * ( 2.0_rp*dpsim*psim*psih - dpsih*psim**2 ) &
2130  / psih**2 - 1.0_rp
2131 
2132  xxxp = xxx
2133  xxx = xxxp - f / df
2134 
2135  xxx = max( xxx, -10.0_rp )
2136 
2137  end do
2138 
2139  else if( rib >= 0.142857_rp ) then
2140 
2141  xxx = 0.714_rp
2142  psim = lnz + 7.0_rp * xxx
2143  psih = psim + 0.4_rp * b1
2144 
2145  else
2146 
2147  al = lnz
2148  xkb = 0.4_rp * b1
2149  dd = -4.0_rp * rib * 7.0_rp * xkb * al + (al+xkb)**2
2150  dd = max( dd, 0.0_rp )
2151 
2152  xxx = ( al + xkb - 2.0_rp*rib*7.0_rp*al - sqrt(dd) ) / ( 2.0_rp * ( rib*7.0_rp**2 - 7.0_rp ) )
2153  xxx = min( xxx, 0.714_rp )
2154  psim = lnz + 7.0_rp * xxx
2155  psih = psim + 0.4_rp * b1
2156 
2157  endif
2158 
2159  us = 0.4_rp * ua / psim ! u*
2160  if( us <= 0.01_rp ) us = 0.01_rp
2161  !TS = 0.4_RP * (TA-TSF) / PSIH ! T*
2162 
2163  cd = us * us / (ua+eps)**2 ! CD
2164  ch = 0.4_rp * us / psih / (ua+eps) ! CH
2165  !ALPHA = RHO * CPdry * 0.4_RP * US / PSIH ! RHO*CP*CH*U
2166 
2167  return
2168  end subroutine mos
2169 
2170  !-------------------------------------------------------------------
2171 !OCL SERIAL
2172  subroutine multi_layer( &
2173  KM,BOUND, &
2174  A, B, C, D, P, Q, &
2175  G0,CAP,AKS,TSL,DZ,DELT,TSLEND)
2176  !
2177  ! calculate temperature in roof/building/road
2178  ! multi-layer heat equation model
2179  ! Solving Heat Equation by Tri Diagonal Matrix Algorithm
2180  !-------------------------------------------------------------------
2181  !$acc routine seq
2182  implicit none
2183  integer, intent(in) :: km
2184  integer, intent(in) :: bound
2185  real(rp), intent(in) :: g0
2186  real(rp), intent(in) :: cap
2187  real(rp), intent(in) :: aks
2188  real(rp), intent(inout) :: tsl(km)
2189  real(rp), intent(in) :: dz(km)
2190  real(rp), intent(in) :: delt ! Tim setep [ s ]
2191  real(rp), intent(in) :: tslend
2192  real(rp), intent(out) :: a(km), b(km), c(km), d(km), p(km), q(km)
2193 
2194  real(rp) :: dzend
2195  integer :: k
2196 
2197  dzend = dz(km)
2198 
2199  a(1) = 0.0_rp
2200 
2201  b(1) = cap * dz(1) / delt &
2202  + 2.0_rp * aks / (dz(1)+dz(2))
2203  c(1) = -2.0_rp * aks / (dz(1)+dz(2))
2204  d(1) = cap * dz(1) / delt * tsl(1) + g0
2205 
2206  do k = 2, km-1
2207  a(k) = -2.0_rp * aks / (dz(k-1)+dz(k))
2208  b(k) = cap * dz(k) / delt + 2.0_rp * aks / (dz(k-1)+dz(k)) + 2.0_rp * aks / (dz(k)+dz(k+1))
2209  c(k) = -2.0_rp * aks / (dz(k)+dz(k+1))
2210  d(k) = cap * dz(k) / delt * tsl(k)
2211  end do
2212 
2213  if( bound == 1 ) then ! Flux=0
2214  a(km) = -2.0_rp * aks / (dz(km-1)+dz(km))
2215  b(km) = cap * dz(km) / delt + 2.0_rp * aks / (dz(km-1)+dz(km))
2216  c(km) = 0.0_rp
2217  d(km) = cap * dz(km) / delt * tsl(km)
2218  else ! T=constant
2219  a(km) = -2.0_rp * aks / (dz(km-1)+dz(km))
2220  b(km) = cap * dz(km) / delt + 2.0_rp * aks / (dz(km-1)+dz(km)) + 2.0_rp * aks / (dz(km)+dzend)
2221  c(km) = 0.0_rp
2222  d(km) = cap * dz(km) / delt * tsl(km) + 2.0_rp * aks * tslend / (dz(km)+dzend)
2223  end if
2224 
2225  p(1) = -c(1) / b(1)
2226  q(1) = d(1) / b(1)
2227 
2228  !$acc loop seq
2229  do k = 2, km
2230  p(k) = -c(k) / ( a(k) * p(k-1) + b(k) )
2231  q(k) = ( -a(k) * q(k-1) + d(k) ) / ( a(k) * p(k-1) + b(k) )
2232  end do
2233 
2234  tsl(km) = q(km)
2235 
2236  !$acc loop seq
2237  do k = km-1, 1, -1
2238  tsl(k) = p(k) * tsl(k+1) + q(k)
2239  end do
2240 
2241  return
2242  end subroutine multi_layer
2243 
2244 !!$ !-------------------------------------------------------------------
2245 !!$!OCL SERIAL
2246 !!$ subroutine multi_layer2(KM,BOUND,G0,CAP,AKS,TSL,DZ,DELT,TSLEND,CAP1,AKS1)
2247 !!$ !
2248 !!$ ! calculate temperature in roof/building/road
2249 !!$ ! multi-layer heat equation model
2250 !!$ ! Solving Heat Equation by Tri Diagonal Matrix Algorithm
2251 !!$ !-------------------------------------------------------------------
2252 !!$
2253 !!$ implicit none
2254 !!$
2255 !!$ real(RP), intent(in) :: G0
2256 !!$ real(RP), intent(in) :: CAP
2257 !!$ real(RP), intent(in) :: AKS
2258 !!$ real(RP), intent(in) :: CAP1 ! for 1st layer
2259 !!$ real(RP), intent(in) :: AKS1 ! for 1st layer
2260 !!$ real(DP), intent(in) :: DELT ! Time step [ s ]
2261 !!$ real(RP), intent(in) :: TSLEND
2262 !!$ integer, intent(in) :: KM
2263 !!$ integer, intent(in) :: BOUND
2264 !!$ real(RP), intent(in) :: DZ(KM)
2265 !!$ real(RP), intent(inout) :: TSL(KM)
2266 !!$ real(RP) :: A(KM), B(KM), C(KM), D(KM), X(KM), P(KM), Q(KM)
2267 !!$ real(RP) :: DZEND
2268 !!$ integer :: K
2269 !!$
2270 !!$ DZEND = DZ(KM)
2271 !!$
2272 !!$ A(1) = 0.0_RP
2273 !!$
2274 !!$ B(1) = CAP1 * DZ(1) / DELT &
2275 !!$ + 2.0_RP * AKS1 / (DZ(1)+DZ(2))
2276 !!$ C(1) = -2.0_RP * AKS1 / (DZ(1)+DZ(2))
2277 !!$ D(1) = CAP1 * DZ(1) / DELT * TSL(1) + G0
2278 !!$
2279 !!$ do K = 2, KM-1
2280 !!$ A(K) = -2.0_RP * AKS / (DZ(K-1)+DZ(K))
2281 !!$ B(K) = CAP * DZ(K) / DELT + 2.0_RP * AKS / (DZ(K-1)+DZ(K)) + 2.0_RP * AKS / (DZ(K)+DZ(K+1))
2282 !!$ C(K) = -2.0_RP * AKS / (DZ(K)+DZ(K+1))
2283 !!$ D(K) = CAP * DZ(K) / DELT * TSL(K)
2284 !!$ end do
2285 !!$
2286 !!$ if( BOUND == 1 ) then ! Flux=0
2287 !!$ A(KM) = -2.0_RP * AKS / (DZ(KM-1)+DZ(KM))
2288 !!$ B(KM) = CAP * DZ(KM) / DELT + 2.0_RP * AKS / (DZ(KM-1)+DZ(KM))
2289 !!$ C(KM) = 0.0_RP
2290 !!$ D(KM) = CAP * DZ(KM) / DELT * TSL(KM)
2291 !!$ else ! T=constant
2292 !!$ A(KM) = -2.0_RP * AKS / (DZ(KM-1)+DZ(KM))
2293 !!$ B(KM) = CAP * DZ(KM) / DELT + 2.0_RP * AKS / (DZ(KM-1)+DZ(KM)) + 2.0_RP * AKS / (DZ(KM)+DZEND)
2294 !!$ C(KM) = 0.0_RP
2295 !!$ D(KM) = CAP * DZ(KM) / DELT * TSL(KM) + 2.0_RP * AKS * TSLEND / (DZ(KM)+DZEND)
2296 !!$ end if
2297 !!$
2298 !!$ P(1) = -C(1) / B(1)
2299 !!$ Q(1) = D(1) / B(1)
2300 !!$
2301 !!$ !$acc loop seq
2302 !!$ do K = 2, KM
2303 !!$ P(K) = -C(K) / ( A(K) * P(K-1) + B(K) )
2304 !!$ Q(K) = ( -A(K) * Q(K-1) + D(K) ) / ( A(K) * P(K-1) + B(K) )
2305 !!$ end do
2306 !!$
2307 !!$ X(KM) = Q(KM)
2308 !!$
2309 !!$ !$acc loop seq
2310 !!$ do K = KM-1, 1, -1
2311 !!$ X(K) = P(K) * X(K+1) + Q(K)
2312 !!$ end do
2313 !!$
2314 !!$ do K = 1, KM
2315 !!$ TSL(K) = X(K)
2316 !!$ enddo
2317 !!$
2318 !!$ return
2319 !!$ end subroutine multi_layer2
2320 
2321  !-----------------------------------------------------------------------------
2323  subroutine urban_param_setup
2324  implicit none
2325 
2326  real(rp) :: dhgt,thgt,vfws,vfgs
2327  integer :: k
2328 
2329  ! initialize
2330  r = 0.0_rp
2331  rw = 0.0_rp
2332  hgt = 0.0_rp
2333  z0hr = 0.0_rp
2334  z0hb = 0.0_rp
2335  z0hg = 0.0_rp
2336  zdc_tbl = 0.0_rp
2337  z0c_tbl = 0.0_rp
2338  z0hc_tbl = 0.0_rp
2339  svf = 0.0_rp
2340 
2341  ! set up other urban parameters
2342  z0hr = 0.1_rp * z0r
2343  z0hb = 0.1_rp * z0b
2344  z0hg = 0.1_rp * z0g
2345  zdc_tbl = zr * 0.3_rp
2346  z0c_tbl = zr * 0.15_rp
2347  z0hc_tbl = 0.1_rp * z0c_tbl
2348 
2349  ! HGT: Normalized height
2350  hgt = zr / ( road_width + roof_width )
2351 
2352  ! R: Normalized Roof Width (a.k.a. "building coverage ratio")
2353  r = roof_width / ( road_width + roof_width )
2354  rw = 1.0_rp - r
2355 
2356  ! Calculate Sky View Factor:
2357  dhgt = hgt / 100.0_rp
2358  thgt = 0.0_rp
2359  vfws = 0.0_rp
2360  thgt = hgt - dhgt / 2.0_rp
2361  do k = 1, 99
2362  thgt = thgt - dhgt
2363  vfws = vfws + 0.25_rp * ( 1.0_rp - thgt / sqrt( thgt**2 + rw**2 ) )
2364  end do
2365 
2366  vfws = vfws / 99.0_rp
2367  vfws = vfws * 2.0_rp
2368  vfgs = 1.0_rp - 2.0_rp * vfws * hgt / rw
2369  svf = vfgs
2370 
2371  return
2372  end subroutine urban_param_setup
2373 
2374  !-----------------------------------------------------------------------------
2376  subroutine read_urban_param_table( INFILENAME )
2377  use scale_prc, only: &
2378  prc_abort
2379  implicit none
2380 
2381  character(*), intent(in) :: infilename
2382 
2383  real(rp) :: betr, betb, betg
2384 
2385  namelist / param_urban_data / &
2386  zr, &
2387  roof_width, &
2388  road_width, &
2389  sigma_zed, &
2390  ah_tbl, &
2391  ahl_tbl, &
2392  betr, &
2393  betb, &
2394  betg, &
2395  strgr, &
2396  strgb, &
2397  strgg, &
2398  capr, &
2399  capb, &
2400  capg, &
2401  aksr, &
2402  aksb, &
2403  aksg, &
2404  albr, &
2405  albb, &
2406  albg, &
2407  epsr, &
2408  epsb, &
2409  epsg, &
2410  z0r, &
2411  z0b, &
2412  z0g, &
2413  trlend, &
2414  tblend, &
2415  tglend
2416 
2417  character(len=H_LONG) :: fname
2418  integer :: fid
2419  integer :: ierr
2420  !---------------------------------------------------------------------------
2421 
2422  fid = io_get_available_fid()
2423  call io_get_fname(fname, infilename)
2424  open( fid, &
2425  file = fname, &
2426  form = 'formatted', &
2427  status = 'old', &
2428  iostat = ierr )
2429  if ( ierr /= 0 ) then
2430  log_newline
2431  log_error("URBAN_DYN_kusaka01_setup",*) 'Failed to open a parameter file : ', trim(fname)
2432  call prc_abort
2433  end if
2434 
2435  log_newline
2436  log_info("URBAN_DYN_kusaka01_setup",*) 'read_urban_param_table: Read urban parameters from file'
2437 
2438  betr = -1.0_rp
2439  betb = -1.0_rp
2440  betg = -1.0_rp
2441 
2442  !--- read namelist
2443  rewind(fid)
2444  read (fid,nml=param_urban_data,iostat=ierr)
2445  if ( ierr < 0 ) then !--- no data
2446  log_info("read_urban_param_table",*) 'Not found namelist of PARAM_URBAN_DATA. Default used.'
2447  elseif( ierr > 0 ) then !--- fatal error
2448  log_error("read_urban_param_table",*) 'Not appropriate names in PARAM_URBAN_DATA of ', trim(infilename)
2449  call prc_abort
2450  endif
2451  log_nml(param_urban_data)
2452 
2453  close( fid )
2454 
2455  if ( zr <= 0.0_rp ) then
2456  log_error("read_urban_param_table",*) 'ZR is not appropriate value; ZR must be larger than 0. ZR=', zr
2457  call prc_abort
2458  endif
2459 
2460  betr_const = betr
2461  betb_const = betb
2462  betg_const = betg
2463 
2464  return
2465  end subroutine read_urban_param_table
2466 
2467  !-----------------------------------------------------------------------------
2469  subroutine read_urban_gridded_data_2d( &
2470  UIA, UJA, &
2471  INFILENAME, &
2472  VARNAME, &
2473  udata )
2474  use scale_file_cartesc, only: &
2475  file_cartesc_open, &
2476  file_cartesc_read, &
2477  file_cartesc_flush, &
2478  file_cartesc_check_coordinates, &
2479  file_cartesc_close
2480  use scale_prc, only: &
2481  prc_abort
2482  implicit none
2483  integer , intent(in) :: uia, uja
2484  character(len=*), intent(in) :: infilename
2485  character(len=*), intent(in) :: varname
2486  real(rp), intent(out) :: udata(uia,uja)
2487 
2488  integer :: fid
2489  !---------------------------------------------------------------------------
2490  log_newline
2491  log_info("URBAN_DYN_kusaka01_setup",*) 'read_urban_gridded_data ',trim(varname)
2492 
2493  fid = io_get_available_fid()
2494  call file_cartesc_open( infilename, fid )
2495  call file_cartesc_read( fid, varname, 'XY', udata(:,:) )
2496 
2497  call file_cartesc_flush( fid )
2498  call file_cartesc_check_coordinates( fid )
2499  call file_cartesc_close( fid )
2500 
2501  return
2502  end subroutine read_urban_gridded_data_2d
2503 
2504  !-----------------------------------------------------------------------------
2506  subroutine read_urban_gridded_data_3d( &
2507  UIA, UJA, &
2508  INFILENAME, &
2509  VARNAME, &
2510  udata )
2511  use scale_file_cartesc, only: &
2512  file_cartesc_open, &
2513  file_cartesc_read, &
2514  file_cartesc_flush, &
2515  file_cartesc_check_coordinates, &
2516  file_cartesc_close
2517  use scale_prc, only: &
2518  prc_abort
2519  implicit none
2520  integer , intent(in) :: uia, uja
2521  character(len=*), intent(in) :: infilename
2522  character(len=*), intent(in) :: varname
2523  real(rp), intent(out) :: udata(uia,uja,24)
2524 
2525  integer :: fid, k
2526  !---------------------------------------------------------------------------
2527  log_newline
2528  log_info("URBAN_DYN_kusaka01_setup",*) 'read_urban_gridded_data ',trim(varname)
2529 
2530  fid = io_get_available_fid()
2531  call file_cartesc_open( infilename, fid )
2532  do k = 1, 24
2533  call file_cartesc_read( fid, varname, 'XY', udata(:,:,k), k )
2534  enddo
2535 
2536  call file_cartesc_flush( fid )
2537  call file_cartesc_check_coordinates( fid )
2538  call file_cartesc_close( fid )
2539 
2540  return
2541  end subroutine read_urban_gridded_data_3d
2542 
2543  !-----------------------------------------------------------------------------
2544  subroutine put_history( &
2545  UIA, UJA, &
2546  SHR, SHB, SHG, &
2547  LHR, LHB, LHG, &
2548  GHR, GHB, GHG, &
2549  RNR, RNB, RNG, &
2550  RNgrd )
2551  use scale_file_history, only: &
2552  file_history_put
2553  integer, intent(in) :: UIA, UJA
2554  real(RP), intent(in) :: SHR(UIA,UJA), SHB(UIA,UJA), SHG(UIA,UJA)
2555  real(RP), intent(in) :: LHR(UIA,UJA), LHB(UIA,UJA), LHG(UIA,UJA)
2556  real(RP), intent(in) :: GHR(UIA,UJA), GHB(UIA,UJA), GHG(UIA,UJA)
2557  real(RP), intent(in) :: RNR(UIA,UJA), RNB(UIA,UJA), RNG(UIA,UJA)
2558  real(RP), intent(in) :: RNgrd(UIA,UJA)
2559 
2560  call file_history_put( i_shr, shr(:,:) )
2561  call file_history_put( i_shb, shb(:,:) )
2562  call file_history_put( i_shg, shg(:,:) )
2563  call file_history_put( i_lhr, lhr(:,:) )
2564  call file_history_put( i_lhb, lhb(:,:) )
2565  call file_history_put( i_lhg, lhg(:,:) )
2566  call file_history_put( i_ghr, ghr(:,:) )
2567  call file_history_put( i_ghb, ghb(:,:) )
2568  call file_history_put( i_ghg, ghg(:,:) )
2569  call file_history_put( i_rnr, rnr(:,:) )
2570  call file_history_put( i_rnb, rnb(:,:) )
2571  call file_history_put( i_rng, rng(:,:) )
2572  call file_history_put( i_rngrd, rngrd(:,:) )
2573 
2574  return
2575  end subroutine put_history
2576 
2577 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:350
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_prc::prc_local_comm_world
integer, public prc_local_comm_world
local communicator
Definition: scale_prc.F90:89
scale_const::const_rvap
real(rp), parameter, public const_rvap
specific gas constant (water vapor) [J/kg/K]
Definition: scale_const.F90:68
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:35
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:91
scale_io::io_get_available_fid
integer function, public io_get_available_fid()
search & get available file ID
Definition: scale_io.F90:373
scale_file_history
module file_history
Definition: scale_file_history.F90:15
scale_prc
module PROCESS
Definition: scale_prc.F90:11
scale_mapprojection
module Map projection
Definition: scale_mapprojection.F90:12
scale_io::io_get_fname
subroutine, public io_get_fname(outstr, instr, rank, ext, len)
generate process specific filename
Definition: scale_io.F90:421
scale_io
module STDIO
Definition: scale_io.F90:10
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:1044
scale_const::const_cpdry
real(rp), public const_cpdry
specific heat (dry air,constant pressure) [J/kg/K]
Definition: scale_const.F90:60
scale_prof
module profiler
Definition: scale_prof.F90:11
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, AH_TOFFSET)
Setup.
Definition: scale_urban_dyn_kusaka01.F90:139
scale_precision::dp
integer, parameter, public dp
Definition: scale_precision.F90:32
scale_mapprojection::mapprojection_basepoint_lon
real(rp), public mapprojection_basepoint_lon
Definition: scale_mapprojection.F90:90
scale_const::const_pi
real(rp), parameter, public const_pi
pi
Definition: scale_const.F90:32
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:1018
scale_urban_dyn_kusaka01::urban_dyn_kusaka01_finalize
subroutine, public urban_dyn_kusaka01_finalize
Finalize.
Definition: scale_urban_dyn_kusaka01.F90:399
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:2551
scale_file_cartesc::file_cartesc_open
subroutine, public file_cartesc_open(basename, fid, single, aggregate)
open a netCDF file for read
Definition: scale_file_cartesC.F90:760
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, 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:456
scale_const::const_rdry
real(rp), public const_rdry
specific gas constant (dry air) [J/kg/K]
Definition: scale_const.F90:59
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:685
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:43
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:57
scale_file_cartesc
module file / cartesianC
Definition: scale_file_cartesC.F90:11
scale_prc::prc_ismaster
logical, public prc_ismaster
master process in local communicator?
Definition: scale_prc.F90:92
scale_bulkflux::bulkflux
procedure(bc), pointer, public bulkflux
Definition: scale_bulkflux.F90:84
Generated by   doxygen 1.8.17