| blob | 048d27dcac4bd4176f122ee13fcf1c2984956d48 |
1 MODULE module_sf_urban
2 #if defined(mpas)
3 use mpas_atmphys_utilities, only: physics_message,physics_error_fatal
4 #define FATAL_ERROR(M) call physics_error_fatal( M )
5 #define WRITE_MESSAGE(M) call physics_message( M )
6 #else
7 use module_wrf_error
8 #define FATAL_ERROR(M) call wrf_error_fatal( M )
9 #define WRITE_MESSAGE(M) call wrf_message( M )
10 #endif
12 USE module_model_constants, ONLY : piconst
14 !===============================================================================
15 ! Single-Layer Urban Canopy Model for WRF Noah-LSM
16 ! Original Version: 2002/11/06 by Hiroyuki Kusaka
17 ! Last Update: 2006/08/24 by Fei Chen and Mukul Tewari (NCAR/RAL)
18 !===============================================================================
20 CHARACTER(LEN=4) :: LU_DATA_TYPE
22 INTEGER :: ICATE
24 REAL, ALLOCATABLE, DIMENSION(:) :: ZR_TBL
25 REAL, ALLOCATABLE, DIMENSION(:) :: Z0C_TBL
26 REAL, ALLOCATABLE, DIMENSION(:) :: Z0HC_TBL
27 REAL, ALLOCATABLE, DIMENSION(:) :: ZDC_TBL
28 REAL, ALLOCATABLE, DIMENSION(:) :: SVF_TBL
29 REAL, ALLOCATABLE, DIMENSION(:) :: R_TBL
30 REAL, ALLOCATABLE, DIMENSION(:) :: RW_TBL
31 REAL, ALLOCATABLE, DIMENSION(:) :: HGT_TBL
32 REAL, ALLOCATABLE, DIMENSION(:) :: AH_TBL
33 REAL, ALLOCATABLE, DIMENSION(:) :: ALH_TBL
34 REAL, ALLOCATABLE, DIMENSION(:) :: BETR_TBL
35 REAL, ALLOCATABLE, DIMENSION(:) :: BETB_TBL
36 REAL, ALLOCATABLE, DIMENSION(:) :: BETG_TBL
37 REAL, ALLOCATABLE, DIMENSION(:) :: FRC_URB_TBL
39 REAL, ALLOCATABLE, DIMENSION(:) :: COP_TBL
40 REAL, ALLOCATABLE, DIMENSION(:) :: BLDAC_FRC_TBL
41 REAL, ALLOCATABLE, DIMENSION(:) :: COOLED_FRC_TBL
42 REAL, ALLOCATABLE, DIMENSION(:) :: PWIN_TBL
43 REAL, ALLOCATABLE, DIMENSION(:) :: BETA_TBL
44 INTEGER, ALLOCATABLE, DIMENSION(:) :: SW_COND_TBL
45 REAL, ALLOCATABLE, DIMENSION(:) :: TIME_ON_TBL
46 REAL, ALLOCATABLE, DIMENSION(:) :: TIME_OFF_TBL
47 REAL, ALLOCATABLE, DIMENSION(:) :: TARGTEMP_TBL
48 REAL, ALLOCATABLE, DIMENSION(:) :: GAPTEMP_TBL
49 REAL, ALLOCATABLE, DIMENSION(:) :: TARGHUM_TBL
50 REAL, ALLOCATABLE, DIMENSION(:) :: GAPHUM_TBL
51 REAL, ALLOCATABLE, DIMENSION(:) :: PERFLO_TBL
52 REAL, ALLOCATABLE, DIMENSION(:) :: PV_FRAC_ROOF_TBL !GRZ
53 REAL, ALLOCATABLE, DIMENSION(:) :: GR_FRAC_ROOF_TBL !GRZ
54 INTEGER :: GR_FLAG_TBL !GRZ
55 INTEGER :: GR_TYPE_TBL !GRZ
56 REAL, DIMENSION(1:24) :: IRHO_TBL
57 REAL, ALLOCATABLE, DIMENSION(:) :: HSESF_TBL
59 REAL, ALLOCATABLE, DIMENSION(:) :: CAPR_TBL, CAPB_TBL, CAPG_TBL
60 REAL, ALLOCATABLE, DIMENSION(:) :: AKSR_TBL, AKSB_TBL, AKSG_TBL
61 REAL, ALLOCATABLE, DIMENSION(:) :: ALBR_TBL, ALBB_TBL, ALBG_TBL
62 REAL, ALLOCATABLE, DIMENSION(:) :: EPSR_TBL, EPSB_TBL, EPSG_TBL
63 REAL, ALLOCATABLE, DIMENSION(:) :: Z0R_TBL, Z0B_TBL, Z0G_TBL
64 REAL, ALLOCATABLE, DIMENSION(:) :: SIGMA_ZED_TBL
65 REAL, ALLOCATABLE, DIMENSION(:) :: Z0HB_TBL, Z0HG_TBL
66 REAL, ALLOCATABLE, DIMENSION(:) :: TRLEND_TBL, TBLEND_TBL, TGLEND_TBL
67 REAL, ALLOCATABLE, DIMENSION(:) :: AKANDA_URBAN_TBL
68 !for BEP
69 ! MAXDIRS :: The maximum number of street directions we're allowed to define
70 INTEGER, PARAMETER :: MAXDIRS = 3
71 ! MAXHGTS :: The maximum number of building height bins we're allowed to define
72 INTEGER, PARAMETER :: MAXHGTS = 50
74 INTEGER, ALLOCATABLE, DIMENSION(:) :: NUMDIR_TBL
75 REAL, ALLOCATABLE, DIMENSION(:,:) :: STREET_DIRECTION_TBL
76 REAL, ALLOCATABLE, DIMENSION(:,:) :: STREET_WIDTH_TBL
77 REAL, ALLOCATABLE, DIMENSION(:,:) :: BUILDING_WIDTH_TBL
78 INTEGER, ALLOCATABLE, DIMENSION(:) :: NUMHGT_TBL
79 REAL, ALLOCATABLE, DIMENSION(:,:) :: HEIGHT_BIN_TBL
80 REAL, ALLOCATABLE, DIMENSION(:,:) :: HPERCENT_BIN_TBL
81 !end BEP
82 INTEGER :: BOUNDR_DATA,BOUNDB_DATA,BOUNDG_DATA
83 INTEGER :: CH_SCHEME_DATA, TS_SCHEME_DATA
84 INTEGER :: ahoption ! Miao, 2007/01/17, cal. ah
85 REAL, DIMENSION(1:24) :: ahdiuprf ! ah diurnal profile, tloc: 1-24
86 REAL, DIMENSION(1:24) :: hsequip_tbl
88 !===Yang, 2014/10/08, urban hydrological processes for single layer UCM===
89 INTEGER :: IMP_SCHEME, IRI_SCHEME
90 INTEGER :: alhoption ! anthropogenic latent heat option
91 INTEGER :: groption ! anthropogenic latent heat option
92 LOGICAL :: distributed_aerodynamics_option
93 REAL :: fgr ! green roof fraction
94 REAL :: oasis ! urban oasis parameter
95 REAL, DIMENSION(1:4) :: DZGR ! Layer depth of green roof
96 REAL, DIMENSION(1:4) :: alhseason ! seasonal variation of alh
97 REAL, DIMENSION(1:48) :: alhdiuprf ! alh diurnal profile, tloc2: 1-48
98 REAL, DIMENSION(1:3) :: porimp ! porosity of pavement over impervious surface
99 REAL, DIMENSION(1:3) :: dengimp ! maximum water-holding depth of pavement
101 !===end hydrological processes===
103 INTEGER :: allocate_status
105 ! INTEGER :: num_roof_layers
106 ! INTEGER :: num_wall_layers
107 ! INTEGER :: num_road_layers
109 CHARACTER (LEN=256) , PRIVATE :: mesg
111 CONTAINS
113 !===============================================================================
114 !
115 ! Author:
116 ! Hiroyuki KUSAKA, PhD
117 ! University of Tsukuba, JAPAN
118 ! (CRIEPI, NCAR/MMM visiting scientist, 2002-2004)
119 ! kusaka@ccs.tsukuba.ac.jp
120 !
121 ! Co-Researchers:
122 ! Fei CHEN, PhD
123 ! NCAR/RAP feichen@ucar.edu
124 ! Mukul TEWARI, PhD
125 ! NCAR/RAP mukul@ucar.edu
126 !
127 ! Purpose:
128 ! Calculate surface temeprature, fluxes, canopy air temperature, and canopy wind
129 !
130 ! Subroutines:
131 ! module_sf_urban
132 ! |- urban
133 ! |- read_param
134 ! |- mos or jurges
135 ! |- multi_layer or force_restore
136 ! |- urban_param_init <-- URBPARM.TBL
137 ! |- urban_var_init
138 !
139 ! Input Data from WRF [MKS unit]:
140 !
141 ! UTYPE [-] : Urban type. 1=Commercial/Industrial; 2=High-intensity residential;
142 ! : 3=low-intensity residential
143 ! TA [K] : Potential temperature at 1st wrf level (absolute temp)
144 ! QA [kg/kg] : Mixing ratio at 1st atmospheric level
145 ! UA [m/s] : Wind speed at 1st atmospheric level
146 ! SSG [W/m/m] : Short wave downward radiation at a flat surface
147 ! Note this is the total of direct and diffusive solar
148 ! downward radiation. If without two components, the
149 ! single solar downward can be used instead.
150 ! SSG = SSGD + SSGQ
151 ! LSOLAR [-] : Indicating the input type of solar downward radiation
152 ! True: both direct and diffusive solar radiation
153 ! are available
154 ! False: only total downward ridiation is available.
155 ! SSGD [W/m/m] : Direct solar radiation at a flat surface
156 ! if SSGD is not available, one can assume a ratio SRATIO
157 ! (e.g., 0.7), so that SSGD = SRATIO*SSG
158 ! SSGQ [W/m/m] : Diffuse solar radiation at a flat surface
159 ! If SSGQ is not available, SSGQ = SSG - SSGD
160 ! LLG [W/m/m] : Long wave downward radiation at a flat surface
161 ! RAIN [mm/h] : Precipitation
162 ! RHOO [kg/m/m/m] : Air density
163 ! ZA [m] : First atmospheric level
164 ! as a lowest boundary condition
165 ! DECLIN [rad] : solar declination
166 ! COSZ : = sin(fai)*sin(del)+cos(fai)*cos(del)*cos(omg)
167 ! OMG [rad] : solar hour angle
168 ! XLAT [deg] : latitude
169 ! DELT [sec] : Time step
170 ! ZNT [m] : Roughnes length
171 !
172 ! Output Data to WRF [MKS unit]:
173 !
174 ! TS [K] : Surface potential temperature (absolute temp)
175 ! QS [-] : Surface humidity
176 !
177 ! SH [W/m/m/] : Sensible heat flux, = FLXTH*RHOO*CPP
178 ! LH [W/m/m] : Latent heat flux, = FLXHUM*RHOO*ELL
179 ! LH_INEMATIC [kg/m/m/sec]: Moisture Kinematic flux, = FLXHUM*RHOO
180 ! SW [W/m/m] : Upward shortwave radiation flux,
181 ! = SSG-SNET*697.7*60. (697.7*60.=100.*100.*4.186)
182 ! ALB [-] : Time-varying albedo
183 ! LW [W/m/m] : Upward longwave radiation flux,
184 ! = LNET*697.7*60.-LLG
185 ! G [W/m/m] : Heat Flux into the Ground
186 ! RN [W/m/m] : Net radiation
187 !
188 ! PSIM [-] : Diagnostic similarity stability function for momentum
189 ! PSIH [-] : Diagnostic similarity stability function for heat
190 !
191 ! TC [K] : Diagnostic canopy air temperature
192 ! QC [-] : Diagnostic canopy humidity
193 !
194 ! TH2 [K] : Diagnostic potential temperature at 2 m
195 ! Q2 [-] : Diagnostic humidity at 2 m
196 ! U10 [m/s] : Diagnostic u wind component at 10 m
197 ! V10 [m/s] : Diagnostic v wind component at 10 m
198 !
199 ! CHS, CHS2 [m/s] : CH*U at ZA, CH*U at 2 m (not used)
200 !
201 ! Important parameters:
202 !
203 ! Morphology of the urban canyon:
204 ! These parameters assigned in the URBPARM.TBL
205 !
206 ! ZR [m] : roof level (building height)
207 ! SIGMA_ZED [m] : Standard Deviation of roof height
208 ! ROOF_WIDTH [m] : roof (i.e., building) width
209 ! ROAD_WIDTH [m] : road width
210 !
211 ! Parameters derived from the morphological terms above.
212 ! These parameters are computed in the code.
213 !
214 ! HGT [-] : normalized building height
215 ! SVF [-] : sky view factor
216 ! R [-] : Normalized roof width (a.k.a. "building coverage ratio")
217 ! RW [-] : = 1 - R
218 ! Z0C [m] : Roughness length above canyon for momentum (1/10 of ZR)
219 ! Z0HC [m] : Roughness length above canyon for heat (1/10 of Z0C)
220 ! ZDC [m] : Zero plane displacement height (1/5 of ZR)
221 !
222 ! Following parameter are assigned in run/URBPARM.TBL
223 !
224 ! AH [ W m{-2} ] : anthropogenic heat ( W m{-2} in the table, converted internally to cal cm{-2} )
225 ! ALH [ W m{-2} ] : anthropogenic latent heat ( W m{-2} in the table, converted internally to cal cm{-2} )
226 ! CAPR[ J m{-3} K{-1} ] : heat capacity of roof ( units converted in code to [ cal cm{-3} deg{-1} ] )
227 ! CAPB[ J m{-3} K{-1} ] : heat capacity of building wall ( units converted in code to [ cal cm{-3} deg{-1} ] )
228 ! CAPG[ J m{-3} K{-1} ] : heat capacity of road ( units converted in code to [ cal cm{-3} deg{-1} ] )
229 ! AKSR [ J m{-1} s{-1} K{-1} ] : thermal conductivity of roof ( units converted in code to [ cal cm{-1} s{-1} deg{-1} ] )
230 ! AKSB [ J m{-1} s{-1} K{-1} ] : thermal conductivity of building wall ( units converted in code to [ cal cm{-1} s{-1} deg{-1} ] )
231 ! AKSG [ J m{-1} s{-1} K{-1} ] : thermal conductivity of road ( units converted in code to [ cal cm{-1} s{-1} deg{-1} ] )
232 ! ALBR [-] : surface albedo of roof
233 ! ALBB [-] : surface albedo of building wall
234 ! ALBG [-] : surface albedo of road
235 ! EPSR [-] : surface emissivity of roof
236 ! EPSB [-] : surface emissivity of building wall
237 ! EPSG [-] : surface emissivity of road
238 ! Z0B [m] : roughness length for momentum of building wall (only for CH_SCHEME = 1)
239 ! Z0G [m] : roughness length for momentum of road (only for CH_SCHEME = 1)
240 ! Z0HB [m] : roughness length for heat of building wall (only for CH_SCHEME = 1)
241 ! Z0HG [m] : roughness length for heat of road
242 ! num_roof_layers : number of layers within roof
243 ! num_wall_layers : number of layers within building walls
244 ! num_road_layers : number of layers within below road surface
245 ! NOTE: for now, these layers are defined as same as the number of soil layers in namelist.input
246 ! DZR [cm] : thickness of each roof layer
247 ! DZB [cm] : thickness of each building wall layer
248 ! DZG [cm] : thickness of each ground layer
249 ! BOUNDR [integer 1 or 2] : Boundary Condition for Roof Layer Temp [1: Zero-Flux, 2: T = Constant]
250 ! BOUNDB [integer 1 or 2] : Boundary Condition for Building Wall Layer Temp [1: Zero-Flux, 2: T = Constant]
251 ! BOUNDG [integer 1 or 2] : Boundary Condition for Road Layer Temp [1: Zero-Flux, 2: T = Constant]
252 ! TRLEND [K] : lower boundary condition of roof temperature
253 ! TBLEND [K] : lower boundary condition of building temperature
254 ! TGLEND [K] : lower boundary condition of ground temperature
255 ! CH_SCHEME [integer 1 or 2] : Sfc exchange scheme used for building wall and road
256 ! [1: M-O Similarity Theory, 2: Empirical Form (recommend)]
257 ! TS_SCHEME [integer 1 or 2] : Scheme for computing surface temperature (for roof, wall, and road)
258 ! [1: 4-layer model, 2: Force-Restore method]
259 ! IMP_SCHEME[integer 1 or 2] : Evaporation scheme for impervious surfaces (roof, wall, and road)
260 ! [1: Hypothesized evaporation during large rainfall events
261 ! [2: Water-holding scheme over impervious surface
262 ! IRI_SCHEME[integer 0 or 1] : Scheme for urban irrigation
263 ! [0: No irrigation, 1: Summertime (May-Sep) irrigation everyday at 9pm]
264 !for BEP
265 ! numdir [ - ] : Number of street directions defined for a particular urban category
266 ! street_direction [ deg ] : Direction of streets for a particular urban category and a particular street direction
267 ! street_width [ m ] : Width of street for a particular urban category and a particular street direction
268 ! building_width [ m ] : Width of buildings for a particular urban category and a particular street direction
269 ! numhgt [ - ] : Number of building height levels defined for a particular urban category
270 ! height_bin [ m ] : Building height bins defined for a particular urban category.
271 ! hpercent_bin [ % ] : Percentage of a particular urban category populated by buildings of particular height bins
272 !end BEP
273 ! Moved from URBPARM.TBL
274 !
275 ! BETR [-] : minimum moisture availability of roof
276 ! BETB [-] : minimum moisture availability of building wall
277 ! BETG [-] : minimum moisture availability of road
278 ! Z0R [m] : roughness length for momentum of roof
279 ! Z0HB [m] : roughness length for heat of building wall (only for CH_SCHEME = 1)
280 ! Z0HG [m] : roughness length for heat of road
281 ! num_roof_layers : number of layers within roof
282 ! num_wall_layers : number of layers within building walls
283 ! num_road_layers : number of layers within below road surface
284 ! NOTE: for now, these layers are defined as same as the number of soil layers in namelist.input
285 !
286 ! References:
287 ! Kusaka and Kimura (2004) J.Appl.Meteor., vol.43, p1899-1910
288 ! Kusaka and Kimura (2004) J.Meteor.Soc.Japan, vol.82, p45-65
289 ! Kusaka et al. (2001) Bound.-Layer Meteor., vol.101, p329-358
290 !
291 ! History:
292 ! 2014/10, modified by Jiachuan Yang (ASU)
293 ! 2006/06 modified by H. Kusaka (Univ. Tsukuba), M. Tewari
294 ! 2005/10/26, modified by Fei Chen, Mukul Tewari
295 ! 2003/07/21 WRF , modified by H. Kusaka of CRIEPI (NCAR/MMM)
296 ! 2001/08/26 PhD , modified by H. Kusaka of CRIEPI (Univ.Tsukuba)
297 ! 1999/08/25 LCM , developed by H. Kusaka of CRIEPI (Univ.Tsukuba)
298 !
299 !===============================================================================
300 !
301 ! subroutine urban:
302 !
303 !===============================================================================
305 SUBROUTINE urban(LSOLAR, & ! L
306 num_roof_layers,num_wall_layers,num_road_layers, & ! I
307 DZR,DZB,DZG, & ! I
308 UTYPE,TA,QA,UA,U1,V1,SSG,SSGD,SSGQ,LLG,RAIN,RHOO, & ! I
309 ZA,DECLIN,COSZ,OMG,XLAT,DELT,ZNT, & ! I
310 CHS, CHS2, & ! I
311 TR, TB, TG, TC, QC, UC, & ! H
312 TRL,TBL,TGL, & ! H
313 XXXR, XXXB, XXXG, XXXC, & ! H
314 TS,QS,SH,LH,LH_KINEMATIC, & ! O
315 SW,ALB,LW,G,RN,PSIM,PSIH, & ! O
316 GZ1OZ0, & ! O
317 CMR_URB,CHR_URB,CMC_URB,CHC_URB, & ! I/O
318 U10,V10,TH2,Q2,UST,mh_urb,stdh_urb,lf_urb, & ! O
319 lp_urb,hgt_urb,frc_urb,lb_urb,zo_check, & ! O
320 CMCR,TGR,TGRL,SMR,CMGR_URB,CHGR_URB,jmonth, & ! H
321 DRELR,DRELB,DRELG,FLXHUMR,FLXHUMB,FLXHUMG, &
322 lf_urb_s, z0_urb, vegfrac_in)
324 IMPLICIT NONE
326 REAL, PARAMETER :: CP=0.24 ! heat capacity of dry air [cgs unit]
327 REAL, PARAMETER :: EL=583. ! latent heat of vaporation [cgs unit]
328 REAL, PARAMETER :: SIG=8.17E-11 ! stefun bolzman constant [cgs unit]
329 REAL, PARAMETER :: SIG_SI=5.67E-8 ! [MKS unit]
330 REAL, PARAMETER :: AK=0.4 ! kalman const. [-]
331 REAL, PARAMETER :: PI=3.14159 ! pi [-]
332 REAL, PARAMETER :: TETENA=7.5 ! const. of Tetens Equation [-]
333 REAL, PARAMETER :: TETENB=237.3 ! const. of Tetens Equation [-]
334 REAL, PARAMETER :: SRATIO=0.75 ! ratio between direct/total solar [-]
336 REAL, PARAMETER :: CPP=1004.5 ! heat capacity of dry air [J/K/kg]
337 REAL, PARAMETER :: ELL=2.442E+06 ! latent heat of vaporization [J/kg]
338 REAL, PARAMETER :: XKA=2.4E-5
340 !-------------------------------------------------------------------------------
341 ! C: configuration variables
342 !-------------------------------------------------------------------------------
344 LOGICAL, INTENT(IN) :: LSOLAR ! logical [true=both, false=SSG only]
346 ! The following variables are also model configuration variables, but are
347 ! defined in the URBAN.TBL and in the contains statement in the top of
348 ! the module_urban_init, so we should not declare them here.
350 INTEGER, INTENT(IN) :: num_roof_layers
351 INTEGER, INTENT(IN) :: num_wall_layers
352 INTEGER, INTENT(IN) :: num_road_layers
355 REAL, INTENT(IN), DIMENSION(1:num_roof_layers) :: DZR ! grid interval of roof layers [cm]
356 REAL, INTENT(IN), DIMENSION(1:num_wall_layers) :: DZB ! grid interval of wall layers [cm]
357 REAL, INTENT(IN), DIMENSION(1:num_road_layers) :: DZG ! grid interval of road layers [cm]
359 !-------------------------------------------------------------------------------
360 ! I: input variables from LSM to Urban
361 !-------------------------------------------------------------------------------
363 INTEGER, INTENT(IN) :: UTYPE ! urban type [1=Commercial/Industrial, 2=High-intensity residential,
364 ! 3=low-intensity residential]
365 INTEGER, INTENT(IN) :: jmonth! current month
366 REAL, INTENT(IN) :: TA ! potential temp at 1st atmospheric level [K]
367 REAL, INTENT(IN) :: QA ! mixing ratio at 1st atmospheric level [kg/kg]
368 REAL, INTENT(IN) :: UA ! wind speed at 1st atmospheric level [m/s]
369 REAL, INTENT(IN) :: U1 ! u at 1st atmospheric level [m/s]
370 REAL, INTENT(IN) :: V1 ! v at 1st atmospheric level [m/s]
371 REAL, INTENT(IN) :: SSG ! downward total short wave radiation [W/m/m]
372 REAL, INTENT(IN) :: LLG ! downward long wave radiation [W/m/m]
373 REAL, INTENT(IN) :: RAIN ! precipitation [mm/h]
374 REAL, INTENT(IN) :: RHOO ! air density [kg/m^3]
375 REAL, INTENT(IN) :: ZA ! first atmospheric level [m]
376 REAL, INTENT(IN) :: DECLIN ! solar declination [rad]
377 REAL, INTENT(IN) :: COSZ ! sin(fai)*sin(del)+cos(fai)*cos(del)*cos(omg)
378 REAL, INTENT(IN) :: OMG ! solar hour angle [rad]
380 REAL, INTENT(IN) :: XLAT ! latitude [deg]
381 REAL, INTENT(IN) :: DELT ! time step [s]
382 REAL, INTENT(IN) :: CHS,CHS2 ! CH*U at za and 2 m [m/s]
384 REAL, INTENT(INOUT) :: SSGD ! downward direct short wave radiation [W/m/m]
385 REAL, INTENT(INOUT) :: SSGQ ! downward diffuse short wave radiation [W/m/m]
386 REAL, INTENT(INOUT) :: CMR_URB
387 REAL, INTENT(INOUT) :: CHR_URB
388 REAL, INTENT(INOUT) :: CMC_URB
389 REAL, INTENT(INOUT) :: CHC_URB
390 REAL, INTENT(INOUT) :: ZNT ! roughness length [m] ! modified by danli
391 !-------------------------------------------------------------------------------
392 ! I: NUDAPT Input Parameters
393 !-------------------------------------------------------------------------------
394 REAL, INTENT(INOUT) :: mh_urb ! mean building height [m]
395 REAL, INTENT(INOUT) :: stdh_urb ! standard deviation of building height [m]
396 REAL, INTENT(INOUT) :: hgt_urb ! area weighted mean building height [m]
397 REAL, INTENT(INOUT) :: lp_urb ! plan area fraction [-]
398 REAL, INTENT(INOUT) :: frc_urb ! urban fraction [-]
399 REAL, INTENT(INOUT) :: lb_urb ! building surface to plan area ratio [-]
400 REAL, INTENT(INOUT), DIMENSION(4) :: lf_urb ! frontal area index [-]
401 REAL, INTENT(INOUT) :: zo_check ! check for printing ZOC
403 !-------------------------------------------------------------------------------
404 ! I: Distributed aerodynamics parameters
405 !-------------------------------------------------------------------------------
406 REAL, INTENT(IN) :: lf_urb_s ! frontal area index [-]
407 REAL, INTENT(IN) :: z0_urb ! roughness length [m]
408 REAL, INTENT(IN) :: vegfrac_in ! vegetation fraction (0 to 1) [-]
410 !-------------------------------------------------------------------------------
411 ! O: output variables from Urban to LSM
412 !-------------------------------------------------------------------------------
414 REAL, INTENT(OUT) :: TS ! surface potential temperature [K]
415 REAL, INTENT(OUT) :: QS ! surface humidity [K]
416 REAL, INTENT(OUT) :: SH ! sensible heat flux [W/m/m]
417 REAL, INTENT(OUT) :: LH ! latent heat flux [W/m/m]
418 REAL, INTENT(OUT) :: LH_KINEMATIC ! latent heat, kinetic [kg/m/m/s]
419 REAL, INTENT(OUT) :: SW ! upward short wave radiation flux [W/m/m]
420 REAL, INTENT(OUT) :: ALB ! time-varying albedo [fraction]
421 REAL, INTENT(OUT) :: LW ! upward long wave radiation flux [W/m/m]
422 REAL, INTENT(OUT) :: G ! heat flux into the ground [W/m/m]
423 REAL, INTENT(OUT) :: RN ! net radition [W/m/m]
424 REAL, INTENT(OUT) :: PSIM ! similality stability shear function for momentum
425 REAL, INTENT(OUT) :: PSIH ! similality stability shear function for heat
426 REAL, INTENT(OUT) :: GZ1OZ0
427 REAL, INTENT(OUT) :: U10 ! u at 10m [m/s]
428 REAL, INTENT(OUT) :: V10 ! u at 10m [m/s]
429 REAL, INTENT(OUT) :: TH2 ! potential temperature at 2 m [K]
430 REAL, INTENT(OUT) :: Q2 ! humidity at 2 m [-]
431 !m REAL, INTENT(OUT) :: CHS,CHS2 ! CH*U at za and 2 m [m/s]
432 REAL, INTENT(OUT) :: UST ! friction velocity [m/s]
435 !-------------------------------------------------------------------------------
436 ! H: Historical (state) variables of Urban : LSM <--> Urban
437 !-------------------------------------------------------------------------------
439 ! TR: roof temperature [K]; TRP: at previous time step [K]
440 ! TB: building wall temperature [K]; TBP: at previous time step [K]
441 ! TG: road temperature [K]; TGP: at previous time step [K]
442 ! TC: urban-canopy air temperature [K]; TCP: at previous time step [K]
443 ! (absolute temperature)
444 ! QC: urban-canopy air mixing ratio [kg/kg]; QCP: at previous time step [kg/kg]
445 !
446 ! XXXR: Monin-Obkhov length for roof [dimensionless]
447 ! XXXB: Monin-Obkhov length for building wall [dimensionless]
448 ! XXXG: Monin-Obkhov length for road [dimensionless]
449 ! XXXC: Monin-Obkhov length for urban-canopy [dimensionless]
450 !
451 ! TRL, TBL, TGL: layer temperature [K] (absolute temperature)
453 REAL, INTENT(INOUT):: TR, TB, TG, TC, QC, UC
454 REAL, INTENT(INOUT):: XXXR, XXXB, XXXG, XXXC
456 REAL, DIMENSION(1:num_roof_layers), INTENT(INOUT) :: TRL
457 REAL, DIMENSION(1:num_wall_layers), INTENT(INOUT) :: TBL
458 REAL, DIMENSION(1:num_road_layers), INTENT(INOUT) :: TGL
460 !===Yang,2014/10/08, urban hydrological variables for single layer UCM===
461 ! FLXHUMR: evaporation over roof [m/s]; FLXHUMRP: at previous time step [m/s]
462 ! FLXHUMB: evaporation over wall [m/s]; FLXHUMBP: at previous time step [m/s]
463 ! FLXHUMG: evaporation over road [m/s]; FLXHUMGP: at previous time step [m/s]
465 ! DRELR: water retention depth on roof [m]; DRELRP: at previous time stp [m]
466 ! DRELB: water retention depth on wall [m]; DRELBP: at previous time stp [m]
467 ! DRELG: water retention depth on road [m]; DRELGP: at previous time stp [m]
469 ! TGR: green roof surface temperature [K]; TGRP: at previous time step [K]
470 ! CMCR: Canopy intercepted water on green roof; CMCRP: at previous time step
471 ! SMR: soil moisture at each layer on roof [-]; SMRP: at previous time step
472 ! TGRL:layer temperature on green roof [K]
474 REAL, INTENT(INOUT):: FLXHUMR,FLXHUMB,FLXHUMG,DRELR,DRELB,DRELG
475 REAL, INTENT(INOUT):: TGR,CMCR,CHGR_URB,CMGR_URB
476 REAL, DIMENSION(1:num_roof_layers), INTENT(INOUT) :: SMR
477 REAL, DIMENSION(1:num_roof_layers), INTENT(INOUT) :: TGRL
479 !-------------------------------------------------------------------------------
480 ! L: Local variables from read_param
481 !-------------------------------------------------------------------------------
483 REAL :: ZR, Z0C, Z0HC, ZDC, SVF, R, RW, HGT, AH, ALH
484 REAL :: SIGMA_ZED
485 REAL :: CAPR, CAPB, CAPG, AKSR, AKSB, AKSG, ALBR, ALBB, ALBG
486 REAL :: EPSR, EPSB, EPSG, Z0R, Z0B, Z0G, Z0HB, Z0HG
487 REAL :: TRLEND,TBLEND,TGLEND
488 REAL :: T1VR, T1VC,TH2V
489 REAL :: RLMO_URB
490 REAL :: AKANDA_URBAN
492 REAL :: TH2X !m
494 INTEGER :: BOUNDR, BOUNDB, BOUNDG
495 INTEGER :: CH_SCHEME, TS_SCHEME
497 LOGICAL :: SHADOW ! [true=consider svf and shadow effects, false=consider svf effect only]
499 !for BEP
500 INTEGER :: NUMDIR
501 REAL, DIMENSION ( MAXDIRS ) :: STREET_DIRECTION
502 REAL, DIMENSION ( MAXDIRS ) :: STREET_WIDTH
503 REAL, DIMENSION ( MAXDIRS ) :: BUILDING_WIDTH
504 INTEGER :: NUMHGT
505 REAL, DIMENSION ( MAXHGTS ) :: HEIGHT_BIN
506 REAL, DIMENSION ( MAXHGTS ) :: HPERCENT_BIN
507 !end BEP
508 !-------------------------------------------------------------------------------
509 ! L: Local variables
510 !-------------------------------------------------------------------------------
512 REAL :: BETR, BETB, BETG
513 REAL :: SX, SD, SQ, RX
514 REAL :: UR, ZC, XLB, BB
515 REAL :: Z, RIBB, RIBG, RIBC, BHR, BHB, BHG, BHC
516 REAL :: TSC, LNET, SNET, FLXUV, THG, FLXTH, FLXHUM, FLXG
517 REAL :: W, VFGS, VFGW, VFWG, VFWS, VFWW
518 REAL :: HOUI1, HOUI2, HOUI3, HOUI4, HOUI5, HOUI6, HOUI7, HOUI8
519 REAL :: SLX, SLX1, SLX2, SLX3, SLX4, SLX5, SLX6, SLX7, SLX8
520 REAL :: FLXTHR, FLXTHB, FLXTHG
521 REAL :: SR, SB, SG, RR, RB, RG
522 REAL :: SR1, SR2, SB1, SB2, SG1, SG2, RR1, RR2, RB1, RB2, RG1, RG2
523 REAL :: HR, HB, HG, ELER, ELEB, ELEG, G0R, G0B, G0G
524 REAL :: ALPHAC, ALPHAR, ALPHAB, ALPHAG
525 REAL :: CHC, CHR, CHB, CHG, CDC, CDR, CDB, CDG, CDGR
526 REAL :: C1R, C1B, C1G, TE, TC1, TC2, QC1, QC2, QS0R, QS0B, QS0G,RHO,ES
528 REAL :: DESDT
529 REAL :: F
530 REAL :: DQS0RDTR
531 REAL :: DRRDTR, DHRDTR, DELERDTR, DG0RDTR
532 REAL :: DTR, DFDT
533 REAL :: FX, FY, GF, GX, GY
534 REAL :: DTCDTB, DTCDTG
535 REAL :: DQCDTB, DQCDTG
536 REAL :: DRBDTB1, DRBDTG1, DRBDTB2, DRBDTG2
537 REAL :: DRGDTB1, DRGDTG1, DRGDTB2, DRGDTG2
538 REAL :: DRBDTB, DRBDTG, DRGDTB, DRGDTG
539 REAL :: DHBDTB, DHBDTG, DHGDTB, DHGDTG
540 REAL :: DELEBDTB, DELEBDTG, DELEGDTG, DELEGDTB
541 REAL :: DG0BDTB, DG0BDTG, DG0GDTG, DG0GDTB
542 REAL :: DQS0BDTB, DQS0GDTG
543 REAL :: DTB, DTG, DTC
545 REAL :: THEATAZ ! Solar Zenith Angle [rad]
546 REAL :: THEATAS ! = PI/2. - THETAZ
547 REAL :: FAI ! Latitude [rad]
548 REAL :: CNT,SNT
549 REAL :: PS ! Surface Pressure [hPa]
550 REAL :: TAV ! Vertial Temperature [K]
552 REAL :: XXX, X, Z0, Z0H, CD, CH
553 REAL :: XXX2, PSIM2, PSIH2, XXX10, PSIM10, PSIH10
554 REAL :: PSIX, PSIT, PSIX2, PSIT2, PSIX10, PSIT10
556 REAL :: TRP, TBP, TGP, TCP, QCP, TST, QST
557 REAL :: TSP, CHS_LOCAL, CHS2_LOCAL
559 REAL :: WDR,HGT2,BW,DHGT
560 REAL, parameter :: VonK = 0.4
561 REAL :: lambda_f,alpha_macd,beta_macd,lambda_fr
562 REAL :: lambda_p, vegfrac
564 INTEGER :: iteration, K, NUDAPT
565 INTEGER :: tloc, tloc2, Kalh
567 !===Yang,2014/10/08, urban hydrological variables for single layer UCM===
568 REAL :: FLXHUMRP, FLXHUMBP, FLXHUMGP
569 REAL :: DRELRP, DRELBP, DRELGP
570 REAL :: TGRP, CMCRP
571 REAL, DIMENSION(1:num_roof_layers) :: ZSOILR, ETR, SMRP
572 !===Define parameters for green roof===
573 INTEGER :: KZ
574 REAL :: RUNOFF1, RUNOFF2, RUNOFF3
575 REAL :: SGR, SGR1, T1VGR, CHGR, ALPHAGR
576 REAL :: FLXTHGR, FLXHUMGR, HGR, ELEGR, G0GR
577 REAL :: QS0GR, EPGR, EDIR, ETTR, FV, DTGR, DRIP
578 ! REAL :: DQS0GRDTGR, ETR, ECR,RAIN1, RAINDR, DEW, ETAR, BETGR
579 REAL :: DQS0GRDTGR, ECR,RAIN1, RAINDR, DEW, ETAR, BETGR
580 ! REAL :: DF1, RGR, RGRR, RCH, RR1, RR2, YY, ZZ1, SSOILR
581 REAL :: DF1, RGR, RGRR, RCH, YY, ZZ1, SSOILR
582 REAL :: DRRDTGR, DHRDTGR, DELERDTGR, DG0RDTGR, DFDVT
583 real,parameter :: SHDFAC = 0.80 ! Vegetated area fraction of green roof vegetation
584 real,parameter :: ALBV = 0.20 ! green roof albedo
585 real,parameter :: EPSV = 0.93 ! green roof emissivity
586 real,parameter :: LAI = 1.50 ! leaf area index on green roof
587 real,parameter :: CMCMAX = 0.5E-3 ! Maximum canopy interception capacity
588 real,parameter :: SMCREF = 0.329 ! Reference soil moisture
589 real,parameter :: SMCDRY = 0.066 ! Residual soil moisture
590 real,parameter :: SMCWLT = 0.084 ! Wilting point
591 real,parameter :: SMCMAX = 0.439 ! Saturated soil moisture
592 real,parameter :: RSMAX = 5000 ! Maximum stomatal resistance
593 real,parameter :: RSMIN = 100 ! Minimum stomatal resistance
594 real,parameter :: RGL = 100 ! Radiation limit where photosynthesis begins
595 real,parameter :: CFACTR = 0.5 ! Parameter used in the canopy inteception calculation
596 real,parameter :: DWSAT = 0.143E-4 ! Saturated soil conductivity
597 real,parameter :: DKSAT = 3.38E-6 ! Saturated soil diffusivity
598 real,parameter :: BEXP = 5.25 ! B parameter in soil hydraulic calculation
599 real,parameter :: FXEXP = 2.0 ! Parameter for computing direct soil evaporation
600 real,parameter :: ZBOT = -2.0
601 real,parameter :: QUARTZ = 0.40
602 real,parameter :: CSOIL = 2.0E+6
603 real,parameter :: HS = 36
604 integer,parameter :: NROOT = 2 ! Root depth layer of green roof
605 integer,parameter :: NGR = 4 ! Layer of green roof
606 integer,parameter :: IMPR = 1
607 integer,parameter :: IMPB = 2
608 integer,parameter :: IMPG = 3
610 SHADOW = .false.
611 ! SHADOW = .true.
613 IF (distributed_aerodynamics_option .and. groption == 1) THEN
614 FATAL_ERROR("slucm_distributed_drag is not compatible with groption")
615 END IF
618 !-------------------------------------------------------------------------------
619 ! Set parameters
620 !-------------------------------------------------------------------------------
622 ! Miao, 2007/01/17, cal. ah
623 if(ahoption==1) then
624 tloc=mod(int(OMG/PI*180./15.+12.+0.5 ),24)
625 if(tloc.lt.0) tloc=tloc+24
626 if(tloc==0) tloc=24
627 endif
628 ! Yang, 2014/10/08, cal. alh
629 if(alhoption==1) then
630 tloc2=mod(int((OMG/PI*180./15.+12.)*2.+0.5 ),48)
631 if(tloc2.lt.0) tloc2=tloc2+48
632 if(tloc2==0) tloc2=48
633 endif
635 CALL read_param(UTYPE,ZR,SIGMA_ZED,Z0C,Z0HC,ZDC,SVF,R,RW,HGT, &
636 AH,CAPR,CAPB,CAPG,AKSR,AKSB,AKSG,ALBR,ALBB, &
637 ALBG,EPSR,EPSB,EPSG,Z0R,Z0B,Z0G,Z0HB,Z0HG, &
638 BETR,BETB,BETG,TRLEND,TBLEND,TGLEND, &
639 !for BEP
640 NUMDIR, STREET_DIRECTION, STREET_WIDTH, &
641 BUILDING_WIDTH, NUMHGT, HEIGHT_BIN, &
642 HPERCENT_BIN, &
643 !end BEP
644 BOUNDR,BOUNDB,BOUNDG,CH_SCHEME,TS_SCHEME, &
645 AKANDA_URBAN,ALH)
647 ! Glotfelty, 2012/07/05, NUDAPT Modification
649 if (mh_urb.gt.0.0 .and. .not. distributed_aerodynamics_option) THEN
650 !write(mesg,*) 'Mean Height NUDAPT',mh_urb
651 !WRITE_MESSAGE(mesg)
652 !write(mesg,*) 'Mean Height Table',ZR
653 !WRITE_MESSAGE(mesg)
654 if(zo_check.eq.1)THEN
655 write(mesg,*) 'Mean Height NUDAPT',mh_urb
656 WRITE_MESSAGE(mesg)
657 write(mesg,*) 'Mean Height Table',ZR
658 WRITE_MESSAGE(mesg)
659 write(mesg,*) 'Roughness Length Table',Z0C
660 WRITE_MESSAGE(mesg)
661 write(mesg,*) 'Roof Roughness Length Table',Z0R
662 WRITE_MESSAGE(mesg)
663 write(mesg,*) 'Sky View Factor Table',SVF
664 WRITE_MESSAGE(mesg)
665 write(mesg,*) 'Normalized Height Table',HGT
666 WRITE_MESSAGE(mesg)
667 write(mesg,*) 'Plan Area Fraction', lp_urb
668 WRITE_MESSAGE(mesg)
669 write(mesg,*) 'Plan Area Fraction table', R
670 WRITE_MESSAGE(mesg)
671 end if
672 !write(mesg,*) 'Area Weighted Mean Height',hgt_urb
673 !WRITE_MESSAGE(mesg)
674 !write(mesg,*) 'Plan Area Fraction', lp_urb
675 !WRITE_MESSAGE(mesg)
676 !write(mesg,*) 'STD Height', stdh_urb
677 !WRITE_MESSAGE(mesg)
678 !write(mesg,*) 'Frontal Area Index',lf_urb
679 !WRITE_MESSAGE(mesg)
680 !write(mesg,*) 'Urban Fraction',frc_urb
681 !WRITE_MESSAGE(mesg)
682 !write(mesg,*) 'Building Surf Ratio',lb_urb
683 !WRITE_MESSAGE(mesg)
685 !Calculate Building Width and Street Width Based on BEP formulation
686 if(lb_urb.gt.lp_urb)THEN
687 BW=2.*hgt_urb*lp_urb/(lb_urb-lp_urb)
688 SW=2.*hgt_urb*lp_urb*((frc_urb/lp_urb)-1.)/(lb_urb-lp_urb)
689 !write(mesg,*) 'Building Width',BW
690 !WRITE_MESSAGE(mesg)
691 !write(mesg,*) 'Street Width',SW
692 !WRITE_MESSAGE(mesg)
693 elseif (SW.lt.0.0.or.BW.lt.0.0)then
694 BW=BUILDING_WIDTH(1)
695 SW=STREET_WIDTH(1)
696 else
697 BW=BUILDING_WIDTH(1)
698 SW=STREET_WIDTH(1)
699 end if
701 !Assign NUDAPT Parameters
702 ZR = mh_urb
703 R = lp_urb
704 RW = 1.0-R
705 HGT = mh_urb/(BW+SW)
706 SIGMA_ZED = stdh_urb
708 !Calculate Wind Direction and Assign Appropriae lf_urb
709 WDR = (180.0/PI)*ATAN2(U10,V10)
711 IF(WDR.ge.0.0.and.WDR.lt.22.5)THEN
712 lambda_f = lf_urb(1)
713 ELSEIF(WDR.ge.-22.5.and.WDR.lt.0.0)THEN
714 lambda_f = lf_urb(1)
715 ELSEIF(WDR.gt.157.5.and.WDR.le.180.0)THEN
716 lambda_f = lf_urb(1)
717 ELSEIF(WDR.lt.-157.5)THEN
718 lambda_f = lf_urb(1)
719 ELSEIF(WDR.gt.22.5.and.WDR.le.67.5)THEN
720 lambda_f = lf_urb(2)
721 ELSEIF(WDR.ge.-157.5.and.WDR.lt.-112.5)THEN
722 lambda_f = lf_urb(2)
723 ELSEIF(WDR.gt.67.5.and.WDR.le.112.5)THEN
724 lambda_f = lf_urb(3)
725 ELSEIF(WDR.ge.-112.5.and.WDR.lt.-67.5)THEN
726 lambda_f = lf_urb(3)
727 ELSEIF(WDR.gt.112.5.and.WDR.le.157.5)THEN
728 lambda_f = lf_urb(4)
729 ELSEIF(WDR.ge.-67.5.and.WDR.lt.-22.5)THEN
730 lambda_f = lf_urb(4)
731 ELSE
732 lambda_f = lf_urb(1)
733 ENDIF
735 !Calculate the following urban canyon geometry parameters following Macdonald's (1998) formulations
736 Cd = 1.2
737 alpha_macd = 4.43
738 beta_macd = 1.0
741 ZDC = ZR * ( 1.0 + ( alpha_macd ** ( -R ) ) * ( R - 1.0 ) )
743 Z0C = ZR * ( 1.0 - ZDC/ZR ) * &
744 exp (-(0.5 * beta_macd * Cd / (VonK**2) * ( 1.0-ZDC/ZR) * lambda_f )**(-0.5))
746 if(zo_check.eq.1)THEN
747 write(mesg,*) 'Roughness Length NUDAPT',Z0C
748 WRITE_MESSAGE(mesg)
749 end if
751 lambda_fr = stdh_urb/(SW + BW)
753 Z0R = ZR * ( 1.0 - ZDC/ZR) &
754 * exp ( -(0.5 * beta_macd * Cd / (VonK**2) &
755 * ( 1.0-ZDC/ZR) * lambda_fr )**(-0.5))
759 Z0HC = 0.1 * Z0C
761 ! Calculate Sky View Factor
763 DHGT=HGT/100.
764 HGT2=0.
765 VFWS=0.
766 HGT2=HGT-DHGT/2.
767 do NUDAPT=1,99
768 HGT2=HGT2-DHGT
769 VFWS=VFWS+0.25*(1.-HGT2/SQRT(HGT2**2.+RW**2.))
770 end do
772 VFWS=VFWS/99.
773 VFWS=VFWS*2.
775 VFGS=1.-2.*VFWS*HGT/RW
776 SVF=VFGS
778 if(zo_check.eq.1)THEN
779 write(mesg,*) 'Roof Roughness Length NUDAPT',Z0R
780 WRITE_MESSAGE(mesg)
781 write(mesg,*) 'Sky View Factor NUDAPT',SVF
782 WRITE_MESSAGE(mesg)
783 write(mesg,*) 'normalized Height NUDAPT', HGT
784 WRITE_MESSAGE(mesg)
785 end if
788 endif
790 !End NUDAPT Modification
793 ! Miao, 2007/01/17, cal. ah
794 if(ahoption==1) AH=AH*ahdiuprf(tloc)
796 ! Yang, 2014/10/08, cal. alh
797 Kalh=0
798 if(alhoption==1) THEN
799 if(jmonth==3 .or. jmonth==4 .or. jmonth==5) Kalh=1
800 if(jmonth==6 .or. jmonth==7 .or. jmonth==8) Kalh=2
801 if(jmonth==9 .or. jmonth==10.or. jmonth==11)Kalh=3
802 if(jmonth==12.or. jmonth==1 .or. jmonth==2) Kalh=4
803 endif
804 if(alhoption==1) ALH = ALH*alhdiuprf(tloc2)*alhseason(Kalh)
806 IF (distributed_aerodynamics_option) THEN
807 ZDC = 0.
808 IF (Z0_URB > MH_URB) THEN
809 FATAL_ERROR("Z0_URB is larger than MH_URB")
810 END IF
811 ZR = MAX(MH_URB, 3.5)
812 Z0C = MAX(Z0_URB, 0.1)
813 lambda_p = MAX(0.05, MIN(1.0, LP_URB))
814 lambda_f = MAX(0.05, MIN(1.0, LF_URB_S))
816 R = lambda_p
817 RW = 1 - R
818 SVF = kanda_kawai_svf(lambda_p, lambda_f)
820 vegfrac = MIN(0.9, MAX(0.1, vegfrac_in))
822 HGT = lambda_f
823 END IF
825 IF( ZDC+Z0C+2. >= ZA) THEN
826 FATAL_ERROR("ZDC + Z0C + 2m is larger than the 1st WRF level - Stop in subroutine urban - change ZDC and Z0C" )
827 END IF
829 IF(.NOT.LSOLAR) THEN
830 SSGD = SRATIO*SSG
831 SSGQ = SSG - SSGD
832 ENDIF
833 SSGD = SRATIO*SSG ! No radiation scheme has SSGD and SSGQ.
834 SSGQ = SSG - SSGD
836 W=2.*1.*HGT
837 VFGS=SVF
838 VFGW=1.-SVF
839 VFWG=(1.-SVF)*(1.-R)/W
840 VFWS=VFWG
841 VFWW=1.-2.*VFWG
843 !-------------------------------------------------------------------------------
844 ! Convert unit from MKS to cgs
845 ! Renew surface and layer temperatures
846 !-------------------------------------------------------------------------------
848 SX=(SSGD+SSGQ)/697.7/60. ! downward short wave radition [ly/min]
849 SD=SSGD/697.7/60. ! downward direct short wave radiation
850 SQ=SSGQ/697.7/60. ! downward diffuse short wave radiation
851 RX=LLG/697.7/60. ! downward long wave radiation
852 RHO=RHOO*0.001 ! air density at first atmospheric level
854 TRP=TR
855 TBP=TB
856 TGP=TG
857 TCP=TC
858 QCP=QC
860 TSP = (TR * R + TB * W + TG * RW) / (R + RW + W)
862 !===Yang,2014/10/08, urban hydrological variables for single layer UCM===
863 FLXHUMRP = FLXHUMR
864 FLXHUMBP = FLXHUMB
865 FLXHUMGP = FLXHUMG
866 DRELRP = DRELR
867 DRELBP = DRELB
868 DRELGP = DRELG
869 TGRP = TGR
870 CMCRP = CMCR
871 SMRP = SMR
873 !===Yang,2014/10/08, urban irrigation, May-Sep, 9-10pm
874 IF(IRI_SCHEME==1) THEN
875 IF (tloc==21 .or. tloc==22) THEN
876 IF(jmonth==5 .or. jmonth==6 .or. jmonth ==7 .or. &
877 jmonth==8 .or. jmonth==9 ) THEN
878 DO KZ = 1,2
879 SMRP(KZ)= SMCREF
880 END DO
881 ENDIF
882 ENDIF
883 ENDIF
885 TAV=TA*(1.+0.61*QA)
886 PS=RHOO*287.*TAV/100. ![hPa]
888 !-------------------------------------------------------------------------------
889 ! Canopy wind
890 !-------------------------------------------------------------------------------
892 IF ( ZR + 2. < ZA ) THEN
893 UR=UA*LOG((ZR-ZDC)/Z0C)/LOG((ZA-ZDC)/Z0C)
894 ZC=0.7*ZR
895 XLB=0.4*(ZR-ZDC)
896 ! BB formulation from Inoue (1963)
897 BB = 0.4 * ZR / ( XLB * alog( ( ZR - ZDC ) / Z0C ) )
898 UC=UR*EXP(-BB*(1.-ZC/ZR))
899 ELSE
900 ! PRINT *, 'Warning ZR + 2m is larger than the 1st WRF level'
901 ZC=ZA/2.
902 UC=UA/2.
903 END IF
905 !-------------------------------------------------------------------------------
906 ! Net Short Wave Radiation at roof, wall, and road
907 !-------------------------------------------------------------------------------
909 IF (SSG > 0.0) THEN
911 IF(.NOT.SHADOW) THEN ! no shadow effects model
913 SR1=SX*(1.-ALBR)
914 SGR1=SX*(1.-ALBV)
915 SG1=SX*VFGS*(1.-ALBG)
916 SB1=SX*VFWS*(1.-ALBB)
917 SG2=SB1*ALBB/(1.-ALBB)*VFGW*(1.-ALBG)
918 SB2=SG1*ALBG/(1.-ALBG)*VFWG*(1.-ALBB) + SB1*ALBB*VFWW
920 ELSE ! shadow effects model
922 FAI=XLAT*PI/180.
924 THEATAS=ABS(ASIN(COSZ))
925 THEATAZ=ABS(ACOS(COSZ))
927 SNT=COS(DECLIN)*SIN(OMG)/COS(THEATAS)
928 CNT=(COSZ*SIN(FAI)-SIN(DECLIN))/COS(THEATAS)/COS(FAI)
930 HOUI1=(SNT*COS(PI/8.) -CNT*SIN(PI/8.))
931 HOUI2=(SNT*COS(2.*PI/8.) -CNT*SIN(2.*PI/8.))
932 HOUI3=(SNT*COS(3.*PI/8.) -CNT*SIN(3.*PI/8.))
933 HOUI4=(SNT*COS(4.*PI/8.) -CNT*SIN(4.*PI/8.))
934 HOUI5=(SNT*COS(5.*PI/8.) -CNT*SIN(5.*PI/8.))
935 HOUI6=(SNT*COS(6.*PI/8.) -CNT*SIN(6.*PI/8.))
936 HOUI7=(SNT*COS(7.*PI/8.) -CNT*SIN(7.*PI/8.))
937 HOUI8=(SNT*COS(8.*PI/8.) -CNT*SIN(8.*PI/8.))
939 SLX1=HGT*ABS(TAN(THEATAZ))*ABS(HOUI1)
940 SLX2=HGT*ABS(TAN(THEATAZ))*ABS(HOUI2)
941 SLX3=HGT*ABS(TAN(THEATAZ))*ABS(HOUI3)
942 SLX4=HGT*ABS(TAN(THEATAZ))*ABS(HOUI4)
943 SLX5=HGT*ABS(TAN(THEATAZ))*ABS(HOUI5)
944 SLX6=HGT*ABS(TAN(THEATAZ))*ABS(HOUI6)
945 SLX7=HGT*ABS(TAN(THEATAZ))*ABS(HOUI7)
946 SLX8=HGT*ABS(TAN(THEATAZ))*ABS(HOUI8)
948 IF(SLX1 > RW) SLX1=RW
949 IF(SLX2 > RW) SLX2=RW
950 IF(SLX3 > RW) SLX3=RW
951 IF(SLX4 > RW) SLX4=RW
952 IF(SLX5 > RW) SLX5=RW
953 IF(SLX6 > RW) SLX6=RW
954 IF(SLX7 > RW) SLX7=RW
955 IF(SLX8 > RW) SLX8=RW
957 SLX=(SLX1+SLX2+SLX3+SLX4+SLX5+SLX6+SLX7+SLX8)/8.
959 SR1=SD*(1.-ALBR)+SQ*(1.-ALBR)
960 SGR1=SD*(1.-ALBV)+SQ*(1.-ALBV)
961 SG1=SD*(RW-SLX)/RW*(1.-ALBG)+SQ*VFGS*(1.-ALBG)
962 SB1=SD*SLX/W*(1.-ALBB)+SQ*VFWS*(1.-ALBB)
963 SG2=SB1*ALBB/(1.-ALBB)*VFGW*(1.-ALBG)
964 SB2=SG1*ALBG/(1.-ALBG)*VFWG*(1.-ALBB) + SB1*ALBB*VFWW
966 END IF
968 SR=SR1
969 SGR=SGR1
970 SG=SG1+SG2
971 SB=SB1+SB2
973 IF (GROPTION ==1) THEN
974 SNET=R*FGR*SGR+R*(1.-FGR)*SR+W*SB+RW*SG
975 ELSE
976 SNET=R*SR+W*SB+RW*SG
977 ENDIF
979 ELSE
981 SR=0.
982 SG=0.
983 SGR=0.
984 SB=0.
985 SNET=0.
987 END IF
989 !-------------------------------------------------------------------------------
990 ! Roof
991 !-------------------------------------------------------------------------------
993 !-------------------------------------------------------------------------------
994 ! CHR, CDR, BETR
995 !-------------------------------------------------------------------------------
997 ! Z=ZA-ZDC
998 ! BHR=LOG(Z0R/Z0HR)/0.4
999 ! RIBR=(9.8*2./(TA+TRP))*(TA-TRP)*(Z+Z0R)/(UA*UA)
1000 ! CALL mos(XXXR,ALPHAR,CDR,BHR,RIBR,Z,Z0R,UA,TA,TRP,RHO)
1002 ! Alternative option for MOST using SFCDIF routine from Noah
1003 ! Virtual temperatures needed by SFCDIF
1004 T1VR = TRP* (1.0+ 0.61 * QA)
1005 TH2V = (TA + ( 0.0098 * ZA)) * (1.0+ 0.61 * QA)
1007 ! note that CHR_URB contains UA (=CHR_MOS*UA)
1008 RLMO_URB=0.0
1009 IF (distributed_aerodynamics_option) THEN
1010 T1VC = TSP* (1.0+ 0.61 * QA)
1011 CALL SFCDIF_URB (ZA,Z0C,T1VC,TH2V,UA,AKANDA_URBAN,CMC_URB,CHC_URB,RLMO_URB,CDC,Z0HC,vegfrac)
1012 CHC = CHC_URB / UA ! canopy bulk transfer coef.
1013 ALPHAC = RHO * CP * CHC_URB
1014 CHR = CHC * R / (R + W + RW) ! local bulk transfer coef for roof
1015 CHB = CHC * W / (R + W + RW) ! local bulk transfer coef for building wall
1016 CHG = CHC * RW / (R + W + RW) ! local bulk transfer coef for floor
1017 ALPHAR = RHO * CP * CHR * UA
1018 ALPHAB = RHO * CP * CHB * UA
1019 ALPHAG = RHO * CP * CHG * UA
1020 ELSE
1021 CALL SFCDIF_URB (ZA,Z0R,T1VR,TH2V,UA,AKANDA_URBAN,CMR_URB,CHR_URB,RLMO_URB,CDR)
1022 ALPHAR = RHO*CP*CHR_URB
1023 CHR=ALPHAR/RHO/CP/UA
1024 END IF
1026 ! Yang, 03/12/2014 -- LH for impervious roof surface
1027 RAIN1 = RAIN * 0.001 /3600 ! CONVERT FROM mm/hr to m/s
1028 IF (IMP_SCHEME==1) then
1029 IF (RAIN > 1.) BETR=0.7
1030 ENDIF
1032 IF (IMP_SCHEME==2) then
1033 IF (FLXHUMRP <= 0.) FLXHUMRP = 0.
1034 ! Compute water retention depth from previous time step
1035 DrelR = DrelRP+(RAIN1-FLXHUMRP)*DELT/porimp(IMPR)
1036 IF (RAIN > 0. .AND. DrelR < DrelRP) DrelR = DrelRP
1038 IF (DrelR <= 0.) then
1039 DrelR = 0.0
1040 BETR = 0.0
1041 ELSEIf (DrelR <= dengimp(IMPR)) then
1042 BETR = DrelR/dengimp(IMPR)*porimp(IMPR)
1043 ELSE
1044 DrelR = dengimp(IMPR)
1045 BETR = porimp(IMPR)
1046 ENDIF
1048 IF ( BETR < 1.E-5 ) BETR = 0.0
1049 ENDIF
1051 IF (TS_SCHEME == 1) THEN
1053 !-------------------------------------------------------------------------------
1054 ! TR Solving Non-Linear Equation by Newton-Rapson
1055 ! TRL Solving Heat Equation by Tri Diagonal Matrix Algorithm
1056 !-------------------------------------------------------------------------------
1057 ! TSC=TRP-273.15
1058 ! ES=EXP(19.482-4303.4/(TSC+243.5)) ! WMO
1059 ! ES=6.11*10.**(TETENA*TSC/(TETENB+TSC)) ! Tetens
1060 ! DESDT=( 6.1078*(2500.-2.4*TSC)/ & ! Tetens
1061 ! (0.46151*(TSC+273.15)**2.) )*10.**(7.5*TSC/(237.3+TSC))
1062 ! ES=6.11*EXP((2.5*10.**6./461.51)*(TRP-273.15)/(273.15*TRP) ) ! Clausius-Clapeyron
1063 ! DESDT=(2.5*10.**6./461.51)*ES/(TRP**2.) ! Clausius-Clapeyron
1064 ! QS0R=0.622*ES/(PS-0.378*ES)
1065 ! DQS0RDTR = DESDT*0.622*PS/((PS-0.378*ES)**2.)
1066 ! DQS0RDTR = 17.269*(273.15-35.86)/((TRP-35.86)**2.)*QS0R
1068 ! TRP=350.
1070 DO ITERATION=1,20
1072 ES=6.11*EXP( (2.5*10.**6./461.51)*(TRP-273.15)/(273.15*TRP) )
1073 DESDT=(2.5*10.**6./461.51)*ES/(TRP**2.)
1074 QS0R=0.622*ES/(PS-0.378*ES)
1075 DQS0RDTR = DESDT*0.622*PS/((PS-0.378*ES)**2.)
1077 RR=EPSR*(RX-SIG*(TRP**4.)/60.)
1078 HR=RHO*CP*CHR*UA*(TRP-TA)*100.
1079 ELER=RHO*EL*CHR*UA*BETR*(QS0R-QA)*100.
1080 G0R=AKSR*(TRP-TRL(1))/(DZR(1)/2.)
1082 F = SR + RR - HR - ELER - G0R
1084 DRRDTR = (-4.*EPSR*SIG*TRP**3.)/60.
1085 DHRDTR = RHO*CP*CHR*UA*100.
1086 DELERDTR = RHO*EL*CHR*UA*BETR*DQS0RDTR*100.
1087 DG0RDTR = 2.*AKSR/DZR(1)
1089 DFDT = DRRDTR - DHRDTR - DELERDTR - DG0RDTR
1090 DTR = F/DFDT
1092 TR = TRP - DTR
1093 TRP = TR
1095 IF( ABS(F) < 0.000001 .AND. ABS(DTR) < 0.000001 ) EXIT
1097 END DO
1099 ! multi-layer heat equation model
1101 CALL multi_layer(num_roof_layers,BOUNDR,G0R,CAPR,AKSR,TRL,DZR,DELT,TRLEND)
1103 ELSE
1105 ES=6.11*EXP( (2.5*10.**6./461.51)*(TRP-273.15)/(273.15*TRP) )
1106 QS0R=0.622*ES/(PS-0.378*ES)
1108 RR=EPSR*(RX-SIG*(TRP**4.)/60.)
1109 HR=RHO*CP*CHR*UA*(TRP-TA)*100.
1110 ELER=RHO*EL*CHR*UA*BETR*(QS0R-QA)*100.
1111 G0R=SR+RR-HR-ELER
1113 CALL force_restore(CAPR,AKSR,DELT,SR,RR,HR,ELER,TRLEND,TRP,TR)
1115 TRP=TR
1117 END IF
1119 FLXTHR=HR/RHO/CP/100.
1120 FLXHUMR=ELER/RHO/EL/100.
1122 !-------------------------------------------------------------------------------
1123 ! Green Roof
1124 ! Must use multiple layers scheme (TS_SCHEME=1)
1125 !-------------------------------------------------------------------------------
1126 IF (GROPTION == 1) THEN
1127 T1VGR = TGRP* (1.0+ 0.61 * QA)
1128 RLMO_URB=0.0
1129 CALL SFCDIF_URB (ZA,Z0R,T1VGR,TH2V,UA,AKANDA_URBAN,CMGR_URB,CHGR_URB,RLMO_URB,CDGR)
1130 ALPHAGR = RHO*CP*CHGR_URB
1131 CHGR=ALPHAGR/RHO/CP/UA
1132 RUNOFF1 = 0.0
1133 RUNOFF2 = 0.0
1134 RUNOFF3 = 0.0
1136 KZ = 1
1137 ZSOILR (KZ) = - DZGR (KZ)
1138 DO KZ = 2,NGR
1139 ZSOILR (KZ) = - DZGR(KZ) + ZSOILR (KZ -1)
1140 END DO
1142 DO ITERATION=1,100
1143 KZ=1
1144 ES=6.11*EXP( (2.5*10.**6./461.51)*(TGRP-273.15)/(273.15*TGRP) )
1145 DESDT=(2.5*10.**6./461.51)*ES/(TGRP**2.)
1146 QS0GR=0.622*ES/(PS-0.378*ES)
1147 DQS0GRDTGR = DESDT*0.622*PS/((PS-0.378*ES)**2.)
1148 EPGR=RHOO*CHGR*UA*(QS0GR-QA) ! Potential evaporation [kg/m2/s]
1150 IF (EPGR > 0.0) THEN
1151 ! Direct evaporation from soil on green roof
1152 CALL DIREVAP (EDIR,EPGR,SMRP(KZ),SHDFAC,SMCMAX,SMCDRY,FXEXP)
1153 ! Evapotranspiration and canopy intercepted evaporation
1154 CALL TRANSP (ETTR,ETR,ECR,SHDFAC,EPGR,CMCRP,CFACTR,CMCMAX,LAI,RSMIN,RSMAX,RGL,SX, &
1155 TGRP,TA,QA,SMRP,SMCWLT,SMCREF,CPP,PS,CHGR,EPSV,DELT,NROOT,NGR,DZGR, &
1156 ZSOILR,HS)
1157 ! Update moisture in soil layers
1158 CALL SMFLX (SMRP,SMR,NGR,CMCRP,CMCR,DELT,RAIN,ZSOILR,SMCMAX,BEXP,SMCWLT,DKSAT,&
1159 DWSAT,SHDFAC,CMCMAX,RUNOFF1,RUNOFF2,RUNOFF3,EDIR,ECR,ETR,DRIP)
1160 else
1161 DEW = - EPGR
1162 RAINDR = RAIN + DEW * 3600.
1163 EDIR=0.0
1164 ECR =0.0
1165 ETTR=0.0
1166 CALL SMFLX (SMRP,SMR,NGR,CMCRP,CMCR,DELT,RAINDR,ZSOILR,SMCMAX,BEXP,SMCWLT,DKSAT,&
1167 DWSAT,SHDFAC,CMCMAX,RUNOFF1,RUNOFF2,RUNOFF3,EDIR,ECR,ETR,DRIP)
1168 END IF
1169 ! ----------------------------------------------------------------------
1170 ! CONVERT MODELED EVAPOTRANSPIRATION FROM M S-1 TO KG M-2 S-1.
1171 ! ----------------------------------------------------------------------
1172 EDIR = EDIR * 1000.0
1173 ETTR = ETTR * 1000.0
1174 ECR = ECR * 1000.0
1175 ETAR = EDIR + ETTR + ECR
1176 IF (ETAR < 1.E-20) ETAR = 0.0
1178 IF ( EPGR <= 0.0 ) THEN
1179 BETGR = 0.0
1180 ELSE
1181 BETGR = ETAR / EPGR
1182 END IF
1183 ELEGR= ETAR* RHO * EL /RHOO * 100
1185 CALL TDFCND (DF1,SMR(KZ), QUARTZ, SMCMAX )
1186 DF1 = DF1 * EXP(-2.0 * SHDFAC)
1187 RGR = EPSV*(RX-SIG*(TA**4.)/60.)
1188 RGRR= (SGR+RGR) * 697.7 * 60.
1189 RCH = RHOO*CPP*CHGR
1190 RR1 = EPSV*(TA**4) * 6.48E-8 / (PS* CHGR) + 1.0
1191 IF (RAIN > 0.0) then
1192 RR2 = RR1 + RAIN / 3600 * 4.218E+3 / RCH
1193 else
1194 RR2 = RR1
1195 end if
1196 YY = TA + (RGRR / RCH - BETGR * EPGR * ELL/ RCH) / RR2
1197 ZZ1 = DF1 / (-0.5 * ZSOILR (KZ) * RCH * RR2 ) + 1.0
1200 HGR=RHO*CP*CHGR*UA*(TGRP-TA)*100.
1201 RUNOFF3 = RUNOFF3/ DELT
1202 RUNOFF2 = RUNOFF2+ RUNOFF3
1203 G0GR = DF1*(TGRP-TGRL(1))/(DZGR(1)/2.)/697.7/60
1205 FV = SGR + RGR - HGR - ELEGR - G0GR
1206 DRRDTGR = (-4.*EPSV*SIG*TGRP**3.)/60.
1207 DHRDTGR = RHO*CP*CHGR*UA*100.
1208 DELERDTGR = RHO*EL*CHGR*UA*BETGR*DQS0GRDTGR*100.
1209 DG0RDTGR = 2.*DF1/ DZGR(KZ) * ( 1.0 / 4.1868 ) * 1.E-4
1210 DFDVT = DRRDTGR - DHRDTGR - DELERDTGR - DG0RDTGR
1211 DTGR = FV/DFDVT/ 6
1212 TGR = TGRP - DTGR
1213 TGRP = TGR
1215 IF( ABS(FV) < 0.0001 .AND. ABS(DTGR) < 0.001 ) then
1216 EXIT
1217 ENDIF
1218 END DO
1219 ! Update temperature in soil layer
1220 CALL SHFLX (SSOILR,TGRL,SMR,SMCMAX,NGR,TGRP,DELT,YY,ZZ1,ZSOILR, &
1221 TRLEND,ZBOT,SMCWLT,DF1,QUARTZ,CSOIL,CAPR)
1222 FLXTHGR=HGR/RHO/CP/100.
1223 FLXHUMGR=ELEGR/RHO/EL/100.
1224 ELSE
1225 FLXTHGR=0.
1226 FLXHUMGR=0.
1227 ENDIF
1229 !-------------------------------------------------------------------------------
1230 ! Wall and Road
1231 !-------------------------------------------------------------------------------
1233 !-------------------------------------------------------------------------------
1234 ! CHC, CHB, CDB, BETB, CHG, CDG, BETG
1235 !-------------------------------------------------------------------------------
1237 ! Z=ZA-ZDC
1238 ! BHC=LOG(Z0C/Z0HC)/0.4
1239 ! RIBC=(9.8*2./(TA+TCP))*(TA-TCP)*(Z+Z0C)/(UA*UA)
1240 !
1241 ! CALL mos(XXXC,ALPHAC,CDC,BHC,RIBC,Z,Z0C,UA,TA,TCP,RHO)
1242 ! Virtual temperatures needed by SFCDIF routine from Noah
1244 IF (.not. distributed_aerodynamics_option) THEN
1246 T1VC = TCP* (1.0+ 0.61 * QA)
1247 RLMO_URB=0.0
1248 CALL SFCDIF_URB(ZA,Z0C,T1VC,TH2V,UA,AKANDA_URBAN,CMC_URB,CHC_URB,RLMO_URB,CDC)
1249 ALPHAC = RHO*CP*CHC_URB
1251 IF (CH_SCHEME == 1) THEN
1253 Z=ZDC
1254 BHB=LOG(Z0B/Z0HB)/0.4
1255 BHG=LOG(Z0G/Z0HG)/0.4
1256 RIBB=(9.8*2./(TCP+TBP))*(TCP-TBP)*(Z+Z0B)/(UC*UC)
1257 RIBG=(9.8*2./(TCP+TGP))*(TCP-TGP)*(Z+Z0G)/(UC*UC)
1259 CALL mos(XXXB,ALPHAB,CDB,BHB,RIBB,Z,Z0B,UC,TCP,TBP,RHO)
1260 CALL mos(XXXG,ALPHAG,CDG,BHG,RIBG,Z,Z0G,UC,TCP,TGP,RHO)
1262 ELSE
1264 ALPHAB=RHO*CP*(6.15+4.18*UC)/1200.
1265 IF(UC > 5.) ALPHAB=RHO*CP*(7.51*UC**0.78)/1200.
1266 ALPHAG=RHO*CP*(6.15+4.18*UC)/1200.
1267 IF(UC > 5.) ALPHAG=RHO*CP*(7.51*UC**0.78)/1200.
1269 END IF
1271 CHC=ALPHAC/RHO/CP/UA
1272 CHB=ALPHAB/RHO/CP/UC
1273 CHG=ALPHAG/RHO/CP/UC
1275 END IF
1277 !Yang 10/10/2013 -- LH from impervious wall and ground
1278 IF (IMP_SCHEME==1) then
1279 BETB=0.0
1280 IF(RAIN > 1.) BETG=0.7
1281 ENDIF
1283 IF (IMP_SCHEME==2) then
1284 IF (FLXHUMBP <= 0.) FLXHUMBP = 0.
1285 IF (FLXHUMGP <= 0.) FLXHUMGP = 0.
1286 ! Compute water retention from previous time step for wall and ground
1287 DrelB = DrelBP+(RAIN1-FLXHUMBP)*DELT/porimp(IMPB)
1288 IF (RAIN > 0. .AND. DrelB < DrelBP) DrelB = DrelBP
1289 DrelG = DrelGP+(RAIN1-FLXHUMGP)*DELT/porimp(IMPG)
1290 IF (RAIN > 0. .AND. DrelG < DrelGP) DrelG = DrelGP
1292 IF (DrelB <= 0.) then
1293 DrelB = 0.0
1294 BETB = 0.0
1295 ELSEIf (DrelB <= dengimp(IMPB)) then
1296 BETB = DrelB/dengimp(IMPB)*porimp(IMPB)
1297 ELSE
1298 DrelB = dengimp(IMPB)
1299 BETB = porimp(IMPB)
1300 ENDIF
1302 IF (DrelG <= 0.) then
1303 DrelG = 0.0
1304 BETG = 0.0
1305 ELSEIf (DrelG <= dengimp(IMPG)) then
1306 BETG = DrelG/dengimp(IMPG)*porimp(IMPG)
1307 ELSE
1308 DrelG = dengimp(IMPG)
1309 BETG = porimp(IMPG)
1310 ENDIF
1312 if ( BETG < 1.E-5 ) BETG = 0.0
1313 if ( BETB < 1.E-5 ) BETB = 0.0
1315 ENDIF
1317 IF (TS_SCHEME == 1) THEN
1319 !-------------------------------------------------------------------------------
1320 ! TB, TG Solving Non-Linear Simultaneous Equation by Newton-Rapson
1321 ! TBL,TGL Solving Heat Equation by Tri Diagonal Matrix Algorithm
1322 !-------------------------------------------------------------------------------
1324 ! TBP=350.
1325 ! TGP=350.
1327 DO ITERATION=1,20
1329 ES=6.11*EXP( (2.5*10.**6./461.51)*(TBP-273.15)/(273.15*TBP) )
1330 DESDT=(2.5*10.**6./461.51)*ES/(TBP**2.)
1331 QS0B=0.622*ES/(PS-0.378*ES)
1332 DQS0BDTB=DESDT*0.622*PS/((PS-0.378*ES)**2.)
1334 ES=6.11*EXP( (2.5*10.**6./461.51)*(TGP-273.15)/(273.15*TGP) )
1335 DESDT=(2.5*10.**6./461.51)*ES/(TGP**2.)
1336 QS0G=0.622*ES/(PS-0.378*ES)
1337 DQS0GDTG=DESDT*0.622*PS/((PS-0.378*ES)**2.)
1339 RG1=EPSG*( RX*VFGS &
1340 +EPSB*VFGW*SIG*TBP**4./60. &
1341 -SIG*TGP**4./60. )
1343 RB1=EPSB*( RX*VFWS &
1344 +EPSG*VFWG*SIG*TGP**4./60. &
1345 +EPSB*VFWW*SIG*TBP**4./60. &
1346 -SIG*TBP**4./60. )
1348 RG2=EPSG*( (1.-EPSB)*(1.-SVF)*VFWS*RX &
1349 +(1.-EPSB)*(1.-SVF)*VFWG*EPSG*SIG*TGP**4./60. &
1350 +EPSB*(1.-EPSB)*(1.-SVF)*(1.-2.*VFWS)*SIG*TBP**4./60. )
1352 RB2=EPSB*( (1.-EPSG)*VFWG*VFGS*RX &
1353 +(1.-EPSG)*EPSB*VFGW*VFWG*SIG*(TBP**4.)/60. &
1354 +(1.-EPSB)*VFWS*(1.-2.*VFWS)*RX &
1355 +(1.-EPSB)*VFWG*(1.-2.*VFWS)*SIG*EPSG*TGP**4./60. &
1356 +EPSB*(1.-EPSB)*(1.-2.*VFWS)*(1.-2.*VFWS)*SIG*TBP**4./60. )
1358 RG=RG1+RG2
1359 RB=RB1+RB2
1361 DRBDTB1=EPSB*(4.*EPSB*SIG*TB**3.*VFWW-4.*SIG*TB**3.)/60.
1362 DRBDTG1=EPSB*(4.*EPSG*SIG*TG**3.*VFWG)/60.
1363 DRBDTB2=EPSB*(4.*(1.-EPSG)*EPSB*SIG*TB**3.*VFGW*VFWG &
1364 +4.*EPSB*(1.-EPSB)*SIG*TB**3.*VFWW*VFWW)/60.
1365 DRBDTG2=EPSB*(4.*(1.-EPSB)*EPSG*SIG*TG**3.*VFWG*VFWW)/60.
1367 DRGDTB1=EPSG*(4.*EPSB*SIG*TB**3.*VFGW)/60.
1368 DRGDTG1=EPSG*(-4.*SIG*TG**3.)/60.
1369 DRGDTB2=EPSG*(4.*EPSB*(1.-EPSB)*SIG*TB**3.*VFWW*VFGW)/60.
1370 DRGDTG2=EPSG*(4.*(1.-EPSB)*EPSG*SIG*TG**3.*VFWG*VFGW)/60.
1372 DRBDTB=DRBDTB1+DRBDTB2
1373 DRBDTG=DRBDTG1+DRBDTG2
1374 DRGDTB=DRGDTB1+DRGDTB2
1375 DRGDTG=DRGDTG1+DRGDTG2
1377 IF (distributed_aerodynamics_option) THEN
1378 HB=RHO*CP*CHB*UA*(TBP-TA)*100.
1379 HG=RHO*CP*CHG*UA*(TGP-TA)*100.
1381 DHBDTB=RHO*CP*CHB*UA*100.
1382 DHBDTG=0.
1383 DHGDTG=RHO*CP*CHG*UA*100.
1384 DHGDTB=0.
1386 ELEB=RHO*EL*CHB*UA*BETB*(QS0B-QA)*100.
1387 ELEG=RHO*EL*CHG*UA*BETG*(QS0G-QA)*100.
1389 DELEBDTB=RHO*EL*CHB*UA*BETB*DQS0BDTB*100.
1390 DELEBDTG=0.
1391 DELEGDTG=RHO*EL*CHG*UA*BETG*DQS0GDTG*100.
1392 DELEGDTB=0.
1393 ELSE
1394 HB=RHO*CP*CHB*UC*(TBP-TCP)*100.
1395 HG=RHO*CP*CHG*UC*(TGP-TCP)*100.
1397 DTCDTB=W*ALPHAB/(RW*ALPHAC+RW*ALPHAG+W*ALPHAB)
1398 DTCDTG=RW*ALPHAG/(RW*ALPHAC+RW*ALPHAG+W*ALPHAB)
1400 DHBDTB=RHO*CP*CHB*UC*(1.-DTCDTB)*100.
1401 DHBDTG=RHO*CP*CHB*UC*(0.-DTCDTG)*100.
1402 DHGDTG=RHO*CP*CHG*UC*(1.-DTCDTG)*100.
1403 DHGDTB=RHO*CP*CHG*UC*(0.-DTCDTB)*100.
1405 ELEB=RHO*EL*CHB*UC*BETB*(QS0B-QCP)*100.
1406 ELEG=RHO*EL*CHG*UC*BETG*(QS0G-QCP)*100.
1408 DQCDTB=W*ALPHAB*BETB*DQS0BDTB/(RW*ALPHAC+RW*ALPHAG*BETG+W*ALPHAB*BETB)
1409 DQCDTG=RW*ALPHAG*BETG*DQS0GDTG/(RW*ALPHAC+RW*ALPHAG*BETG+W*ALPHAB*BETB)
1411 DELEBDTB=RHO*EL*CHB*UC*BETB*(DQS0BDTB-DQCDTB)*100.
1412 DELEBDTG=RHO*EL*CHB*UC*BETB*(0.-DQCDTG)*100.
1413 DELEGDTG=RHO*EL*CHG*UC*BETG*(DQS0GDTG-DQCDTG)*100.
1414 DELEGDTB=RHO*EL*CHG*UC*BETG*(0.-DQCDTB)*100.
1415 ENDIF
1417 G0B=AKSB*(TBP-TBL(1))/(DZB(1)/2.)
1418 G0G=AKSG*(TGP-TGL(1))/(DZG(1)/2.)
1420 DG0BDTB=2.*AKSB/DZB(1)
1421 DG0BDTG=0.
1422 DG0GDTG=2.*AKSG/DZG(1)
1423 DG0GDTB=0.
1425 F = SB + RB - HB - ELEB - G0B
1426 FX = DRBDTB - DHBDTB - DELEBDTB - DG0BDTB
1427 FY = DRBDTG - DHBDTG - DELEBDTG - DG0BDTG
1429 GF = SG + RG - HG - ELEG - G0G
1430 GX = DRGDTB - DHGDTB - DELEGDTB - DG0GDTB
1431 GY = DRGDTG - DHGDTG - DELEGDTG - DG0GDTG
1433 DTB = (GF*FY-F*GY)/(FX*GY-GX*FY)
1434 DTG = -(GF+GX*DTB)/GY
1436 TB = TBP + DTB
1437 TG = TGP + DTG
1439 TBP = TB
1440 TGP = TG
1442 IF (distributed_aerodynamics_option) THEN
1443 DTC = 0.0
1444 ELSE
1445 TC1=RW*ALPHAC+RW*ALPHAG+W*ALPHAB
1446 TC2=RW*ALPHAC*TA+RW*ALPHAG*TGP+W*ALPHAB*TBP
1447 TC=TC2/TC1
1449 QC1=RW*ALPHAC+RW*ALPHAG*BETG+W*ALPHAB*BETB
1450 QC2=RW*ALPHAC*QA+RW*ALPHAG*BETG*QS0G+W*ALPHAB*BETB*QS0B
1451 QC=QC2/QC1
1453 DTC=TCP - TC
1454 TCP=TC
1455 QCP=QC
1456 END IF
1458 IF( ABS(F) < 0.000001 .AND. ABS(DTB) < 0.000001 &
1459 .AND. ABS(GF) < 0.000001 .AND. ABS(DTG) < 0.000001 &
1460 .AND. ABS(DTC) < 0.000001) EXIT
1462 END DO
1464 CALL multi_layer(num_wall_layers,BOUNDB,G0B,CAPB,AKSB,TBL,DZB,DELT,TBLEND)
1466 CALL multi_layer(num_road_layers,BOUNDG,G0G,CAPG,AKSG,TGL,DZG,DELT,TGLEND)
1468 ELSE
1470 !-------------------------------------------------------------------------------
1471 ! TB, TG by Force-Restore Method
1472 !-------------------------------------------------------------------------------
1474 ES=6.11*EXP((2.5*10.**6./461.51)*(TBP-273.15)/(273.15*TBP) )
1475 QS0B=0.622*ES/(PS-0.378*ES)
1477 ES=6.11*EXP((2.5*10.**6./461.51)*(TGP-273.15)/(273.15*TGP) )
1478 QS0G=0.622*ES/(PS-0.378*ES)
1480 RG1=EPSG*( RX*VFGS &
1481 +EPSB*VFGW*SIG*TBP**4./60. &
1482 -SIG*TGP**4./60. )
1484 RB1=EPSB*( RX*VFWS &
1485 +EPSG*VFWG*SIG*TGP**4./60. &
1486 +EPSB*VFWW*SIG*TBP**4./60. &
1487 -SIG*TBP**4./60. )
1489 RG2=EPSG*( (1.-EPSB)*(1.-SVF)*VFWS*RX &
1490 +(1.-EPSB)*(1.-SVF)*VFWG*EPSG*SIG*TGP**4./60. &
1491 +EPSB*(1.-EPSB)*(1.-SVF)*(1.-2.*VFWS)*SIG*TBP**4./60. )
1493 RB2=EPSB*( (1.-EPSG)*VFWG*VFGS*RX &
1494 +(1.-EPSG)*EPSB*VFGW*VFWG*SIG*(TBP**4.)/60. &
1495 +(1.-EPSB)*VFWS*(1.-2.*VFWS)*RX &
1496 +(1.-EPSB)*VFWG*(1.-2.*VFWS)*SIG*EPSG*TGP**4./60. &
1497 +EPSB*(1.-EPSB)*(1.-2.*VFWS)*(1.-2.*VFWS)*SIG*TBP**4./60. )
1499 RG=RG1+RG2
1500 RB=RB1+RB2
1502 HB=RHO*CP*CHB*UC*(TBP-TCP)*100.
1503 ELEB=RHO*EL*CHB*UC*BETB*(QS0B-QCP)*100.
1504 G0B=SB+RB-HB-ELEB
1506 HG=RHO*CP*CHG*UC*(TGP-TCP)*100.
1507 ELEG=RHO*EL*CHG*UC*BETG*(QS0G-QCP)*100.
1508 G0G=SG+RG-HG-ELEG
1510 CALL force_restore(CAPB,AKSB,DELT,SB,RB,HB,ELEB,TBLEND,TBP,TB)
1511 CALL force_restore(CAPG,AKSG,DELT,SG,RG,HG,ELEG,TGLEND,TGP,TG)
1513 TBP=TB
1514 TGP=TG
1516 TC1=RW*ALPHAC+RW*ALPHAG+W*ALPHAB
1517 TC2=RW*ALPHAC*TA+RW*ALPHAG*TGP+W*ALPHAB*TBP
1518 TC=TC2/TC1
1520 QC1=RW*ALPHAC+RW*ALPHAG*BETG+W*ALPHAB*BETB
1521 QC2=RW*ALPHAC*QA+RW*ALPHAG*BETG*QS0G+W*ALPHAB*BETB*QS0B
1522 QC=QC2/QC1
1524 TCP=TC
1525 QCP=QC
1527 END IF
1530 FLXTHB=HB/RHO/CP/100.
1531 FLXHUMB=ELEB/RHO/EL/100.
1532 FLXTHG=HG/RHO/CP/100.
1533 FLXHUMG=ELEG/RHO/EL/100.
1535 !-------------------------------------------------------------------------------
1536 ! Total Fluxes from Urban Canopy
1537 !-------------------------------------------------------------------------------
1538 !===Yang, 2014/10/08, cal. ah. alh. green roof===
1539 if(groption==1) then
1540 if(ahoption==1) then
1541 FLXTH = ((1.-FGR)*R*FLXTHR + FGR*R*FLXTHGR + W*FLXTHB + RW*FLXTHG)+ AH/RHOO/CPP
1542 else
1543 FLXTH = ((1.-FGR)*R*FLXTHR + FGR*R*FLXTHGR + W*FLXTHB + RW*FLXTHG)
1544 endif
1545 if(alhoption==1) then
1546 FLXHUM = ((1.-FGR)*R*FLXHUMR + FGR*R*FLXHUMGR + W*FLXHUMB + RW*FLXHUMG)+ ALH/RHOO/ELL
1547 else
1548 FLXHUM = ((1.-FGR)*R*FLXHUMR + FGR*R*FLXHUMGR + W*FLXHUMB + RW*FLXHUMG)
1549 endif
1550 FLXUV = ((1.-FGR)*R*CDR + FGR*R*CDGR + RW*CDC )*UA*UA
1551 FLXG = ((1.-FGR)*R*G0R + FGR*R*G0GR+ W*G0B + RW*G0G)
1552 LNET = (1.-FGR) * R * RR + FGR *R* RGR + W * RB + RW * RG
1553 else
1554 if(ahoption==1) then
1555 FLXTH = ( R*FLXTHR + W*FLXTHB + RW*FLXTHG ) + AH/RHOO/CPP
1556 else
1557 FLXTH = ( R*FLXTHR + W*FLXTHB + RW*FLXTHG )
1558 endif
1559 if(alhoption==1) then
1560 FLXHUM = ( R*FLXHUMR + W*FLXHUMB + RW*FLXHUMG )+ ALH/RHOO/ELL
1561 else
1562 FLXHUM = ( R*FLXHUMR + W*FLXHUMB + RW*FLXHUMG )
1563 endif
1564 IF (distributed_aerodynamics_option) THEN
1565 FLXUV = CDC * UA * UA
1566 ELSE
1567 FLXUV = ( R*CDR + RW*CDC )*UA*UA
1568 END IF
1569 FLXG = ( R*G0R + W*G0B + RW*G0G )
1570 LNET = R*RR + W*RB + RW*RG
1571 endif
1573 !----------------------------------------------------------------------------
1574 ! Convert Unit: FLUXES and u* T* q* --> WRF
1575 !----------------------------------------------------------------------------
1577 SH = FLXTH * RHOO * CPP ! Sensible heat flux [W/m/m]
1578 LH = FLXHUM * RHOO * ELL ! Latent heat flux [W/m/m]
1579 LH_KINEMATIC = FLXHUM * RHOO ! Latent heat, Kinematic [kg/m/m/s]
1580 LW = LLG - (LNET*697.7*60.) ! Upward longwave radiation [W/m/m]
1581 SW = SSG - (SNET*697.7*60.) ! Upward shortwave radiation [W/m/m]
1582 ALB = 0.
1583 IF( ABS(SSG) > 0.0001) ALB = SW/SSG ! Effective albedo [-]
1584 G = -FLXG*697.7*60. ! [W/m/m]
1585 RN = (SNET+LNET)*697.7*60. ! Net radiation [W/m/m]
1587 UST = SQRT(FLXUV) ! u* [m/s]
1588 TST = -FLXTH/UST ! T* [K]
1589 QST = -FLXHUM/UST ! q* [-]
1591 !------------------------------------------------------
1592 ! diagnostic GRID AVERAGED PSIM PSIH TS QS --> WRF
1593 !------------------------------------------------------
1595 Z0 = Z0C
1596 Z0H = Z0HC
1597 Z = ZA - ZDC
1598 ZNT = Z0 ! add by Dan Li
1600 XXX = 0.4*9.81*Z*TST/TA/UST/UST
1602 IF ( XXX >= 1. ) XXX = 1.
1603 IF ( XXX <= -5. ) XXX = -5.
1605 IF ( XXX > 0 ) THEN
1606 PSIM = -5. * XXX
1607 PSIH = -5. * XXX
1608 ELSE
1609 X = (1.-16.*XXX)**0.25
1610 PSIM = 2.*ALOG((1.+X)/2.) + ALOG((1.+X*X)/2.) - 2.*ATAN(X) + PI/2.
1611 PSIH = 2.*ALOG((1.+X*X)/2.)
1612 END IF
1614 GZ1OZ0 = ALOG(Z/Z0)
1615 CD = 0.4**2./(ALOG(Z/Z0)-PSIM)**2.
1616 CHS_LOCAL = 0.4 * UST / (ALOG(Z / Z0H) - PSIH)
1617 !
1618 !m CH = 0.4**2./(ALOG(Z/Z0)-PSIM)/(ALOG(Z/Z0H)-PSIH)
1619 !m CHS = 0.4*UST/(ALOG(Z/Z0H)-PSIH)
1620 !m TS = TA + FLXTH/CH/UA ! surface potential temp (flux temp)
1621 !m QS = QA + FLXHUM/CH/UA ! surface humidity
1622 !
1623 IF (distributed_aerodynamics_option) THEN
1624 TS = TA + FLXTH / (ALPHAC / (RHO * CP)) ! surface potential temp (flux temp)
1625 QS = QA + FLXHUM / (ALPHAC / (RHO * CP)) ! surface humidity
1626 ELSE
1627 TS = TA + FLXTH/CHS ! surface potential temp (flux temp)
1628 QS = QA + FLXHUM/CHS ! surface humidity
1629 END IF
1631 !-------------------------------------------------------
1632 ! diagnostic GRID AVERAGED U10 V10 TH2 Q2 --> WRF
1633 !-------------------------------------------------------
1635 XXX2 = (2./Z)*XXX
1636 IF ( XXX2 >= 1. ) XXX2 = 1.
1637 IF ( XXX2 <= -5. ) XXX2 = -5.
1639 IF ( XXX2 > 0 ) THEN
1640 PSIM2 = -5. * XXX2
1641 PSIH2 = -5. * XXX2
1642 ELSE
1643 X = (1.-16.*XXX2)**0.25
1644 PSIM2 = 2.*ALOG((1.+X)/2.) + ALOG((1.+X*X)/2.) - 2.*ATAN(X) + 2.*ATAN(1.)
1645 PSIH2 = 2.*ALOG((1.+X*X)/2.)
1646 END IF
1647 !
1648 !m CHS2 = 0.4*UST/(ALOG(2./Z0H)-PSIH2)
1649 !
1651 XXX10 = (10./Z)*XXX
1652 IF ( XXX10 >= 1. ) XXX10 = 1.
1653 IF ( XXX10 <= -5. ) XXX10 = -5.
1655 IF ( XXX10 > 0 ) THEN
1656 PSIM10 = -5. * XXX10
1657 PSIH10 = -5. * XXX10
1658 ELSE
1659 X = (1.-16.*XXX10)**0.25
1660 PSIM10 = 2.*ALOG((1.+X)/2.) + ALOG((1.+X*X)/2.) - 2.*ATAN(X) + 2.*ATAN(1.)
1661 PSIH10 = 2.*ALOG((1.+X*X)/2.)
1662 END IF
1664 PSIX = ALOG(Z/Z0) - PSIM
1665 PSIT = ALOG(Z/Z0H) - PSIH
1667 PSIX2 = ALOG(2./Z0) - PSIM2
1668 PSIT2 = ALOG(2./Z0H) - PSIH2
1670 PSIX10 = ALOG(10./Z0) - PSIM10
1671 PSIT10 = ALOG(10./Z0H) - PSIH10
1673 U10 = U1 * (PSIX10/PSIX) ! u at 10 m [m/s]
1674 V10 = V1 * (PSIX10/PSIX) ! v at 10 m [m/s]
1676 ! TH2 = TS + (TA-TS)*(PSIT2/PSIT) ! potential temp at 2 m [K]
1677 ! TH2 = TS + (TA-TS)*(PSIT2/PSIT) ! Fei: this seems to be temp (not potential) at 2 m [K]
1678 !Fei: consistant with M-O theory
1679 IF (distributed_aerodynamics_option) THEN
1680 CHS2_LOCAL = 0.4 * UST / (ALOG(2. / Z0H) - PSIH2)
1681 TH2 = TS + (TA - TS) * (CHS_LOCAL / CHS2_LOCAL)
1682 Q2 = QS + (QA - QS) * (CHS_LOCAL / CHS2_LOCAL)
1683 ELSE
1684 TH2 = TS + (TA-TS) *(CHS/CHS2)
1685 Q2 = QS + (QA-QS)*(PSIT2/PSIT) ! humidity at 2 m [-]
1686 END IF
1688 ! TS = (LW/SIG_SI/0.88)**0.25 ! Radiative temperature [K]
1690 END SUBROUTINE urban
1691 !===============================================================================
1692 !
1693 ! mos
1694 !
1695 !===============================================================================
1696 SUBROUTINE mos(XXX,ALPHA,CD,B1,RIB,Z,Z0,UA,TA,TSF,RHO)
1698 ! XXX: z/L (requires iteration by Newton-Rapson method)
1699 ! B1: Stanton number
1700 ! PSIM: = PSIX of LSM
1701 ! PSIH: = PSIT of LSM
1703 IMPLICIT NONE
1705 REAL, PARAMETER :: CP=0.24
1706 REAL, INTENT(IN) :: B1, Z, Z0, UA, TA, TSF, RHO
1707 REAL, INTENT(OUT) :: ALPHA, CD
1708 REAL, INTENT(INOUT) :: XXX, RIB
1709 REAL :: XXX0, X, X0, FAIH, DPSIM, DPSIH
1710 REAL :: F, DF, XXXP, US, TS, AL, XKB, DD, PSIM, PSIH
1711 INTEGER :: NEWT
1712 INTEGER, PARAMETER :: NEWT_END=10
1714 IF(RIB <= -15.) RIB=-15.
1716 IF(RIB < 0.) THEN
1718 DO NEWT=1,NEWT_END
1720 IF(XXX >= 0.) XXX=-1.E-3
1722 XXX0=XXX*Z0/(Z+Z0)
1724 X=(1.-16.*XXX)**0.25
1725 X0=(1.-16.*XXX0)**0.25
1727 PSIM=ALOG((Z+Z0)/Z0) &
1728 -ALOG((X+1.)**2.*(X**2.+1.)) &
1729 +2.*ATAN(X) &
1730 +ALOG((X0+1.)**2.*(X0**2.+1.)) &
1731 -2.*ATAN(X0)
1732 FAIH=1./SQRT(1.-16.*XXX)
1733 PSIH=ALOG((Z+Z0)/Z0)+0.4*B1 &
1734 -2.*ALOG(SQRT(1.-16.*XXX)+1.) &
1735 +2.*ALOG(SQRT(1.-16.*XXX0)+1.)
1737 DPSIM=(1.-16.*XXX)**(-0.25)/XXX &
1738 -(1.-16.*XXX0)**(-0.25)/XXX
1739 DPSIH=1./SQRT(1.-16.*XXX)/XXX &
1740 -1./SQRT(1.-16.*XXX0)/XXX
1742 F=RIB*PSIM**2./PSIH-XXX
1744 DF=RIB*(2.*DPSIM*PSIM*PSIH-DPSIH*PSIM**2.) &
1745 /PSIH**2.-1.
1747 XXXP=XXX
1748 XXX=XXXP-F/DF
1749 IF(XXX <= -10.) XXX=-10.
1751 END DO
1753 ELSE IF(RIB >= 0.142857) THEN
1755 XXX=0.714
1756 PSIM=ALOG((Z+Z0)/Z0)+7.*XXX
1757 PSIH=PSIM+0.4*B1
1759 ELSE
1761 AL=ALOG((Z+Z0)/Z0)
1762 XKB=0.4*B1
1763 DD=-4.*RIB*7.*XKB*AL+(AL+XKB)**2.
1764 IF(DD <= 0.) DD=0.
1765 XXX=(AL+XKB-2.*RIB*7.*AL-SQRT(DD))/(2.*(RIB*7.**2-7.))
1766 PSIM=ALOG((Z+Z0)/Z0)+7.*MIN(XXX,0.714)
1767 PSIH=PSIM+0.4*B1
1769 END IF
1771 US=0.4*UA/PSIM ! u*
1772 IF(US <= 0.01) US=0.01
1773 TS=0.4*(TA-TSF)/PSIH ! T*
1775 CD=US*US/UA**2. ! CD
1776 ALPHA=RHO*CP*0.4*US/PSIH ! RHO*CP*CH*U
1778 RETURN
1779 END SUBROUTINE mos
1780 !===============================================================================
1781 !
1782 ! louis79
1783 !
1784 !===============================================================================
1785 SUBROUTINE louis79(ALPHA,CD,RIB,Z,Z0,UA,RHO)
1787 IMPLICIT NONE
1789 REAL, PARAMETER :: CP=0.24
1790 REAL, INTENT(IN) :: Z, Z0, UA, RHO
1791 REAL, INTENT(OUT) :: ALPHA, CD
1792 REAL, INTENT(INOUT) :: RIB
1793 REAL :: A2, XX, CH, CMB, CHB
1795 A2=(0.4/ALOG(Z/Z0))**2.
1797 IF(RIB <= -15.) RIB=-15.
1799 IF(RIB >= 0.0) THEN
1800 IF(RIB >= 0.142857) THEN
1801 XX=0.714
1802 ELSE
1803 XX=RIB*LOG(Z/Z0)/(1.-7.*RIB)
1804 END IF
1805 CH=0.16/0.74/(LOG(Z/Z0)+7.*MIN(XX,0.714))**2.
1806 CD=0.16/(LOG(Z/Z0)+7.*MIN(XX,0.714))**2.
1807 ELSE
1808 CMB=7.4*A2*9.4*SQRT(Z/Z0)
1809 CHB=5.3*A2*9.4*SQRT(Z/Z0)
1810 CH=A2/0.74*(1.-9.4*RIB/(1.+CHB*SQRT(-RIB)))
1811 CD=A2*(1.-9.4*RIB/(1.+CHB*SQRT(-RIB)))
1812 END IF
1814 ALPHA=RHO*CP*CH*UA
1816 RETURN
1817 END SUBROUTINE louis79
1818 !===============================================================================
1819 !
1820 ! louis82
1821 !
1822 !===============================================================================
1823 SUBROUTINE louis82(ALPHA,CD,RIB,Z,Z0,UA,RHO)
1825 IMPLICIT NONE
1827 REAL, PARAMETER :: CP=0.24
1828 REAL, INTENT(IN) :: Z, Z0, UA, RHO
1829 REAL, INTENT(OUT) :: ALPHA, CD
1830 REAL, INTENT(INOUT) :: RIB
1831 REAL :: A2, FM, FH, CH, CHH
1833 A2=(0.4/ALOG(Z/Z0))**2.
1835 IF(RIB <= -15.) RIB=-15.
1837 IF(RIB >= 0.0) THEN
1838 FM=1./((1.+(2.*5.*RIB)/SQRT(1.+5.*RIB)))
1839 FH=1./(1.+(3.*5.*RIB)*SQRT(1.+5.*RIB))
1840 CH=A2*FH
1841 CD=A2*FM
1842 ELSE
1843 CHH=5.*3.*5.*A2*SQRT(Z/Z0)
1844 FM=1.-(2.*5.*RIB)/(1.+3.*5.*5.*A2*SQRT(Z/Z0+1.)*(-RIB))
1845 FH=1.-(3.*5.*RIB)/(1.+CHH*SQRT(-RIB))
1846 CH=A2*FH
1847 CD=A2*FM
1848 END IF
1850 ALPHA=RHO*CP*CH*UA
1852 RETURN
1853 END SUBROUTINE louis82
1854 !===============================================================================
1855 !
1856 ! multi_layer
1857 !
1858 !===============================================================================
1859 SUBROUTINE multi_layer(KM,BOUND,G0,CAP,AKS,TSL,DZ,DELT,TSLEND)
1861 IMPLICIT NONE
1863 REAL, INTENT(IN) :: G0
1865 REAL, INTENT(IN) :: CAP
1867 REAL, INTENT(IN) :: AKS
1869 REAL, INTENT(IN) :: DELT ! Time step [ s ]
1871 REAL, INTENT(IN) :: TSLEND
1873 INTEGER, INTENT(IN) :: KM
1875 INTEGER, INTENT(IN) :: BOUND
1877 REAL, DIMENSION(KM), INTENT(IN) :: DZ
1879 REAL, DIMENSION(KM), INTENT(INOUT) :: TSL
1881 REAL, DIMENSION(KM) :: A, B, C, D, X, P, Q
1883 REAL :: DZEND
1885 INTEGER :: K
1887 DZEND=DZ(KM)
1889 A(1) = 0.0
1891 B(1) = CAP*DZ(1)/DELT &
1892 +2.*AKS/(DZ(1)+DZ(2))
1893 C(1) = -2.*AKS/(DZ(1)+DZ(2))
1894 D(1) = CAP*DZ(1)/DELT*TSL(1) + G0
1896 DO K=2,KM-1
1897 A(K) = -2.*AKS/(DZ(K-1)+DZ(K))
1898 B(K) = CAP*DZ(K)/DELT + 2.*AKS/(DZ(K-1)+DZ(K)) + 2.*AKS/(DZ(K)+DZ(K+1))
1899 C(K) = -2.*AKS/(DZ(K)+DZ(K+1))
1900 D(K) = CAP*DZ(K)/DELT*TSL(K)
1901 END DO
1903 IF(BOUND == 1) THEN ! Flux=0
1904 A(KM) = -2.*AKS/(DZ(KM-1)+DZ(KM))
1905 B(KM) = CAP*DZ(KM)/DELT + 2.*AKS/(DZ(KM-1)+DZ(KM))
1906 C(KM) = 0.0
1907 D(KM) = CAP*DZ(KM)/DELT*TSL(KM)
1908 ELSE ! T=constant
1909 A(KM) = -2.*AKS/(DZ(KM-1)+DZ(KM))
1910 B(KM) = CAP*DZ(KM)/DELT + 2.*AKS/(DZ(KM-1)+DZ(KM)) + 2.*AKS/(DZ(KM)+DZEND)
1911 C(KM) = 0.0
1912 D(KM) = CAP*DZ(KM)/DELT*TSL(KM) + 2.*AKS*TSLEND/(DZ(KM)+DZEND)
1913 END IF
1915 P(1) = -C(1)/B(1)
1916 Q(1) = D(1)/B(1)
1918 DO K=2,KM
1919 P(K) = -C(K)/(A(K)*P(K-1)+B(K))
1920 Q(K) = (-A(K)*Q(K-1)+D(K))/(A(K)*P(K-1)+B(K))
1921 END DO
1923 X(KM) = Q(KM)
1925 DO K=KM-1,1,-1
1926 X(K) = P(K)*X(K+1)+Q(K)
1927 END DO
1929 DO K=1,KM
1930 TSL(K) = X(K)
1931 END DO
1933 RETURN
1934 END SUBROUTINE multi_layer
1935 !===============================================================================
1936 !
1937 ! subroutine read_param
1938 !
1939 !===============================================================================
1940 SUBROUTINE read_param(UTYPE, & ! in
1941 ZR,SIGMA_ZED,Z0C,Z0HC,ZDC,SVF,R,RW,HGT,AH, & ! out
1942 CAPR,CAPB,CAPG,AKSR,AKSB,AKSG,ALBR,ALBB,ALBG, & ! out
1943 EPSR,EPSB,EPSG,Z0R,Z0B,Z0G,Z0HB,Z0HG, & ! out
1944 BETR,BETB,BETG,TRLEND,TBLEND,TGLEND, & ! out
1945 !for BEP
1946 NUMDIR, STREET_DIRECTION, STREET_WIDTH, & ! out
1947 BUILDING_WIDTH, NUMHGT, HEIGHT_BIN, & ! out
1948 HPERCENT_BIN, & ! out
1949 !end BEP
1950 BOUNDR,BOUNDB,BOUNDG,CH_SCHEME,TS_SCHEME, & ! out
1951 AKANDA_URBAN,ALH) ! out
1953 INTEGER, INTENT(IN) :: UTYPE
1955 REAL, INTENT(OUT) :: ZR,Z0C,Z0HC,ZDC,SVF,R,RW,HGT,AH,ALH, &
1956 CAPR,CAPB,CAPG,AKSR,AKSB,AKSG,ALBR,ALBB,ALBG, &
1957 SIGMA_ZED, &
1958 EPSR,EPSB,EPSG,Z0R,Z0B,Z0G,Z0HB,Z0HG, &
1959 BETR,BETB,BETG,TRLEND,TBLEND,TGLEND
1960 REAL, INTENT(OUT) :: AKANDA_URBAN
1961 !for BEP
1962 INTEGER, INTENT(OUT) :: NUMDIR
1963 REAL, DIMENSION(MAXDIRS), INTENT(OUT) :: STREET_DIRECTION
1964 REAL, DIMENSION(MAXDIRS), INTENT(OUT) :: STREET_WIDTH
1965 REAL, DIMENSION(MAXDIRS), INTENT(OUT) :: BUILDING_WIDTH
1966 INTEGER, INTENT(OUT) :: NUMHGT
1967 REAL, DIMENSION(MAXHGTS), INTENT(OUT) :: HEIGHT_BIN
1968 REAL, DIMENSION(MAXHGTS), INTENT(OUT) :: HPERCENT_BIN
1970 !end BEP
1972 INTEGER, INTENT(OUT) :: BOUNDR,BOUNDB,BOUNDG,CH_SCHEME,TS_SCHEME
1974 ZR = ZR_TBL(UTYPE)
1975 SIGMA_ZED = SIGMA_ZED_TBL(UTYPE)
1976 Z0C= Z0C_TBL(UTYPE)
1977 Z0HC= Z0HC_TBL(UTYPE)
1978 ZDC= ZDC_TBL(UTYPE)
1979 SVF= SVF_TBL(UTYPE)
1980 R= R_TBL(UTYPE)
1981 RW= RW_TBL(UTYPE)
1982 HGT= HGT_TBL(UTYPE)
1983 AH= AH_TBL(UTYPE)
1984 ALH= ALH_TBL(UTYPE)
1985 BETR= BETR_TBL(UTYPE)
1986 BETB= BETB_TBL(UTYPE)
1987 BETG= BETG_TBL(UTYPE)
1989 !m FRC_URB= FRC_URB_TBL(UTYPE)
1991 CAPR= CAPR_TBL(UTYPE)
1992 CAPB= CAPB_TBL(UTYPE)
1993 CAPG= CAPG_TBL(UTYPE)
1994 AKSR= AKSR_TBL(UTYPE)
1995 AKSB= AKSB_TBL(UTYPE)
1996 AKSG= AKSG_TBL(UTYPE)
1997 ALBR= ALBR_TBL(UTYPE)
1998 ALBB= ALBB_TBL(UTYPE)
1999 ALBG= ALBG_TBL(UTYPE)
2000 EPSR= EPSR_TBL(UTYPE)
2001 EPSB= EPSB_TBL(UTYPE)
2002 EPSG= EPSG_TBL(UTYPE)
2003 Z0R= Z0R_TBL(UTYPE)
2004 Z0B= Z0B_TBL(UTYPE)
2005 Z0G= Z0G_TBL(UTYPE)
2006 Z0HB= Z0HB_TBL(UTYPE)
2007 Z0HG= Z0HG_TBL(UTYPE)
2008 TRLEND= TRLEND_TBL(UTYPE)
2009 TBLEND= TBLEND_TBL(UTYPE)
2010 TGLEND= TGLEND_TBL(UTYPE)
2011 BOUNDR= BOUNDR_DATA
2012 BOUNDB= BOUNDB_DATA
2013 BOUNDG= BOUNDG_DATA
2014 CH_SCHEME = CH_SCHEME_DATA
2015 TS_SCHEME = TS_SCHEME_DATA
2016 AKANDA_URBAN = AKANDA_URBAN_TBL(UTYPE)
2018 !for BEP
2020 STREET_DIRECTION = -1.E36
2021 STREET_WIDTH = -1.E36
2022 BUILDING_WIDTH = -1.E36
2023 HEIGHT_BIN = -1.E36
2024 HPERCENT_BIN = -1.E36
2026 NUMDIR = NUMDIR_TBL ( UTYPE )
2027 STREET_DIRECTION(1:NUMDIR) = STREET_DIRECTION_TBL( 1:NUMDIR, UTYPE )
2028 STREET_WIDTH (1:NUMDIR) = STREET_WIDTH_TBL ( 1:NUMDIR, UTYPE )
2029 BUILDING_WIDTH (1:NUMDIR) = BUILDING_WIDTH_TBL ( 1:NUMDIR, UTYPE )
2030 NUMHGT = NUMHGT_TBL ( UTYPE )
2031 HEIGHT_BIN (1:NUMHGT) = HEIGHT_BIN_TBL ( 1:NUMHGT , UTYPE )
2032 HPERCENT_BIN (1:NUMHGT) = HPERCENT_BIN_TBL ( 1:NUMHGT , UTYPE )
2034 !end BEP
2035 END SUBROUTINE read_param
2036 !===============================================================================
2037 !
2038 ! subroutine urban_param_init: Read parameters from URBPARM.TBL
2039 !
2040 !===============================================================================
2041 SUBROUTINE urban_param_init(DZR,DZB,DZG,num_soil_layers, &
2042 sf_urban_physics,use_wudapt_lcz, slucm_distributed_drag)
2043 ! num_roof_layers,num_wall_layers,num_road_layers)
2045 IMPLICIT NONE
2047 INTEGER, INTENT(IN) :: num_soil_layers
2049 ! REAL, DIMENSION(1:num_roof_layers), INTENT(INOUT) :: DZR
2050 ! REAL, DIMENSION(1:num_wall_layers), INTENT(INOUT) :: DZB
2051 ! REAL, DIMENSION(1:num_road_layers), INTENT(INOUT) :: DZG
2052 REAL, DIMENSION(1:num_soil_layers), INTENT(INOUT) :: DZR
2053 REAL, DIMENSION(1:num_soil_layers), INTENT(INOUT) :: DZB
2054 REAL, DIMENSION(1:num_soil_layers), INTENT(INOUT) :: DZG
2055 INTEGER, INTENT(IN) :: SF_URBAN_PHYSICS
2056 INTEGER, INTENT(IN) :: USE_WUDAPT_LCZ !AndreaLCZ
2057 LOGICAL, INTENT(IN) :: slucm_distributed_drag
2059 INTEGER :: LC, K
2060 INTEGER :: IOSTATUS, ALLOCATE_STATUS
2061 INTEGER :: num_roof_layers
2062 INTEGER :: num_wall_layers
2063 INTEGER :: num_road_layers
2064 INTEGER :: dummy
2065 REAL :: DHGT, HGT, VFWS, VFGS
2067 REAL, allocatable, dimension(:) :: ROOF_WIDTH
2068 REAL, allocatable, dimension(:) :: ROAD_WIDTH
2070 character(len=512) :: string
2071 character(len=128) :: name
2072 integer :: indx
2074 real, parameter :: VonK = 0.4
2075 real :: lambda_p
2076 real :: lambda_f
2077 real :: Cd
2078 real :: alpha_macd
2079 real :: beta_macd
2080 real :: lambda_fr
2083 !for BEP
2084 real :: dummy_hgt
2085 real :: dummy_pct
2086 real :: pctsum
2087 !end BEP
2088 num_roof_layers = num_soil_layers
2089 num_wall_layers = num_soil_layers
2090 num_road_layers = num_soil_layers
2093 ICATE=0
2095 distributed_aerodynamics_option = slucm_distributed_drag
2098 if(USE_WUDAPT_LCZ.eq.0)then !AndreaLCZ
2099 OPEN (UNIT=11, &
2100 FILE='URBPARM.TBL', &
2101 ACCESS='SEQUENTIAL', &
2102 STATUS='OLD', &
2103 ACTION='READ', &
2104 POSITION='REWIND', &
2105 IOSTAT=IOSTATUS)
2107 IF (IOSTATUS > 0) THEN
2108 FATAL_ERROR('Error opening URBPARM.TBL. Make sure URBPARM.TBL (found in run/) is linked to the running directory.')
2109 ENDIF
2111 else
2112 OPEN (UNIT=11, &
2113 FILE='URBPARM_LCZ.TBL', &
2114 ACCESS='SEQUENTIAL', &
2115 STATUS='OLD', &
2116 ACTION='READ', &
2117 POSITION='REWIND', &
2118 IOSTAT=IOSTATUS)
2120 IF (IOSTATUS > 0) THEN
2121 FATAL_ERROR('Error opening URBPARM_LCZ.TBL. Make sure URBPARM_LCZ.TBL (found in run/) is linked to the running directory.')
2122 ENDIF
2123 endif
2126 READLOOP : do
2127 read(11,'(A512)', iostat=iostatus) string
2128 if (iostatus /= 0) exit READLOOP
2129 if (string(1:1) == "#") cycle READLOOP
2130 if (trim(string) == "") cycle READLOOP
2131 indx = index(string,":")
2132 if (indx <= 0) cycle READLOOP
2133 name = trim(adjustl(string(1:indx-1)))
2135 ! Here are the variables we expect to be defined in the URBPARM.TBL:
2136 if (name == "Number of urban categories") then
2137 read(string(indx+1:),*) icate
2138 IF (.not. ALLOCATED(ZR_TBL)) then
2139 ALLOCATE( ZR_TBL(ICATE), stat=allocate_status )
2140 if(allocate_status /= 0) FATAL_ERROR('Error allocating ZR_TBL in urban_param_init')
2141 ALLOCATE( SIGMA_ZED_TBL(ICATE), stat=allocate_status )
2142 if(allocate_status /= 0) FATAL_ERROR('Error allocating SIGMA_ZED_TBL in urban_param_init')
2143 ALLOCATE( Z0C_TBL(ICATE), stat=allocate_status )
2144 if(allocate_status /= 0) FATAL_ERROR('Error allocating Z0C_TBL in urban_param_init')
2145 ALLOCATE( Z0HC_TBL(ICATE), stat=allocate_status )
2146 if(allocate_status /= 0) FATAL_ERROR('Error allocating Z0HC_TBL in urban_param_init')
2147 ALLOCATE( ZDC_TBL(ICATE), stat=allocate_status )
2148 if(allocate_status /= 0) FATAL_ERROR('Error allocating ZDC_TBL in urban_param_init')
2149 ALLOCATE( SVF_TBL(ICATE), stat=allocate_status )
2150 if(allocate_status /= 0) FATAL_ERROR('Error allocating SVF_TBL in urban_param_init')
2151 ALLOCATE( R_TBL(ICATE), stat=allocate_status )
2152 if(allocate_status /= 0) FATAL_ERROR('Error allocating R_TBL in urban_param_init')
2153 ALLOCATE( RW_TBL(ICATE), stat=allocate_status )
2154 if(allocate_status /= 0) FATAL_ERROR('Error allocating RW_TBL in urban_param_init')
2155 ALLOCATE( HGT_TBL(ICATE), stat=allocate_status )
2156 if(allocate_status /= 0) FATAL_ERROR('Error allocating HGT_TBL in urban_param_init')
2157 ALLOCATE( AH_TBL(ICATE), stat=allocate_status )
2158 if(allocate_status /= 0) FATAL_ERROR('Error allocating AH_TBL in urban_param_init')
2159 ALLOCATE( ALH_TBL(ICATE), stat=allocate_status )
2160 if(allocate_status /= 0) FATAL_ERROR('Error allocating ALH_TBL in urban_param_init')
2161 ALLOCATE( BETR_TBL(ICATE), stat=allocate_status )
2162 if(allocate_status /= 0) FATAL_ERROR('Error allocating BETR_TBL in urban_param_init')
2163 ALLOCATE( BETB_TBL(ICATE), stat=allocate_status )
2164 if(allocate_status /= 0) FATAL_ERROR('Error allocating BETB_TBL in urban_param_init')
2165 ALLOCATE( BETG_TBL(ICATE), stat=allocate_status )
2166 if(allocate_status /= 0) FATAL_ERROR('Error allocating BETG_TBL in urban_param_init')
2167 ALLOCATE( CAPR_TBL(ICATE), stat=allocate_status )
2168 if(allocate_status /= 0) FATAL_ERROR('Error allocating CAPR_TBL in urban_param_init')
2169 ALLOCATE( CAPB_TBL(ICATE), stat=allocate_status )
2170 if(allocate_status /= 0) FATAL_ERROR('Error allocating CAPB_TBL in urban_param_init')
2171 ALLOCATE( CAPG_TBL(ICATE), stat=allocate_status )
2172 if(allocate_status /= 0) FATAL_ERROR('Error allocating CAPG_TBL in urban_param_init')
2173 ALLOCATE( AKSR_TBL(ICATE), stat=allocate_status )
2174 if(allocate_status /= 0) FATAL_ERROR('Error allocating AKSR_TBL in urban_param_init')
2175 ALLOCATE( AKSB_TBL(ICATE), stat=allocate_status )
2176 if(allocate_status /= 0) FATAL_ERROR('Error allocating AKSB_TBL in urban_param_init')
2177 ALLOCATE( AKSG_TBL(ICATE), stat=allocate_status )
2178 if(allocate_status /= 0) FATAL_ERROR('Error allocating AKSG_TBL in urban_param_init')
2179 ALLOCATE( ALBR_TBL(ICATE), stat=allocate_status )
2180 if(allocate_status /= 0) FATAL_ERROR('Error allocating ALBR_TBL in urban_param_init')
2181 ALLOCATE( ALBB_TBL(ICATE), stat=allocate_status )
2182 if(allocate_status /= 0) FATAL_ERROR('Error allocating ALBB_TBL in urban_param_init')
2183 ALLOCATE( ALBG_TBL(ICATE), stat=allocate_status )
2184 if(allocate_status /= 0) FATAL_ERROR('Error allocating ALBG_TBL in urban_param_init')
2185 ALLOCATE( EPSR_TBL(ICATE), stat=allocate_status )
2186 if(allocate_status /= 0) FATAL_ERROR('Error allocating EPSR_TBL in urban_param_init')
2187 ALLOCATE( EPSB_TBL(ICATE), stat=allocate_status )
2188 if(allocate_status /= 0) FATAL_ERROR('Error allocating EPSB_TBL in urban_param_init')
2189 ALLOCATE( EPSG_TBL(ICATE), stat=allocate_status )
2190 if(allocate_status /= 0) FATAL_ERROR('Error allocating EPSG_TBL in urban_param_init')
2191 ALLOCATE( Z0R_TBL(ICATE), stat=allocate_status )
2192 if(allocate_status /= 0) FATAL_ERROR('Error allocating Z0R_TBL in urban_param_init')
2193 ALLOCATE( Z0B_TBL(ICATE), stat=allocate_status )
2194 if(allocate_status /= 0) FATAL_ERROR('Error allocating Z0B_TBL in urban_param_init')
2195 ALLOCATE( Z0G_TBL(ICATE), stat=allocate_status )
2196 if(allocate_status /= 0) FATAL_ERROR('Error allocating Z0G_TBL in urban_param_init')
2197 ALLOCATE( AKANDA_URBAN_TBL(ICATE), stat=allocate_status )
2198 if(allocate_status /= 0) FATAL_ERROR('Error allocating AKANDA_URBAN_TBL in urban_param_init')
2199 ALLOCATE( Z0HB_TBL(ICATE), stat=allocate_status )
2200 if(allocate_status /= 0) FATAL_ERROR('Error allocating Z0HB_TBL in urban_param_init')
2201 ALLOCATE( Z0HG_TBL(ICATE), stat=allocate_status )
2202 if(allocate_status /= 0) FATAL_ERROR('Error allocating Z0HG_TBL in urban_param_init')
2203 ALLOCATE( TRLEND_TBL(ICATE), stat=allocate_status )
2204 if(allocate_status /= 0) FATAL_ERROR('Error allocating TRLEND_TBL in urban_param_init')
2205 ALLOCATE( TBLEND_TBL(ICATE), stat=allocate_status )
2206 if(allocate_status /= 0) FATAL_ERROR('Error allocating TBLEND_TBL in urban_param_init')
2207 ALLOCATE( TGLEND_TBL(ICATE), stat=allocate_status )
2208 if(allocate_status /= 0) FATAL_ERROR('Error allocating TGLEND_TBL in urban_param_init')
2209 ALLOCATE( FRC_URB_TBL(ICATE), stat=allocate_status )
2210 if(allocate_status /= 0) FATAL_ERROR('Error allocating FRC_URB_TBL in urban_param_init')
2211 ! ALLOCATE( ROOF_WIDTH(ICATE), stat=allocate_status )
2212 ! if(allocate_status /= 0) FATAL_ERROR('Error allocating ROOF_WIDTH in urban_param_init')
2213 ! ALLOCATE( ROAD_WIDTH(ICATE), stat=allocate_status )
2214 ! if(allocate_status /= 0) FATAL_ERROR('Error allocating ROAD_WIDTH in urban_param_init')
2215 !for BEP
2216 ALLOCATE( NUMDIR_TBL(ICATE), stat=allocate_status )
2217 if(allocate_status /= 0) FATAL_ERROR('Error allocating NUMDIR_TBL in urban_param_init')
2218 ALLOCATE( STREET_DIRECTION_TBL(MAXDIRS , ICATE), stat=allocate_status )
2219 if(allocate_status /= 0) FATAL_ERROR('Error allocating STREET_DIRECTION_TBL in urban_param_init')
2220 ALLOCATE( STREET_WIDTH_TBL(MAXDIRS , ICATE), stat=allocate_status )
2221 if(allocate_status /= 0) FATAL_ERROR('Error allocating STREET_WIDTH_TBL in urban_param_init')
2222 ALLOCATE( BUILDING_WIDTH_TBL(MAXDIRS , ICATE), stat=allocate_status )
2223 if(allocate_status /= 0) FATAL_ERROR('Error allocating BUILDING_WIDTH_TBL in urban_param_init')
2224 ALLOCATE( NUMHGT_TBL(ICATE), stat=allocate_status )
2225 if(allocate_status /= 0) FATAL_ERROR('Error allocating NUMHGT_TBL in urban_param_init')
2226 ALLOCATE( HEIGHT_BIN_TBL(MAXHGTS , ICATE), stat=allocate_status )
2227 if(allocate_status /= 0) FATAL_ERROR('Error allocating HEIGHT_BIN_TBL in urban_param_init')
2228 ALLOCATE( HPERCENT_BIN_TBL(MAXHGTS , ICATE), stat=allocate_status )
2229 if(allocate_status /= 0) FATAL_ERROR('Error allocating HPERCENT_BIN_TBL in urban_param_init')
2230 ALLOCATE( COP_TBL(ICATE), stat=allocate_status )
2231 if(allocate_status /= 0) FATAL_ERROR('Error allocating COP_TBL in urban_param_init')
2232 ALLOCATE( BLDAC_FRC_TBL(ICATE), stat=allocate_status )
2233 if(allocate_status /= 0) FATAL_ERROR('Error allocating BLDAC_FRC_TBL in urban_param_init')
2234 ALLOCATE( COOLED_FRC_TBL(ICATE), stat=allocate_status )
2235 if(allocate_status /= 0) FATAL_ERROR('Error allocating COOLED_FRC_TBL in urban_param_init')
2236 ALLOCATE( PWIN_TBL(ICATE), stat=allocate_status )
2237 if(allocate_status /= 0) FATAL_ERROR('Error allocating PWIN_TBL in urban_param_init')
2238 ALLOCATE( BETA_TBL(ICATE), stat=allocate_status )
2239 if(allocate_status /= 0) FATAL_ERROR('Error allocating BETA_TBL in urban_param_init')
2240 ALLOCATE( SW_COND_TBL(ICATE), stat=allocate_status )
2241 if(allocate_status /= 0) FATAL_ERROR('Error allocating SW_COND_TBL in urban_param_init')
2242 ALLOCATE( TIME_ON_TBL(ICATE), stat=allocate_status )
2243 if(allocate_status /= 0) FATAL_ERROR('Error allocating TIME_ON_TBL in urban_param_init')
2244 ALLOCATE( TIME_OFF_TBL(ICATE), stat=allocate_status )
2245 if(allocate_status /= 0) FATAL_ERROR('Error allocating TIME_OFF_TBL in urban_param_init')
2246 ALLOCATE( TARGTEMP_TBL(ICATE), stat=allocate_status )
2247 if(allocate_status /= 0) FATAL_ERROR('Error allocating TARGTEMP_TBL in urban_param_init')
2248 ALLOCATE( GAPTEMP_TBL(ICATE), stat=allocate_status )
2249 if(allocate_status /= 0) FATAL_ERROR('Error allocating GAPTEMP_TBL in urban_param_init')
2250 ALLOCATE( TARGHUM_TBL(ICATE), stat=allocate_status )
2251 if(allocate_status /= 0) FATAL_ERROR('Error allocating TARGHUM_TBL in urban_param_init')
2252 ALLOCATE( GAPHUM_TBL(ICATE), stat=allocate_status )
2253 if(allocate_status /= 0) FATAL_ERROR('Error allocating GAPHUM_TBL in urban_param_init')
2254 ALLOCATE( PERFLO_TBL(ICATE), stat=allocate_status )
2255 if(allocate_status /= 0) FATAL_ERROR('Error allocating PERFLO_TBL in urban_param_init')
2256 ALLOCATE( HSESF_TBL(ICATE), stat=allocate_status )
2257 if(allocate_status /= 0) FATAL_ERROR('Error allocating HSESF_TBL in urban_param_init')
2258 ALLOCATE( PV_FRAC_ROOF_TBL(ICATE), stat=allocate_status )
2259 if(allocate_status /= 0) FATAL_ERROR('Error allocating PV_FRAC_ROOF_TBL in urban_param_init')
2260 ALLOCATE( GR_FRAC_ROOF_TBL(ICATE), stat=allocate_status )
2261 if(allocate_status /= 0) FATAL_ERROR('Error allocating GR_FRAC_ROOF_TBL in urban_param_init')
2263 endif
2264 numdir_tbl = 0
2265 street_direction_tbl = -1.E36
2266 street_width_tbl = 0
2267 building_width_tbl = 0
2268 numhgt_tbl = 0
2269 height_bin_tbl = -1.E36
2270 hpercent_bin_tbl = -1.E36
2271 !end BEP
2273 else if (name == "ZR") then
2274 read(string(indx+1:),*) zr_tbl(1:icate)
2275 else if (name == "SIGMA_ZED") then
2276 read(string(indx+1:),*) sigma_zed_tbl(1:icate)
2277 else if (name == "ROOF_WIDTH") then
2278 ALLOCATE( ROOF_WIDTH(ICATE), stat=allocate_status )
2279 if(allocate_status /= 0) FATAL_ERROR('Error allocating ROOF_WIDTH in urban_param_init')
2281 read(string(indx+1:),*) roof_width(1:icate)
2282 else if (name == "ROAD_WIDTH") then
2283 ALLOCATE( ROAD_WIDTH(ICATE), stat=allocate_status )
2284 if(allocate_status /= 0) FATAL_ERROR('Error allocating ROAD_WIDTH in urban_param_init')
2285 read(string(indx+1:),*) road_width(1:icate)
2286 else if (name == "AH") then
2287 read(string(indx+1:),*) ah_tbl(1:icate)
2288 else if (name == "ALH") then
2289 read(string(indx+1:),*) alh_tbl(1:icate)
2290 else if (name == "FRC_URB") then
2291 read(string(indx+1:),*) frc_urb_tbl(1:icate)
2292 else if (name == "CAPR") then
2293 read(string(indx+1:),*) capr_tbl(1:icate)
2294 ! Convert CAPR_TBL from J m{-3} K{-1} to cal cm{-3} deg{-1}
2295 capr_tbl = capr_tbl * ( 1.0 / 4.1868 ) * 1.E-6
2296 else if (name == "CAPB") then
2297 read(string(indx+1:),*) capb_tbl(1:icate)
2298 ! Convert CABR_TBL from J m{-3} K{-1} to cal cm{-3} deg{-1}
2299 capb_tbl = capb_tbl * ( 1.0 / 4.1868 ) * 1.E-6
2300 else if (name == "CAPG") then
2301 read(string(indx+1:),*) capg_tbl(1:icate)
2302 ! Convert CABG_TBL from J m{-3} K{-1} to cal cm{-3} deg{-1}
2303 capg_tbl = capg_tbl * ( 1.0 / 4.1868 ) * 1.E-6
2304 else if (name == "AKSR") then
2305 read(string(indx+1:),*) aksr_tbl(1:icate)
2306 ! Convert AKSR_TBL from J m{-1} s{-1} K{-1} to cal cm{-1} s{-1} deg{-1}
2307 AKSR_TBL = AKSR_TBL * ( 1.0 / 4.1868 ) * 1.E-2
2308 else if (name == "AKSB") then
2309 read(string(indx+1:),*) aksb_tbl(1:icate)
2310 ! Convert AKSB_TBL from J m{-1} s{-1} K{-1} to cal cm{-1} s{-1} deg{-1}
2311 AKSB_TBL = AKSB_TBL * ( 1.0 / 4.1868 ) * 1.E-2
2312 else if (name == "AKSG") then
2313 read(string(indx+1:),*) aksg_tbl(1:icate)
2314 ! Convert AKSG_TBL from J m{-1} s{-1} K{-1} to cal cm{-1} s{-1} deg{-1}
2315 AKSG_TBL = AKSG_TBL * ( 1.0 / 4.1868 ) * 1.E-2
2316 else if (name == "ALBR") then
2317 read(string(indx+1:),*) albr_tbl(1:icate)
2318 else if (name == "ALBB") then
2319 read(string(indx+1:),*) albb_tbl(1:icate)
2320 else if (name == "ALBG") then
2321 read(string(indx+1:),*) albg_tbl(1:icate)
2322 else if (name == "EPSR") then
2323 read(string(indx+1:),*) epsr_tbl(1:icate)
2324 else if (name == "EPSB") then
2325 read(string(indx+1:),*) epsb_tbl(1:icate)
2326 else if (name == "EPSG") then
2327 read(string(indx+1:),*) epsg_tbl(1:icate)
2328 else if (name == "AKANDA_URBAN") then
2329 read(string(indx+1:),*) akanda_urban_tbl(1:icate)
2330 else if (name == "Z0B") then
2331 read(string(indx+1:),*) z0b_tbl(1:icate)
2332 else if (name == "Z0G") then
2333 read(string(indx+1:),*) z0g_tbl(1:icate)
2334 else if (name == "DDZR") then
2335 read(string(indx+1:),*) dzr(1:num_roof_layers)
2336 ! Convert thicknesses from m to cm
2337 dzr = dzr * 100.0
2338 else if (name == "DDZB") then
2339 read(string(indx+1:),*) dzb(1:num_wall_layers)
2340 ! Convert thicknesses from m to cm
2341 dzb = dzb * 100.0
2342 else if (name == "DDZG") then
2343 read(string(indx+1:),*) dzg(1:num_road_layers)
2344 ! Convert thicknesses from m to cm
2345 dzg = dzg * 100.0
2346 else if (name == "BOUNDR") then
2347 read(string(indx+1:),*) boundr_data
2348 else if (name == "BOUNDB") then
2349 read(string(indx+1:),*) boundb_data
2350 else if (name == "BOUNDG") then
2351 read(string(indx+1:),*) boundg_data
2352 else if (name == "TRLEND") then
2353 read(string(indx+1:),*) trlend_tbl(1:icate)
2354 else if (name == "TBLEND") then
2355 read(string(indx+1:),*) tblend_tbl(1:icate)
2356 else if (name == "TGLEND") then
2357 read(string(indx+1:),*) tglend_tbl(1:icate)
2358 else if (name == "CH_SCHEME") then
2359 read(string(indx+1:),*) ch_scheme_data
2360 else if (name == "TS_SCHEME") then
2361 read(string(indx+1:),*) ts_scheme_data
2362 else if (name == "AHOPTION") then
2363 read(string(indx+1:),*) ahoption
2364 else if (name == "AHDIUPRF") then
2365 read(string(indx+1:),*) ahdiuprf(1:24)
2366 else if (name == "ALHOPTION") then
2367 read(string(indx+1:),*) alhoption
2368 else if (name == "ALHSEASON") then
2369 read(string(indx+1:),*) alhseason(1:4)
2370 else if (name == "ALHDIUPRF") then
2371 read(string(indx+1:),*) alhdiuprf(1:48)
2372 else if (name == "PORIMP") then
2373 read(string(indx+1:),*) porimp(1:3)
2374 else if (name == "DENGIMP") then
2375 read(string(indx+1:),*) dengimp(1:3)
2376 else if (name == "IMP_SCHEME") then
2377 read(string(indx+1:),*) imp_scheme
2378 else if (name == "IRI_SCHEME") then
2379 read(string(indx+1:),*) iri_scheme
2380 else if (name == "OASIS") then
2381 read(string(indx+1:),*) oasis
2382 else if (name == "GROPTION") then
2383 read(string(indx+1:),*) groption
2384 else if (name == "FGR") then
2385 read(string(indx+1:),*) fgr
2386 else if (name == "DZGR") then
2387 read(string(indx+1:),*) dzgr(1:4)
2388 !for BEP
2389 else if (name == "STREET PARAMETERS") then
2391 STREETLOOP : do
2392 read(11,'(A512)', iostat=iostatus) string
2393 if (string(1:1) == "#") cycle STREETLOOP
2394 if (trim(string) == "") cycle STREETLOOP
2395 if (string == "END STREET PARAMETERS") exit STREETLOOP
2396 read(string, *) k ! , dirst, ws, bs
2397 numdir_tbl(k) = numdir_tbl(k) + 1
2398 read(string, *) k, street_direction_tbl(numdir_tbl(k),k), &
2399 street_width_tbl(numdir_tbl(k),k), &
2400 building_width_tbl(numdir_tbl(k),k)
2401 enddo STREETLOOP
2403 else if (name == "BUILDING HEIGHTS") then
2405 read(string(indx+1:),*) k
2406 HEIGHTLOOP : do
2407 read(11,'(A512)', iostat=iostatus) string
2408 if (string(1:1) == "#") cycle HEIGHTLOOP
2409 if (trim(string) == "") cycle HEIGHTLOOP
2410 if (string == "END BUILDING HEIGHTS") exit HEIGHTLOOP
2411 read(string,*) dummy_hgt, dummy_pct
2412 numhgt_tbl(k) = numhgt_tbl(k) + 1
2413 height_bin_tbl(numhgt_tbl(k), k) = dummy_hgt
2414 hpercent_bin_tbl(numhgt_tbl(k),k) = dummy_pct
2416 enddo HEIGHTLOOP
2417 pctsum = sum ( hpercent_bin_tbl(:,k) , mask=(hpercent_bin_tbl(:,k)>-1.E25 ) )
2418 if ( pctsum /= 100.) then
2419 write (*,'(//,"Building height percentages for category ", I2, " must sum to 100.0")') k
2420 write (*,'("Currently, they sum to ", F6.2,/)') pctsum
2421 FATAL_ERROR('pctsum is not equal to 100.')
2422 endif
2423 else if ( name == "Z0R") then
2424 read(string(indx+1:),*) Z0R_tbl(1:icate)
2425 else if ( name == "COP") then
2426 read(string(indx+1:),*) cop_tbl(1:icate)
2427 else if ( name == "BLDAC_FRC") then
2428 read(string(indx+1:),*) bldac_frc_tbl(1:icate)
2429 else if ( name == "COOLED_FRC") then
2430 read(string(indx+1:),*) cooled_frc_tbl(1:icate)
2431 else if ( name == "PWIN") then
2432 read(string(indx+1:),*) pwin_tbl(1:icate)
2433 else if ( name == "BETA") then
2434 read(string(indx+1:),*) beta_tbl(1:icate)
2435 else if ( name == "SW_COND") then
2436 read(string(indx+1:),*) sw_cond_tbl(1:icate)
2437 else if ( name == "TIME_ON") then
2438 read(string(indx+1:),*) time_on_tbl(1:icate)
2439 else if ( name == "TIME_OFF") then
2440 read(string(indx+1:),*) time_off_tbl(1:icate)
2441 else if ( name == "TARGTEMP") then
2442 read(string(indx+1:),*) targtemp_tbl(1:icate)
2443 else if ( name == "GAPTEMP") then
2444 read(string(indx+1:),*) gaptemp_tbl(1:icate)
2445 else if ( name == "TARGHUM") then
2446 read(string(indx+1:),*) targhum_tbl(1:icate)
2447 else if ( name == "GAPHUM") then
2448 read(string(indx+1:),*) gaphum_tbl(1:icate)
2449 else if ( name == "PERFLO") then
2450 read(string(indx+1:),*) perflo_tbl(1:icate)
2451 else if (name == "HSEQUIP") then
2452 read(string(indx+1:),*) hsequip_tbl(1:24)
2453 else if (name == "HSEQUIP_SCALE_FACTOR") then
2454 read(string(indx+1:),*) hsesf_tbl(1:icate)
2455 else if (name == "IRHO") then
2456 read(string(indx+1:),*) IRHO_TBL(1:24)
2457 else if ( name == "PV_FRAC_ROOF") then
2458 read(string(indx+1:),*) pv_frac_roof_tbl(1:icate)
2459 else if ( name == "GR_FRAC_ROOF") then
2460 read(string(indx+1:),*) gr_frac_roof_tbl(1:icate)
2461 else if (name == "GR_FLAG") then
2462 read(string(indx+1:),*) gr_flag_tbl
2463 else if (name == "GR_TYPE") then
2464 read(string(indx+1:),*) gr_type_tbl
2466 !end BEP
2467 else
2468 FATAL_ERROR('URBPARM.TBL: Unrecognized NAME = "'//trim(name)//'" in Subr URBAN_PARAM_INIT')
2469 endif
2470 enddo READLOOP
2472 CLOSE(11)
2474 ! Assign a few table values that do not need to come from URBPARM.TBL
2476 Z0HB_TBL = 0.1 * Z0B_TBL
2477 Z0HG_TBL = 0.1 * Z0G_TBL
2479 DO LC = 1, ICATE
2481 ! HGT: Normalized height
2482 HGT_TBL(LC) = ZR_TBL(LC) / ( ROAD_WIDTH(LC) + ROOF_WIDTH(LC) )
2484 ! R: Normalized Roof Width (a.k.a. "building coverage ratio")
2485 R_TBL(LC) = ROOF_WIDTH(LC) / ( ROAD_WIDTH(LC) + ROOF_WIDTH(LC) )
2487 RW_TBL(LC) = 1.0 - R_TBL(LC)
2488 BETR_TBL(LC) = 0.0
2489 BETB_TBL(LC) = 0.0
2490 BETG_TBL(LC) = 0.0
2492 ! The following urban canyon geometry parameters are following Macdonald's (1998) formulations
2494 ! Lambda_P :: Plan areal fraction, which corresponds to R for a 2-d canyon.
2495 ! Lambda_F :: Frontal area index, which corresponds to HGT for a 2-d canyon
2496 ! Cd :: Drag coefficient ( 1.2 from Grimmond and Oke, 1998 )
2497 ! Alpha_macd :: Emperical coefficient ( 4.43 from Macdonald et al., 1998 )
2498 ! Beta_macd :: Correction factor for the drag coefficient ( 1.0 from Macdonald et al., 1998 )
2500 Lambda_P = R_TBL(LC)
2501 Lambda_F = HGT_TBL(LC)
2502 Cd = 1.2
2503 alpha_macd = 4.43
2504 beta_macd = 1.0
2507 ZDC_TBL(LC) = ZR_TBL(LC) * ( 1.0 + ( alpha_macd ** ( -Lambda_P ) ) * ( Lambda_P - 1.0 ) )
2509 Z0C_TBL(LC) = ZR_TBL(LC) * ( 1.0 - ZDC_TBL(LC)/ZR_TBL(LC) ) * &
2510 exp (-(0.5 * beta_macd * Cd / (VonK**2) * ( 1.0-ZDC_TBL(LC)/ZR_TBL(LC) ) * Lambda_F )**(-0.5))
2512 IF (SF_URBAN_PHYSICS == 1) THEN
2513 ! Include roof height variability in Macdonald
2514 ! to parameterize Z0R as a function of ZR_SD (Standard Deviation)
2515 Lambda_FR = SIGMA_ZED_TBL(LC) / ( ROAD_WIDTH(LC) + ROOF_WIDTH(LC) )
2516 Z0R_TBL(LC) = ZR_TBL(LC) * ( 1.0 - ZDC_TBL(LC)/ZR_TBL(LC) ) &
2517 * exp ( -(0.5 * beta_macd * Cd / (VonK**2) &
2518 * ( 1.0-ZDC_TBL(LC)/ZR_TBL(LC) ) * Lambda_FR )**(-0.5))
2519 ENDIF
2521 !
2522 ! Z0HC still one-tenth of Z0C, as before ?
2523 !
2525 Z0HC_TBL(LC) = 0.1 * Z0C_TBL(LC)
2527 !
2528 ! Calculate Sky View Factor:
2529 !
2530 DHGT=HGT_TBL(LC)/100.
2531 HGT=0.
2532 VFWS=0.
2533 HGT=HGT_TBL(LC)-DHGT/2.
2534 do k=1,99
2535 HGT=HGT-DHGT
2536 VFWS=VFWS+0.25*(1.-HGT/SQRT(HGT**2.+RW_TBL(LC)**2.))
2537 end do
2539 VFWS=VFWS/99.
2540 VFWS=VFWS*2.
2542 VFGS=1.-2.*VFWS*HGT_TBL(LC)/RW_TBL(LC)
2543 SVF_TBL(LC)=VFGS
2544 END DO
2546 deallocate(roof_width)
2547 deallocate(road_width)
2549 END SUBROUTINE urban_param_init
2550 !===========================================================================
2551 !
2552 ! subroutine urban_var_init: initialization of urban state variables
2553 !
2554 !===========================================================================
2555 SUBROUTINE urban_var_init(ISURBAN, TSURFACE0_URB,TLAYER0_URB,TDEEP0_URB,IVGTYP, & ! in
2556 ims,ime,jms,jme,kms,kme,num_soil_layers, & ! in
2557 LCZ_1,LCZ_2,LCZ_3,LCZ_4,LCZ_5, &
2558 LCZ_6,LCZ_7,LCZ_8,LCZ_9,LCZ_10,LCZ_11, &
2559 restart,sf_urban_physics, & !in
2560 XXXR_URB2D,XXXB_URB2D,XXXG_URB2D,XXXC_URB2D, & ! inout
2561 TR_URB2D,TB_URB2D,TG_URB2D,TC_URB2D,QC_URB2D, & ! inout
2562 TRL_URB3D,TBL_URB3D,TGL_URB3D, & ! inout
2563 SH_URB2D,LH_URB2D,G_URB2D,RN_URB2D, & ! inout
2564 TS_URB2D, & ! inout
2565 num_urban_ndm, & ! in
2566 urban_map_zrd, & ! in
2567 urban_map_zwd, & ! in
2568 urban_map_gd, & ! in
2569 urban_map_zd, & ! in
2570 urban_map_zdf, & ! in
2571 urban_map_bd, & ! in
2572 urban_map_wd, & ! in
2573 urban_map_gbd, & ! in
2574 urban_map_fbd, & ! in
2575 urban_map_zgrd, & ! in
2576 num_urban_hi, & ! in
2577 TRB_URB4D,TW1_URB4D,TW2_URB4D,TGB_URB4D, & ! inout
2578 TLEV_URB3D,QLEV_URB3D, & ! inout
2579 TW1LEV_URB3D,TW2LEV_URB3D, & ! inout
2580 TGLEV_URB3D,TFLEV_URB3D, & ! inout
2581 SF_AC_URB3D,LF_AC_URB3D,CM_AC_URB3D, & ! inout
2582 SFVENT_URB3D,LFVENT_URB3D, & ! inout
2583 SFWIN1_URB3D,SFWIN2_URB3D, & ! inout
2584 SFW1_URB3D,SFW2_URB3D,SFR_URB3D,SFG_URB3D, & ! inout
2585 EP_PV_URB3D,T_PV_URB3D, & !GRZ
2586 TRV_URB4D,QR_URB4D,QGR_URB3D,TGR_URB3D, & !GRZ
2587 DRAIN_URB4D,DRAINGR_URB3D,SFRV_URB3D, & !GRZ
2588 LFRV_URB3D,DGR_URB3D,DG_URB3D,LFR_URB3D,LFG_URB3D,&!GRZ
2589 SMOIS_URB, &
2590 LP_URB2D,HI_URB2D,LB_URB2D, & ! inout
2591 HGT_URB2D,MH_URB2D,STDH_URB2D, & ! inout
2592 LF_URB2D, & ! inout
2593 CMCR_URB2D,TGR_URB2D,TGRL_URB3D,SMR_URB3D, & ! inout
2594 DRELR_URB2D,DRELB_URB2D,DRELG_URB2D, & ! inout
2595 FLXHUMR_URB2D, FLXHUMB_URB2D, FLXHUMG_URB2D, & ! inout
2596 A_U_BEP,A_V_BEP,A_T_BEP,A_Q_BEP, & ! inout multi-layer urban
2597 A_E_BEP,B_U_BEP,B_V_BEP, & ! inout multi-layer urban
2598 B_T_BEP,B_Q_BEP,B_E_BEP,DLG_BEP, & ! inout multi-layer urban
2599 DL_U_BEP,SF_BEP,VL_BEP, & ! inout multi-layer urban
2600 FRC_URB2D, UTYPE_URB2D,USE_WUDAPT_LCZ) ! inout
2601 IMPLICIT NONE
2603 INTEGER, INTENT(IN) :: ISURBAN, sf_urban_physics,use_wudapt_lcz
2604 INTEGER, INTENT(IN) :: LCZ_1,LCZ_2,LCZ_3,LCZ_4,LCZ_5,LCZ_6,LCZ_7,LCZ_8,LCZ_9,LCZ_10,LCZ_11
2605 INTEGER, INTENT(IN) :: ims,ime,jms,jme,kms,kme,num_soil_layers
2606 INTEGER, INTENT(IN) :: num_urban_ndm
2607 INTEGER, INTENT(IN) :: urban_map_zrd
2608 INTEGER, INTENT(IN) :: urban_map_zwd
2609 INTEGER, INTENT(IN) :: urban_map_gd
2610 INTEGER, INTENT(IN) :: urban_map_zd
2611 INTEGER, INTENT(IN) :: urban_map_zdf
2612 INTEGER, INTENT(IN) :: urban_map_bd
2613 INTEGER, INTENT(IN) :: urban_map_wd
2614 INTEGER, INTENT(IN) :: urban_map_gbd
2615 INTEGER, INTENT(IN) :: urban_map_fbd
2616 INTEGER, INTENT(IN) :: urban_map_zgrd
2617 INTEGER, INTENT(IN) :: num_urban_hi !multi-layer urban
2618 ! INTEGER, INTENT(IN) :: num_roof_layers, num_wall_layers, num_road_layers
2620 REAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN) :: TSURFACE0_URB
2621 REAL, DIMENSION( ims:ime, 1:num_soil_layers, jms:jme ), INTENT(IN) :: TLAYER0_URB
2622 REAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN) :: TDEEP0_URB
2623 INTEGER, DIMENSION( ims:ime, jms:jme ), INTENT(IN) :: IVGTYP
2624 LOGICAL , INTENT(IN) :: restart
2626 REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: TR_URB2D
2627 REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: TB_URB2D
2628 REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: TG_URB2D
2629 REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: TC_URB2D
2630 REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: QC_URB2D
2631 REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: XXXR_URB2D
2632 REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: XXXB_URB2D
2633 REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: XXXG_URB2D
2634 REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: XXXC_URB2D
2636 REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: DRELR_URB2D
2637 REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: DRELB_URB2D
2638 REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: DRELG_URB2D
2639 REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: FLXHUMR_URB2D
2640 REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: FLXHUMB_URB2D
2641 REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: FLXHUMG_URB2D
2642 REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: CMCR_URB2D
2643 REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: TGR_URB2D
2645 ! REAL, DIMENSION(ims:ime, 1:num_roof_layers, jms:jme), INTENT(INOUT) :: TRL_URB3D
2646 ! REAL, DIMENSION(ims:ime, 1:num_wall_layers, jms:jme), INTENT(INOUT) :: TBL_URB3D
2647 ! REAL, DIMENSION(ims:ime, 1:num_road_layers, jms:jme), INTENT(INOUT) :: TGL_URB3D
2648 REAL, DIMENSION(ims:ime, 1:num_soil_layers, jms:jme), INTENT(INOUT) :: TRL_URB3D
2649 REAL, DIMENSION(ims:ime, 1:num_soil_layers, jms:jme), INTENT(INOUT) :: TBL_URB3D
2650 REAL, DIMENSION(ims:ime, 1:num_soil_layers, jms:jme), INTENT(INOUT) :: TGL_URB3D
2651 REAL, DIMENSION(ims:ime, 1:num_soil_layers, jms:jme), INTENT(INOUT) :: TGRL_URB3D
2652 REAL, DIMENSION(ims:ime, 1:num_soil_layers, jms:jme), INTENT(INOUT) :: SMR_URB3D
2654 REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: SH_URB2D
2655 REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: LH_URB2D
2656 REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: G_URB2D
2657 REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: RN_URB2D
2658 REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: TS_URB2D
2660 ! multi-layer UCM variables
2661 REAL, DIMENSION(ims:ime, 1:urban_map_zrd, jms:jme), INTENT(INOUT) :: TRB_URB4D
2662 REAL, DIMENSION(ims:ime, 1:urban_map_zwd, jms:jme), INTENT(INOUT) :: TW1_URB4D
2663 REAL, DIMENSION(ims:ime, 1:urban_map_zwd, jms:jme), INTENT(INOUT) :: TW2_URB4D
2664 REAL, DIMENSION(ims:ime, 1:urban_map_gd , jms:jme), INTENT(INOUT) :: TGB_URB4D
2665 REAL, DIMENSION(ims:ime, 1:urban_map_bd , jms:jme), INTENT(INOUT) :: TLEV_URB3D
2666 REAL, DIMENSION(ims:ime, 1:urban_map_bd , jms:jme), INTENT(INOUT) :: QLEV_URB3D
2667 REAL, DIMENSION(ims:ime, 1:urban_map_wd , jms:jme), INTENT(INOUT) :: TW1LEV_URB3D
2668 REAL, DIMENSION(ims:ime, 1:urban_map_wd , jms:jme), INTENT(INOUT) :: TW2LEV_URB3D
2669 REAL, DIMENSION(ims:ime, 1:urban_map_gbd, jms:jme), INTENT(INOUT) :: TGLEV_URB3D
2670 REAL, DIMENSION(ims:ime, 1:urban_map_fbd, jms:jme), INTENT(INOUT) :: TFLEV_URB3D
2671 REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: LF_AC_URB3D
2672 REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: SF_AC_URB3D
2673 REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: CM_AC_URB3D
2674 REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: SFVENT_URB3D
2675 REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: LFVENT_URB3D
2676 REAL, DIMENSION( ims:ime, 1:urban_map_wd, jms:jme), INTENT(INOUT) :: SFWIN1_URB3D
2677 REAL, DIMENSION( ims:ime, 1:urban_map_wd, jms:jme), INTENT(INOUT) :: SFWIN2_URB3D
2678 REAL, DIMENSION(ims:ime, 1:urban_map_zd , jms:jme), INTENT(INOUT) :: SFW1_URB3D
2679 REAL, DIMENSION(ims:ime, 1:urban_map_zd , jms:jme), INTENT(INOUT) :: SFW2_URB3D
2680 REAL, DIMENSION(ims:ime, 1:urban_map_zdf, jms:jme), INTENT(INOUT) :: SFR_URB3D
2681 REAL, DIMENSION(ims:ime, 1:num_urban_ndm, jms:jme), INTENT(INOUT) :: SFG_URB3D
2682 REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: EP_PV_URB3D!GRZ
2683 REAL, DIMENSION( ims:ime, 1:urban_map_zdf,jms:jme ), INTENT(INOUT) :: T_PV_URB3D!GRZ
2684 REAL, DIMENSION( ims:ime, 1:urban_map_zgrd, jms:jme),INTENT(INOUT) :: TRV_URB4D ! GRZ
2685 REAL, DIMENSION( ims:ime, 1:urban_map_zgrd, jms:jme),INTENT(INOUT) :: QR_URB4D ! GRZ
2686 REAL, DIMENSION( ims:ime,jms:jme), INTENT(INOUT) :: QGR_URB3D ! GRZ
2687 REAL, DIMENSION( ims:ime,jms:jme), INTENT(INOUT) :: TGR_URB3D ! GRZ
2688 REAL, DIMENSION( ims:ime, 1:urban_map_zdf, jms:jme),INTENT(INOUT) :: DRAIN_URB4D !GRZ
2689 REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: DRAINGR_URB3D !GRZ
2690 REAL, DIMENSION( ims:ime, 1:urban_map_zdf, jms:jme),INTENT(INOUT) :: SFRV_URB3D !GRZ
2691 REAL, DIMENSION( ims:ime, 1:urban_map_zdf, jms:jme),INTENT(INOUT) :: LFRV_URB3D ! GRZ
2692 REAL, DIMENSION( ims:ime, 1:urban_map_zdf, jms:jme ),INTENT(INOUT) :: DGR_URB3D !GRZ
2693 REAL, DIMENSION( ims:ime, 1:num_urban_ndm, jms:jme ),INTENT(INOUT) :: DG_URB3D !GRZ
2694 REAL, DIMENSION( ims:ime, 1:urban_map_zdf, jms:jme ),INTENT(INOUT) :: LFR_URB3D !GRZ
2695 REAL, DIMENSION( ims:ime, 1:num_urban_ndm, jms:jme ), INTENT(INOUT) :: LFG_URB3D !GRZ
2696 REAL, DIMENSION( ims:ime, 1:num_soil_layers, jms:jme ), INTENT(IN) ::SMOIS_URB
2697 REAL, DIMENSION( ims:ime,1:num_urban_hi , jms:jme), INTENT(INOUT) :: HI_URB2D
2698 REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: LP_URB2D
2699 REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: LB_URB2D
2700 REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: HGT_URB2D
2701 REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: MH_URB2D
2702 REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: STDH_URB2D
2703 REAL, DIMENSION( ims:ime, 4,jms:jme ), INTENT(INOUT) :: LF_URB2D
2704 REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(INOUT) :: A_U_BEP
2705 REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(INOUT) :: A_V_BEP
2706 REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(INOUT) :: A_T_BEP
2707 REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(INOUT) :: A_Q_BEP
2708 REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(INOUT) :: A_E_BEP
2709 REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(INOUT) :: B_U_BEP
2710 REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(INOUT) :: B_V_BEP
2711 REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(INOUT) :: B_T_BEP
2712 REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(INOUT) :: B_Q_BEP
2713 REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(INOUT) :: B_E_BEP
2714 REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(INOUT) :: VL_BEP
2715 REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(INOUT) :: DLG_BEP
2716 REAL, DIMENSION(ims:ime, kms:kme,jms:jme),INTENT(INOUT) :: SF_BEP
2717 REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(INOUT) :: DL_U_BEP
2718 !
2719 REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: FRC_URB2D
2720 INTEGER, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: UTYPE_URB2D
2721 INTEGER :: UTYPE_URB
2722 !FS
2723 INTEGER :: SWITCH_URB
2725 INTEGER :: I,J,K,CHECK
2727 CHECK = 0
2729 DO I=ims,ime
2730 DO J=jms,jme
2732 ! XXXR_URB2D(I,J)=0.
2733 ! XXXB_URB2D(I,J)=0.
2734 ! XXXG_URB2D(I,J)=0.
2735 ! XXXC_URB2D(I,J)=0.
2737 SH_URB2D(I,J)=0.
2738 LH_URB2D(I,J)=0.
2739 G_URB2D(I,J)=0.
2740 RN_URB2D(I,J)=0.
2742 !m
2743 !FS FRC_URB2D(I,J)=0.
2744 UTYPE_URB2D(I,J)=0
2746 distributed_aerodynamics_check: IF (distributed_aerodynamics_option) THEN
2747 IF (IVGTYP(I, J) == ISURBAN) THEN
2748 UTYPE_URB2D(I, J) = 2
2749 ELSE
2750 UTYPE_URB2D(I, J) = 0
2751 END IF
2752 ELSE
2753 SWITCH_URB=1
2755 IF( IVGTYP(I,J) == ISURBAN) THEN
2756 IF(use_wudapt_lcz==0) THEN
2757 UTYPE_URB2D(I,J) = 2 ! for default. high-intensity
2758 ELSE
2759 UTYPE_URB2D(I,J) = 5 ! for default. high-intensity
2760 ENDIF
2761 ELSE IF( IVGTYP(I,J) == LCZ_1) THEN
2762 UTYPE_URB2D(I,J) = 1
2763 ELSE IF( IVGTYP(I,J) == LCZ_2) THEN
2764 UTYPE_URB2D(I,J) = 2
2765 ELSE IF( IVGTYP(I,J) == LCZ_3) THEN
2766 UTYPE_URB2D(I,J) = 3
2767 ELSE IF( IVGTYP(I,J) == LCZ_4) THEN
2768 UTYPE_URB2D(I,J) = 4
2769 ELSE IF( IVGTYP(I,J) == LCZ_5) THEN
2770 UTYPE_URB2D(I,J) = 5
2771 ELSE IF( IVGTYP(I,J) == LCZ_6) THEN
2772 UTYPE_URB2D(I,J) = 6
2773 ELSE IF( IVGTYP(I,J) == LCZ_7) THEN
2774 UTYPE_URB2D(I,J) = 7
2775 ELSE IF( IVGTYP(I,J) == LCZ_8) THEN
2776 UTYPE_URB2D(I,J) = 8
2777 ELSE IF( IVGTYP(I,J) == LCZ_9) THEN
2778 UTYPE_URB2D(I,J) = 9
2779 ELSE IF( IVGTYP(I,J) == LCZ_10) THEN
2780 UTYPE_URB2D(I,J) = 10
2781 ELSE IF( IVGTYP(I,J) == LCZ_11) THEN
2782 UTYPE_URB2D(I,J) = 11
2783 ELSE
2784 SWITCH_URB=0
2785 ENDIF
2787 IF (SWITCH_URB == 1) THEN
2788 UTYPE_URB = UTYPE_URB2D(I,J) ! for default. high-intensity
2789 IF (HGT_URB2D(I,J)>0.) THEN
2790 CONTINUE
2791 ELSE
2792 WRITE(mesg,*) 'USING DEFAULT URBAN MORPHOLOGY'
2793 WRITE_MESSAGE(mesg)
2794 LP_URB2D(I,J)=0.
2795 LB_URB2D(I,J)=0.
2796 HGT_URB2D(I,J)=0.
2797 IF ( sf_urban_physics == 1 ) THEN
2798 MH_URB2D(I,J)=0.
2799 STDH_URB2D(I,J)=0.
2800 DO K=1,4
2801 LF_URB2D(I,K,J)=0.
2802 ENDDO
2803 ELSE IF ( ( sf_urban_physics == 2 ) .or. ( sf_urban_physics == 3 ) ) THEN
2804 DO K=1,num_urban_hi
2805 HI_URB2D(I,K,J)=0.
2806 ENDDO
2807 ENDIF
2808 ENDIF
2809 IF (FRC_URB2D(I,J)>0.and.FRC_URB2D(I,J)<=1.) THEN
2810 CONTINUE
2811 ELSE
2812 WRITE(mesg,*) 'WARNING, FRC_URB2D = 0 BUT IVGTYP IS URBAN'
2813 WRITE_MESSAGE(mesg)
2814 WRITE(mesg,*) 'WARNING, THE URBAN FRACTION WILL BE READ FROM URBPARM.TBL'
2815 WRITE_MESSAGE(mesg)
2816 FRC_URB2D(I,J) = FRC_URB_TBL(UTYPE_URB)
2817 ENDIF
2818 ELSE
2819 FRC_URB2D(I,J)=0.
2820 LP_URB2D(I,J)=0.
2821 LB_URB2D(I,J)=0.
2822 HGT_URB2D(I,J)=0.
2823 IF ( sf_urban_physics == 1 ) THEN
2824 MH_URB2D(I,J)=0.
2825 STDH_URB2D(I,J)=0.
2826 DO K=1,4
2827 LF_URB2D(I,K,J)=0.
2828 ENDDO
2829 ELSE IF ( ( sf_urban_physics == 2 ) .or. ( sf_urban_physics == 3 ) ) THEN
2830 DO K=1,num_urban_hi
2831 HI_URB2D(I,K,J)=0.
2832 ENDDO
2833 ENDIF
2834 ENDIF
2835 END IF distributed_aerodynamics_check
2838 QC_URB2D(I,J)=0.01
2840 IF (.not.restart) THEN
2842 XXXR_URB2D(I,J)=0.
2843 XXXB_URB2D(I,J)=0.
2844 XXXG_URB2D(I,J)=0.
2845 XXXC_URB2D(I,J)=0.
2847 IF ( sf_urban_physics == 1 ) THEN
2848 DRELR_URB2D(I,J) = 0.
2849 DRELB_URB2D(I,J) = 0.
2850 DRELG_URB2D(I,J) = 0.
2851 FLXHUMR_URB2D(I,J) = 0.
2852 FLXHUMB_URB2D(I,J) = 0.
2853 FLXHUMG_URB2D(I,J) = 0.
2854 CMCR_URB2D(I,J) = 0.
2855 TGR_URB2D(I,J)=TSURFACE0_URB(I,J)+0.
2856 ENDIF
2858 TC_URB2D(I,J)=TSURFACE0_URB(I,J)+0.
2859 TR_URB2D(I,J)=TSURFACE0_URB(I,J)+0.
2860 TB_URB2D(I,J)=TSURFACE0_URB(I,J)+0.
2861 TG_URB2D(I,J)=TSURFACE0_URB(I,J)+0.
2862 !
2863 TS_URB2D(I,J)=TSURFACE0_URB(I,J)+0.
2865 ! DO K=1,num_roof_layers
2866 ! DO K=1,num_soil_layers
2867 ! TRL_URB3D(I,1,J)=TLAYER0_URB(I,1,J)+0.
2868 ! TRL_URB3D(I,2,J)=TLAYER0_URB(I,2,J)+0.
2869 ! TRL_URB3D(I,3,J)=TLAYER0_URB(I,3,J)+0.
2870 ! TRL_URB3D(I,4,J)=TLAYER0_URB(I,4,J)+0.
2872 TRL_URB3D(I,1,J)=TLAYER0_URB(I,1,J)+0.
2873 TRL_URB3D(I,2,J)=0.5*(TLAYER0_URB(I,1,J)+TLAYER0_URB(I,2,J))
2874 TRL_URB3D(I,3,J)=TLAYER0_URB(I,2,J)+0.
2875 TRL_URB3D(I,4,J)=TLAYER0_URB(I,2,J)+(TLAYER0_URB(I,3,J)-TLAYER0_URB(I,2,J))*0.29
2877 IF ( sf_urban_physics == 1 ) THEN
2878 TGRL_URB3D(I,1,J)=TLAYER0_URB(I,1,J)+0.
2879 TGRL_URB3D(I,2,J)=0.5*(TLAYER0_URB(I,1,J)+TLAYER0_URB(I,2,J))
2880 TGRL_URB3D(I,3,J)=TLAYER0_URB(I,2,J)+0.
2881 TGRL_URB3D(I,4,J)=TLAYER0_URB(I,2,J)+(TLAYER0_URB(I,3,J)-TLAYER0_URB(I,2,J))*0.29
2883 SMR_URB3D(I,1,J)=0.2
2884 SMR_URB3D(I,2,J)=0.2
2885 SMR_URB3D(I,3,J)=0.2
2886 SMR_URB3D(I,4,J)=0.
2887 ENDIF
2888 ! END DO
2890 ! DO K=1,num_wall_layers
2891 ! DO K=1,num_soil_layers
2892 !m TBL_URB3D(I,1,J)=TLAYER0_URB(I,1,J)+0.
2893 !m TBL_URB3D(I,2,J)=TLAYER0_URB(I,2,J)+0.
2894 !m TBL_URB3D(I,3,J)=TLAYER0_URB(I,3,J)+0.
2895 !m TBL_URB3D(I,4,J)=TLAYER0_URB(I,4,J)+0.
2897 TBL_URB3D(I,1,J)=TLAYER0_URB(I,1,J)+0.
2898 TBL_URB3D(I,2,J)=0.5*(TLAYER0_URB(I,1,J)+TLAYER0_URB(I,2,J))
2899 TBL_URB3D(I,3,J)=TLAYER0_URB(I,2,J)+0.
2900 TBL_URB3D(I,4,J)=TLAYER0_URB(I,2,J)+(TLAYER0_URB(I,3,J)-TLAYER0_URB(I,2,J))*0.29
2902 ! END DO
2904 ! DO K=1,num_road_layers
2905 DO K=1,num_soil_layers
2906 TGL_URB3D(I,K,J)=TLAYER0_URB(I,K,J)+0.
2907 END DO
2909 ! multi-layer urban
2910 IF((SF_URBAN_PHYSICS.eq.2).OR.(SF_URBAN_PHYSICS.eq.3)) THEN
2911 IF (UTYPE_URB2D(I,J) > 0) THEN
2912 TRB_URB4D(I,:,J)=TBLEND_TBL(UTYPE_URB2D(I,J))
2913 TW1_URB4D(I,:,J)=TBLEND_TBL(UTYPE_URB2D(I,J))
2914 TW2_URB4D(I,:,J)=TBLEND_TBL(UTYPE_URB2D(I,J))
2915 ELSE
2916 TRB_URB4D(I,:,J)=tlayer0_urb(I,1,J)
2917 TW1_URB4D(I,:,J)=tlayer0_urb(I,1,J)
2918 TW2_URB4D(I,:,J)=tlayer0_urb(I,1,J)
2919 ENDIF
2920 TGB_URB4D(I,:,J)=tlayer0_urb(I,1,J)
2921 SFW1_URB3D(I,:,J)=0.
2922 SFW2_URB3D(I,:,J)=0.
2923 SFR_URB3D(I,:,J)=0.
2924 SFG_URB3D(I,:,J)=0.
2926 ENDIF
2928 if (SF_URBAN_PHYSICS.EQ.3) then
2929 LF_AC_URB3D(I,J)=0.
2930 SF_AC_URB3D(I,J)=0.
2931 CM_AC_URB3D(I,J)=0.
2932 SFVENT_URB3D(I,J)=0.
2933 LFVENT_URB3D(I,J)=0.
2934 EP_PV_URB3D(I,J)=0.
2935 T_PV_URB3D(I,:,J)=tlayer0_urb(I,1,J)
2936 TGR_URB3D(I,J)=tlayer0_urb(I,1,J)
2937 QR_URB4D(I,:,J)=SMOIS_URB(I,1,J)
2938 DRAIN_URB4D(I,:,J)=0. !GRZ
2939 SFRV_URB3D(I,:,J)=0. !GRZ
2940 LFRV_URB3D(I,:,J)=0. !GRZ
2941 DGR_URB3D(I,:,J)=0. !GRZ
2942 DG_URB3D(I,:,J)=0.
2943 LFR_URB3D(I,:,J)=0.
2944 LFG_URB3D(I,:,J)=0.
2945 QGR_URB3D(I,J)=SMOIS_URB(I,1,J) !GRZ
2946 TGR_URB3D(I,J)=0.
2947 DRAINGR_URB3D(I,J)=0. !GRZ
2949 IF (UTYPE_URB2D(I,J) > 0) THEN
2950 TRV_URB4D(I,:,J)=TBLEND_TBL(UTYPE_URB2D(I,J))
2951 ELSE
2952 TRV_URB4D(I,:,J)=tlayer0_urb(I,1,J) !GRZ
2953 ENDIF
2955 TLEV_URB3D(I,:,J)=tlayer0_urb(I,1,J)
2956 TW1LEV_URB3D(I,:,J)=tlayer0_urb(I,1,J)
2957 TW2LEV_URB3D(I,:,J)=tlayer0_urb(I,1,J)
2958 TGLEV_URB3D(I,:,J)=tlayer0_urb(I,1,J)
2959 TFLEV_URB3D(I,:,J)=tlayer0_urb(I,1,J)
2960 QLEV_URB3D(I,:,J)=0.01
2961 SFWIN1_URB3D(I,:,J)=0.
2962 SFWIN2_URB3D(I,:,J)=0.
2963 !rm LF_AC_URB3D(I,J)=0.
2964 !rm SF_AC_URB3D(I,J)=0.
2965 !rm CM_AC_URB3D(I,J)=0.
2966 !rm SFVENT_URB3D(I,J)=0.
2967 !rm LFVENT_URB3D(I,J)=0.
2969 endif
2971 ! IF( sf_urban_physics .EQ. 2 )THEN
2972 IF((SF_URBAN_PHYSICS.eq.2).OR.(SF_URBAN_PHYSICS.eq.3)) THEN
2973 DO K= KMS,KME
2974 SF_BEP(I,K,J)=1.
2975 VL_BEP(I,K,J)=1.
2976 A_U_BEP(I,K,J)=0.
2977 A_V_BEP(I,K,J)=0.
2978 A_T_BEP(I,K,J)=0.
2979 A_E_BEP(I,K,J)=0.
2980 A_Q_BEP(I,K,J)=0.
2981 B_U_BEP(I,K,J)=0.
2982 B_V_BEP(I,K,J)=0.
2983 B_T_BEP(I,K,J)=0.
2984 B_E_BEP(I,K,J)=0.
2985 B_Q_BEP(I,K,J)=0.
2986 DLG_BEP(I,K,J)=0.
2987 DL_U_BEP(I,K,J)=0.
2988 END DO
2989 ENDIF !sf_urban_physics=2
2990 ENDIF !restart
2993 IF (CHECK.EQ.0)THEN
2994 IF(IVGTYP(I,J).EQ.1)THEN
2995 write(mesg,*) 'Sample of Urban settings'
2996 WRITE_MESSAGE(mesg)
2997 write(mesg,*) 'TSURFACE0_URB',TSURFACE0_URB(I,J)
2998 WRITE_MESSAGE(mesg)
2999 write(mesg,*) 'TDEEP0_URB', TDEEP0_URB(I,J)
3000 WRITE_MESSAGE(mesg)
3001 write(mesg,*) 'IVGTYP',IVGTYP(I,J)
3002 WRITE_MESSAGE(mesg)
3003 write(mesg,*) 'TR_URB2D',TR_URB2D(I,J)
3004 WRITE_MESSAGE(mesg)
3005 write(mesg,*) 'TB_URB2D',TB_URB2D(I,J)
3006 WRITE_MESSAGE(mesg)
3007 write(mesg,*) 'TG_URB2D',TG_URB2D(I,J)
3008 WRITE_MESSAGE(mesg)
3009 write(mesg,*) 'TC_URB2D',TC_URB2D(I,J)
3010 WRITE_MESSAGE(mesg)
3011 write(mesg,*) 'QC_URB2D',QC_URB2D(I,J)
3012 WRITE_MESSAGE(mesg)
3013 write(mesg,*) 'XXXR_URB2D',XXXR_URB2D(I,J)
3014 WRITE_MESSAGE(mesg)
3015 write(mesg,*) 'SH_URB2D',SH_URB2D(I,J)
3016 WRITE_MESSAGE(mesg)
3017 write(mesg,*) 'LH_URB2D',LH_URB2D(I,J)
3018 WRITE_MESSAGE(mesg)
3019 write(mesg,*) 'G_URB2D',G_URB2D(I,J)
3020 WRITE_MESSAGE(mesg)
3021 write(mesg,*) 'RN_URB2D',RN_URB2D(I,J)
3022 WRITE_MESSAGE(mesg)
3023 write(mesg,*) 'TS_URB2D',TS_URB2D(I,J)
3024 WRITE_MESSAGE(mesg)
3025 write(mesg,*) 'LF_AC_URB3D', LF_AC_URB3D(I,J)
3026 WRITE_MESSAGE(mesg)
3027 write(mesg,*) 'SF_AC_URB3D', SF_AC_URB3D(I,J)
3028 WRITE_MESSAGE(mesg)
3029 write(mesg,*) 'CM_AC_URB3D', CM_AC_URB3D(I,J)
3030 WRITE_MESSAGE(mesg)
3031 write(mesg,*) 'SFVENT_URB3D', SFVENT_URB3D(I,J)
3032 WRITE_MESSAGE(mesg)
3033 write(mesg,*) 'LFVENT_URB3D', LFVENT_URB3D(I,J)
3034 WRITE_MESSAGE(mesg)
3035 write(mesg,*) 'FRC_URB2D', FRC_URB2D(I,J)
3036 WRITE_MESSAGE(mesg)
3037 write(mesg,*) 'UTYPE_URB2D', UTYPE_URB2D(I,J)
3038 WRITE_MESSAGE(mesg)
3039 write(mesg,*) 'I',I,'J',J
3040 WRITE_MESSAGE(mesg)
3041 write(mesg,*) 'num_urban_hi', num_urban_hi
3042 WRITE_MESSAGE(mesg)
3043 CHECK = 1
3044 END IF
3045 END IF
3047 END DO
3048 END DO
3049 RETURN
3050 END SUBROUTINE urban_var_init
3051 !===========================================================================
3052 !
3053 ! force_restore
3054 !
3055 !===========================================================================
3056 SUBROUTINE force_restore(CAP,AKS,DELT,S,R,H,LE,TSLEND,TSP,TS)
3058 REAL, INTENT(IN) :: CAP,AKS,DELT,S,R,H,LE,TSLEND,TSP
3059 REAL, INTENT(OUT) :: TS
3060 REAL :: C1,C2
3062 C2=24.*3600./2./3.14159
3063 C1=SQRT(0.5*C2*CAP*AKS)
3065 TS = TSP + DELT*( (S+R-H-LE)/C1 -(TSP-TSLEND)/C2 )
3067 END SUBROUTINE force_restore
3068 !===========================================================================
3069 !
3070 ! bisection (not used)
3071 !
3072 !==============================================================================
3073 SUBROUTINE bisection(TSP,PS,S,EPS,RX,SIG,RHO,CP,CH,UA,QA,TA,EL,BET,AKS,TSL,DZ,TS)
3075 REAL, INTENT(IN) :: TSP,PS,S,EPS,RX,SIG,RHO,CP,CH,UA,QA,TA,EL,BET,AKS,TSL,DZ
3076 REAL, INTENT(OUT) :: TS
3077 REAL :: ES,QS0,R,H,ELE,G0,F1,F
3079 TS1 = TSP - 5.
3080 TS2 = TSP + 5.
3082 DO ITERATION = 1,22
3084 ES=6.11*EXP( (2.5*10.**6./461.51)*(TS1-273.15)/(273.15*TS1) )
3085 QS0=0.622*ES/(PS-0.378*ES)
3086 R=EPS*(RX-SIG*(TS1**4.)/60.)
3087 H=RHO*CP*CH*UA*(TS1-TA)*100.
3088 ELE=RHO*EL*CH*UA*BET*(QS0-QA)*100.
3089 G0=AKS*(TS1-TSL)/(DZ/2.)
3090 F1= S + R - H - ELE - G0
3092 TS=0.5*(TS1+TS2)
3094 ES=6.11*EXP( (2.5*10.**6./461.51)*(TS-273.15)/(273.15*TS) )
3095 QS0=0.622*ES/(PS-0.378*ES)
3096 R=EPS*(RX-SIG*(TS**4.)/60.)
3097 H=RHO*CP*CH*UA*(TS-TA)*100.
3098 ELE=RHO*EL*CH*UA*BET*(QS0-QA)*100.
3099 G0=AKS*(TS-TSL)/(DZ/2.)
3100 F = S + R - H - ELE - G0
3102 IF (F1*F > 0.0) THEN
3103 TS1=TS
3104 ELSE
3105 TS2=TS
3106 END IF
3108 END DO
3110 RETURN
3111 END SUBROUTINE bisection
3112 !===========================================================================
3114 SUBROUTINE SFCDIF_URB (ZLM,Z0,THZ0,THLM,SFCSPD,AKANDA,AKMS,AKHS,RLMO,CD,ZT_OUT,VEGFRAC)
3116 ! ----------------------------------------------------------------------
3117 ! SUBROUTINE SFCDIF_URB (Urban version of SFCDIF_off)
3118 ! ----------------------------------------------------------------------
3119 ! CALCULATE SURFACE LAYER EXCHANGE COEFFICIENTS VIA ITERATIVE PROCESS.
3120 ! SEE CHEN ET AL (1997, BLM)
3121 ! ----------------------------------------------------------------------
3123 IMPLICIT NONE
3124 REAL WWST, WWST2, G, VKRM, EXCM, BETA, BTG, ELFC, WOLD, WNEW
3125 REAL PIHF, EPSU2, EPSUST, EPSIT, EPSA, ZTMIN, ZTMAX, HPBL, &
3126 & SQVISC
3127 REAL RIC, RRIC, FHNEU, RFC,RLMO_THR, RFAC, ZZ, PSLMU, PSLMS, PSLHU, &
3128 & PSLHS
3129 REAL XX, PSPMU, YY, PSPMS, PSPHU, PSPHS, ZLM, Z0, THZ0, THLM
3130 REAL SFCSPD, AKANDA, AKMS, AKHS, ZU, ZT, RDZ, CXCH
3131 REAL DTHV, DU2, BTGH, WSTAR2, USTAR, ZSLU, ZSLT, RLOGU, RLOGT
3132 REAL RLMO, ZETALT, ZETALU, ZETAU, ZETAT, XLU4, XLT4, XU4, XT4
3133 REAL, INTENT(OUT), OPTIONAL :: ZT_OUT
3134 REAL, INTENT(IN), OPTIONAL :: VEGFRAC
3135 !CC ......REAL ZTFC
3137 REAL XLU, XLT, XU, XT, PSMZ, SIMM, PSHZ, SIMH, USTARK, RLMN, &
3138 & RLMA
3140 INTEGER ITRMX, ILECH, ITR
3141 REAL, INTENT(OUT) :: CD
3142 PARAMETER &
3143 & (WWST = 1.2,WWST2 = WWST * WWST,G = 9.8,VKRM = 0.40, &
3144 & EXCM = 0.001 &
3145 & ,BETA = 1./270.,BTG = BETA * G,ELFC = VKRM * BTG &
3146 & ,WOLD =.15,WNEW = 1. - WOLD,ITRMX = 05, &
3147 & PIHF = 3.14159265/2.)
3148 PARAMETER &
3149 & (EPSU2 = 1.E-4,EPSUST = 0.07,EPSIT = 1.E-4,EPSA = 1.E-8 &
3150 & ,ZTMIN = -5.,ZTMAX = 1.,HPBL = 1000.0 &
3151 & ,SQVISC = 258.2)
3152 PARAMETER &
3153 & (RIC = 0.183,RRIC = 1.0/ RIC,FHNEU = 0.8,RFC = 0.191 &
3154 & ,RLMO_THR = 0.001,RFAC = RIC / (FHNEU * RFC * RFC))
3156 ! ----------------------------------------------------------------------
3157 ! NOTE: THE TWO CODE BLOCKS BELOW DEFINE FUNCTIONS
3158 ! ----------------------------------------------------------------------
3159 ! LECH'S SURFACE FUNCTIONS
3160 ! ----------------------------------------------------------------------
3161 PSLMU (ZZ)= -0.96* log (1.0-4.5* ZZ)
3162 PSLMS (ZZ)= (ZZ/RFC) -2.076* (1. -1./ (ZZ +1.))
3163 PSLHU (ZZ)= -0.96* log (1.0-4.5* ZZ)
3165 ! ----------------------------------------------------------------------
3166 ! PAULSON'S SURFACE FUNCTIONS
3167 ! ----------------------------------------------------------------------
3168 PSLHS (ZZ)= ZZ * RFAC -2.076* (1. - exp(-1.2 * ZZ))
3169 PSPMU (XX)= -2.* log ( (XX +1.)*0.5) - log ( (XX * XX +1.)*0.5) &
3170 & +2.* ATAN (XX) &
3171 &- PIHF
3172 PSPMS (YY)= 5.* YY
3173 PSPHU (XX)= -2.* log ( (XX * XX +1.)*0.5)
3175 ! ----------------------------------------------------------------------
3176 ! THIS ROUTINE SFCDIF CAN HANDLE BOTH OVER OPEN WATER (SEA, OCEAN) AND
3177 ! OVER SOLID SURFACE (LAND, SEA-ICE).
3178 ! ----------------------------------------------------------------------
3179 PSPHS (YY)= 5.* YY
3181 ! ----------------------------------------------------------------------
3182 ! ZTFC: RATIO OF ZOH/ZOM LESS OR EQUAL THAN 1
3183 ! C......ZTFC=0.1
3184 ! CZIL: CONSTANT C IN Zilitinkevich, S. S.1995,:NOTE ABOUT ZT
3185 ! ----------------------------------------------------------------------
3186 ILECH = 0
3188 ! ----------------------------------------------------------------------
3189 ! ZILFC = - CZIL * VKRM * SQVISC
3190 ! C.......ZT=Z0*ZTFC
3191 ZU = Z0
3192 RDZ = 1./ ZLM
3193 CXCH = EXCM * RDZ
3194 DTHV = THLM - THZ0
3196 ! ----------------------------------------------------------------------
3197 ! BELJARS CORRECTION OF USTAR
3198 ! ----------------------------------------------------------------------
3199 DU2 = MAX (SFCSPD * SFCSPD,EPSU2)
3200 !cc If statements to avoid TANGENT LINEAR problems near zero
3201 BTGH = BTG * HPBL
3202 IF (BTGH * AKHS * DTHV .ne. 0.0) THEN
3203 WSTAR2 = WWST2* ABS (BTGH * AKHS * DTHV)** (2./3.)
3204 ELSE
3205 WSTAR2 = 0.0
3206 END IF
3208 ! ----------------------------------------------------------------------
3209 ! ZILITINKEVITCH APPROACH FOR ZT
3210 ! ----------------------------------------------------------------------
3211 USTAR = MAX (SQRT (AKMS * SQRT (DU2+ WSTAR2)),EPSUST)
3213 ! ----------------------------------------------------------------------
3214 ! KCL/TL Try Kanda approach instead (Kanda et al. 2007, JAMC)
3215 ! ZT = EXP (ZILFC * SQRT (USTAR * Z0))* Z0
3216 IF (PRESENT(VEGFRAC)) THEN
3217 ! Kawai et al. (2009) JAMC
3218 ZT = EXP (2.0-(AKANDA-0.9*VEGFRAC**0.29)*(SQVISC**2 * USTAR * Z0)**0.25)* Z0
3219 ELSE
3220 ZT = EXP (2.0-AKANDA*(SQVISC**2 * USTAR * Z0)**0.25)* Z0
3221 END IF
3223 ZSLU = ZLM + ZU
3225 ZSLT = ZLM + ZT
3226 RLOGU = log (ZSLU / ZU)
3228 RLOGT = log (ZSLT / ZT)
3230 RLMO = ELFC * AKHS * DTHV / USTAR **3
3231 ! ----------------------------------------------------------------------
3232 ! 1./MONIN-OBUKKHOV LENGTH-SCALE
3233 ! ----------------------------------------------------------------------
3234 DO ITR = 1,ITRMX
3235 ZETALT = MAX (ZSLT * RLMO,ZTMIN)
3236 RLMO = ZETALT / ZSLT
3237 ZETALU = ZSLU * RLMO
3238 ZETAU = ZU * RLMO
3240 ZETAT = ZT * RLMO
3241 IF (ILECH .eq. 0) THEN
3242 IF (RLMO .lt. 0.0)THEN
3243 XLU4 = 1. -16.* ZETALU
3244 XLT4 = 1. -16.* ZETALT
3245 XU4 = 1. -16.* ZETAU
3247 XT4 = 1. -16.* ZETAT
3248 XLU = SQRT (SQRT (XLU4))
3249 XLT = SQRT (SQRT (XLT4))
3250 XU = SQRT (SQRT (XU4))
3252 XT = SQRT (SQRT (XT4))
3254 PSMZ = PSPMU (XU)
3255 SIMM = PSPMU (XLU) - PSMZ + RLOGU
3256 PSHZ = PSPHU (XT)
3257 SIMH = PSPHU (XLT) - PSHZ + RLOGT
3258 ELSE
3259 ZETALU = MIN (ZETALU,ZTMAX)
3260 ZETALT = MIN (ZETALT,ZTMAX)
3261 PSMZ = PSPMS (ZETAU)
3262 SIMM = PSPMS (ZETALU) - PSMZ + RLOGU
3263 PSHZ = PSPHS (ZETAT)
3264 SIMH = PSPHS (ZETALT) - PSHZ + RLOGT
3265 END IF
3266 ! ----------------------------------------------------------------------
3267 ! LECH'S FUNCTIONS
3268 ! ----------------------------------------------------------------------
3269 ELSE
3270 IF (RLMO .lt. 0.)THEN
3271 PSMZ = PSLMU (ZETAU)
3272 SIMM = PSLMU (ZETALU) - PSMZ + RLOGU
3273 PSHZ = PSLHU (ZETAT)
3274 SIMH = PSLHU (ZETALT) - PSHZ + RLOGT
3275 ELSE
3276 ZETALU = MIN (ZETALU,ZTMAX)
3277 ZETALT = MIN (ZETALT,ZTMAX)
3278 PSMZ = PSLMS (ZETAU)
3279 SIMM = PSLMS (ZETALU) - PSMZ + RLOGU
3280 PSHZ = PSLHS (ZETAT)
3281 SIMH = PSLHS (ZETALT) - PSHZ + RLOGT
3282 END IF
3283 ! ----------------------------------------------------------------------
3284 ! BELJAARS CORRECTION FOR USTAR
3285 ! ----------------------------------------------------------------------
3286 END IF
3287 USTAR = MAX (SQRT (AKMS * SQRT (DU2+ WSTAR2)),EPSUST)
3288 !KCL/TL
3289 !ZT = EXP (ZILFC * SQRT (USTAR * Z0))* Z0
3290 IF (PRESENT(VEGFRAC)) THEN
3291 ! Kawai et al. (2009) JAMC
3292 ZT = EXP (2.0-(AKANDA-0.9*VEGFRAC**0.29)*(SQVISC**2 * USTAR * Z0)**0.25)* Z0
3293 ELSE
3294 ZT = EXP (2.0-AKANDA*(SQVISC**2 * USTAR * Z0)**0.25)* Z0
3295 END IF
3296 ZSLT = ZLM + ZT
3297 RLOGT = log (ZSLT / ZT)
3298 USTARK = USTAR * VKRM
3299 AKMS = MAX (USTARK / SIMM,CXCH)
3300 AKHS = MAX (USTARK / SIMH,CXCH)
3301 !
3302 IF (BTGH * AKHS * DTHV .ne. 0.0) THEN
3303 WSTAR2 = WWST2* ABS (BTGH * AKHS * DTHV)** (2./3.)
3304 ELSE
3305 WSTAR2 = 0.0
3306 END IF
3307 !-----------------------------------------------------------------------
3308 RLMN = ELFC * AKHS * DTHV / USTAR **3
3309 !-----------------------------------------------------------------------
3310 ! IF(ABS((RLMN-RLMO)/RLMA).LT.EPSIT) GO TO 110
3311 !-----------------------------------------------------------------------
3312 RLMA = RLMO * WOLD+ RLMN * WNEW
3313 !-----------------------------------------------------------------------
3314 RLMO = RLMA
3316 END DO
3318 CD = USTAR*USTAR/SFCSPD**2
3320 IF (PRESENT(ZT_OUT)) ZT_OUT = ZT
3321 ! ----------------------------------------------------------------------
3322 END SUBROUTINE SFCDIF_URB
3323 ! ----------------------------------------------------------------------
3324 !===========================================================================
3325 ! DIREVAP
3326 ! CALCULATE DIRECT SOIL EVAPORATION
3327 !===========================================================================
3328 SUBROUTINE DIREVAP (EDIR,ETP,SMC,SHDFAC,SMCMAX,SMCDRY,FXEXP)
3330 REAL, INTENT(IN) :: ETP,SMC,SHDFAC,SMCMAX,SMCDRY,FXEXP
3331 REAL, INTENT(OUT) :: EDIR
3332 REAL :: FX, SRATIO
3334 ! ----------------------------------------------------------------------
3335 ! FX > 1 REPRESENTS DEMAND CONTROL
3336 ! FX < 1 REPRESENTS FLUX CONTROL
3337 ! ----------------------------------------------------------------------
3338 SRATIO = (SMC - SMCDRY) / (SMCMAX - SMCDRY)
3339 IF (SRATIO > 0.) THEN
3340 FX = SRATIO**FXEXP
3341 FX = MAX ( MIN ( FX, 1. ) ,0. )
3342 ELSE
3343 FX = 0.
3344 ENDIF
3345 EDIR = FX * ( 1.0- SHDFAC ) * ETP * 0.001
3347 END SUBROUTINE DIREVAP
3348 !===========================================================================
3349 ! TRANSP
3350 ! CALCULATE EVAPOTRANSPIRATION FOR VEGETATIO SURFACE
3351 !===========================================================================
3353 SUBROUTINE TRANSP (ETT,ET,EC,SHDFAC,ETP1,CMC,CFACTR,CMCMAX,LAI,RSMIN,RSMAX,RGL,SX, &
3354 TS,TA,QA,SMC,SMCWLT,SMCREF,CPP,PS,CH,EPSV,DELT, NROOT,NSOIL, &
3355 DZVR, ZSOIL, HS)
3356 INTEGER, INTENT(IN) :: NROOT, NSOIL
3357 REAL, INTENT(IN) :: SHDFAC,ETP1,CMC,CFACTR,CMCMAX,LAI,RSMIN,RSMAX,RGL,SX,TA
3358 REAL, INTENT(IN) :: TS,QA, SMCWLT, SMCREF, CPP, PS,CH, EPSV, DELT, HS
3359 REAL, DIMENSION(1:NSOIL), INTENT(IN) :: ZSOIL, DZVR, SMC
3360 REAL, DIMENSION(1:NSOIL), INTENT(INOUT):: ET
3361 REAL, INTENT(OUT) :: EC, ETT
3362 REAL :: RC, RCS, RCT, RCQ, RCSOIL, FF, WS, SLV, DESDT
3363 REAL :: SIGMA, PC, CMC2MS, SGX, DENOM, RTX, ETT1
3364 INTEGER :: K
3365 REAL, DIMENSION(1:NROOT) :: PART, GX
3367 SLV = 2.501E+6
3368 SIGMA = 5.67E-8
3369 ETT = 0.0
3370 DO K = 1, NSOIL
3371 ET(K) = 0.
3372 END DO
3374 ! ----------------------------------------------------------------------
3375 ! INITIALIZE CANOPY RESISTANCE MULTIPLIER TERMS.
3376 ! ----------------------------------------------------------------------
3377 RCS = 0.0
3378 RCT = 0.0
3379 RCQ = 0.0
3380 RCSOIL = 0.0
3382 ! ----------------------------------------------------------------------
3383 ! CONTRIBUTION DUE TO INCOMING SOLAR RADIATION
3384 ! ----------------------------------------------------------------------
3385 FF = 0.55*2.0* SX*697.7 * 60/ (RGL * LAI)
3386 RCS = (FF + RSMIN / RSMAX) / (1.0+ FF)
3387 RCS = MAX (RCS,0.0001)
3388 ! ----------------------------------------------------------------------
3389 ! CONTRIBUTION DUE TO AIR TEMPERATURE AT FIRST MODEL LEVEL ABOVE GROUND
3390 ! RCT EXPRESSION FROM NOILHAN AND PLANTON (1989, MWR).
3391 ! ----------------------------------------------------------------------
3392 RCT = 1.0- 0.0016* ( (298 - TA)**2.0)
3393 RCT = MAX (RCT,0.0001)
3394 ! ----------------------------------------------------------------------
3395 ! CONTRIBUTION DUE TO VAPOR PRESSURE DEFICIT AT FIRST MODEL LEVEL.
3396 ! RCQ EXPRESSION FROM SSIB (Niyogi and Raman, 1997)
3397 ! ----------------------------------------------------------------------
3398 EA = 6.11*EXP((2.5*10.**6./461.51)*(TA-273.15)/(273.15*TA) )
3399 WS = 0.622*EA/1013
3400 RCQ = 1.0/ (1.0+ HS * (WS - QA))
3401 RCQ = MAX (RCQ,0.01)
3402 ! ----------------------------------------------------------------------
3403 ! CONTRIBUTION DUE TO SOIL MOISTURE AVAILABILITY.
3404 ! DETERMINE CONTRIBUTION FROM EACH SOIL LAYER, THEN ADD THEM UP.
3405 ! ----------------------------------------------------------------------
3406 DO K = 1, NROOT
3407 GX(K) = (SMC(K) - SMCWLT) / (SMCREF - SMCWLT)
3408 IF (GX(K) > 1.) GX(K) = 1.
3409 IF (GX(K) < 0.) GX(K) = 0.
3410 PART (K) = ( -DZVR (K)/ ZSOIL (3)) * GX(K)
3411 END DO
3413 SGX =0.0
3414 DO K = 1, NROOT
3415 SGX = SGX + GX (K)
3416 RCSOIL = RCSOIL + PART (K)
3417 END DO
3418 SGX =SGX / NROOT
3420 RCSOIL = MAX (RCSOIL,0.0001)
3422 RC = RSMIN / (LAI * RCS * RCT * RCQ * RCSOIL)
3423 DESDT = 0.622*SLV*EA/461.51/TA/TA/1013
3424 DELTA = (SLV / CPP)* DESDT
3425 RR = (4.* EPSV *SIGMA * 287.04 / CPP)* (TA **4.)/ (TS * CH) + 1.0
3426 PC = (RR + DELTA)/ (RR * (1. + RC * CH) + DELTA)
3428 IF (CMC .ne. 0.0) THEN
3429 ETT1 = SHDFAC * PC * ETP1 * (1.0- (CMC / CMCMAX) ** CFACTR) * 0.001
3430 ELSE
3431 ETT1 = SHDFAC * PC * ETP1 * 0.001
3432 ENDIF
3434 DENOM = 0.
3435 DO K = 1, NROOT
3436 RTX= (-DZVR (K)/ ZSOIL (3)) + GX(K) - SGX
3437 GX (K) = GX (K) * MAX ( RTX, 0. )
3438 DENOM = DENOM + GX (K)
3439 END DO
3440 IF (DENOM .le. 0.0) DENOM =1.
3442 DO K = 1, NROOT
3443 ET(K) = ETT1 * GX (K) / DENOM
3444 ETT = ETT + ET (K)
3445 END DO
3448 IF (CMC > 0.0) THEN
3449 EC = SHDFAC * ( ( CMC / CMCMAX ) ** CFACTR ) * ETP1 * 0.001
3450 ELSE
3451 EC = 0.0
3452 END IF
3453 CMC2MS = CMC / DELT
3454 EC = MIN ( CMC2MS, EC )
3456 END SUBROUTINE TRANSP
3457 ! ----------------------------------------------------------------------
3458 ! SUBROUTINE SMFLX
3459 ! ----------------------------------------------------------------------
3461 SUBROUTINE SMFLX (SMCP,SMC,NSOIL,CMCP,CMC,DT,PRCP1,ZSOIL, &
3462 & SMCMAX,BEXP,SMCWLT,DKSAT,DWSAT, &
3463 & SHDFAC,CMCMAX,RUNOFF1,RUNOFF2,RUNOFF3, &
3464 EDIR,EC,ET,DRIP)
3466 ! CALCULATE SOIL MOISTURE FLUX. THE SOIL MOISTURE CONTENT IS UPDATED WITH
3467 ! PROGNOSTIC EQNS. THE CANOPY MOISTURE CONTENT (CMC) IS ALSO UPDATED.
3468 ! ----------------------------------------------------------------------
3469 IMPLICIT NONE
3471 INTEGER, INTENT(IN) :: NSOIL
3472 INTEGER :: I,K
3473 REAL, INTENT(IN) :: BEXP, CMCMAX, DKSAT,DWSAT, DT, EC, EDIR, &
3474 PRCP1, SHDFAC, SMCMAX, SMCWLT
3475 REAL, INTENT(OUT) :: DRIP, RUNOFF1, RUNOFF2, RUNOFF3
3476 REAL, INTENT(IN) :: CMCP
3477 REAL, INTENT(OUT) :: CMC
3478 REAL, DIMENSION(1:NSOIL), INTENT(IN) :: ZSOIL, ET
3479 REAL, DIMENSION(1:NSOIL), INTENT(IN) :: SMCP
3480 REAL, DIMENSION(1:NSOIL), INTENT(OUT) :: SMC
3481 REAL, DIMENSION(1:NSOIL) :: AI, BI, CI, STCF,RHSTS, RHSTT
3482 REAL :: EXCESS,PCPDRP,RHSCT,TRHSCT
3485 ! ----------------------------------------------------------------------
3486 ! ADD PRECIPITATION TO EXISTING CMC.IF RESULTING AMT EXCEEDS MAX CAPACITY,
3487 ! IT BECOMES DRIP AND WILL FALL TO THE GRND.
3488 ! ----------------------------------------------------------------------
3489 RHSCT = SHDFAC * PRCP1 * 0.001 /3600. - EC
3490 DRIP = 0.
3491 TRHSCT = DT * RHSCT
3492 EXCESS = CMCP + TRHSCT
3494 ! ----------------------------------------------------------------------
3495 ! PCPDRP IS THE COMBINED PRCP1 AND DRIP (FROM CMCP) THAT GOES INTO THE
3496 ! SOIL
3497 ! ----------------------------------------------------------------------
3498 IF (EXCESS > CMCMAX) DRIP = EXCESS - CMCMAX
3499 PCPDRP = (1. - SHDFAC) * PRCP1 * 0.001 /3600. + DRIP / DT
3501 ! ----------------------------------------------------------------------
3502 ! CALL SUBROUTINES SRT AND SSTEP TO SOLVE THE SOIL MOISTURE
3503 ! TENDENCY EQUATIONS.
3504 ! ----------------------------------------------------------------------
3505 CALL SRT (RHSTT,EDIR,ET,SMCP,NSOIL,PCPDRP,ZSOIL,DWSAT,DKSAT, &
3506 SMCMAX,BEXP,RUNOFF1,RUNOFF2,DT,SMCWLT,AI,BI,CI)
3508 CALL SSTEP (SMCP,SMC,CMCP,CMC,RHSTT,RHSCT,DT,NSOIL,SMCMAX, &
3509 CMCMAX,RUNOFF3,ZSOIL,AI,BI,CI)
3510 ! ----------------------------------------------------------------------
3511 END SUBROUTINE SMFLX
3512 ! ----------------------------------------------------------------------
3514 SUBROUTINE SRT (RHSTT,EDIR,ET,SMCP,NSOIL,PCPDRP,ZSOIL,DWSAT, &
3515 DKSAT,SMCMAX,BEXP,RUNOFF1, &
3516 RUNOFF2,DT,SMCWLT,AI,BI,CI)
3518 ! ----------------------------------------------------------------------
3519 ! CALCULATE THE RIGHT HAND SIDE OF THE TIME TENDENCY TERM OF THE SOIL
3520 ! WATER DIFFUSION EQUATION. ALSO TO COMPUTE ( PREPARE ) THE MATRIX
3521 ! COEFFICIENTS FOR THE TRI-DIAGONAL MATRIX OF THE IMPLICIT TIME SCHEME.
3522 ! ----------------------------------------------------------------------
3523 IMPLICIT NONE
3524 INTEGER, INTENT(IN) :: NSOIL
3525 INTEGER :: K, KS
3527 REAL, INTENT(IN) :: BEXP, DKSAT, DT, DWSAT, EDIR, &
3528 PCPDRP, SMCMAX, SMCWLT
3529 REAL, INTENT(OUT) :: RUNOFF1, RUNOFF2
3530 REAL, DIMENSION(1:NSOIL), INTENT(IN) :: SMCP, ZSOIL, ET
3531 REAL, DIMENSION(1:NSOIL), INTENT(OUT) :: RHSTT
3532 REAL, DIMENSION(1:NSOIL), INTENT(OUT) :: AI, BI, CI
3533 REAL, DIMENSION(1:NSOIL) :: DDMAX
3534 REAL :: DD, DDT, DDZ, DDZ2, DENOM, &
3535 DENOM2, DSMDZ, DSMDZ2, DT1, &
3536 INFMAX,MXSMC,MXSMC2,NUMER,PDDUM, &
3537 PX,SMCAV, SSTT, PAR, &
3538 VAL, WCND, WCND2, WDF, WDF2,KDT
3540 ! ----------------------------------------------------------------------
3541 ! DETERMINE RAINFALL INFILTRATION RATE AND RUNOFF. INCLUDE THE
3542 ! INFILTRATION FORMULE FROM SCHAAKE AND KOREN MODEL.
3543 ! MODIFIED BY Q DUAN
3544 ! ----------------------------------------------------------------------
3546 PDDUM = PCPDRP
3547 RUNOFF1 = 0.0
3548 PAR = 2.0E-6
3550 IF (PCPDRP /= 0.0) THEN
3551 SMCAV = SMCMAX - SMCWLT
3552 DDMAX (1) = - ZSOIL (1)* SMCAV
3553 DDMAX (1) = DDMAX (1)* (1.0- (SMCP (1) - SMCWLT)/ SMCAV)
3554 DDMAX (2) = (ZSOIL (1) - ZSOIL (2))* SMCAV
3555 DDMAX (2) = DDMAX (2)* (1.0- (SMCP (2) - SMCWLT)/ SMCAV)
3556 DDMAX (3) = (ZSOIL (2) - ZSOIL (3))* SMCAV
3557 DDMAX (3) = DDMAX (3)* (1.0- (SMCP (3) - SMCWLT)/ SMCAV)
3559 DD = DDMAX(1)+DDMAX(2)+DDMAX(3)
3560 DT1 = DT/86400
3561 KDT = 3.0 * DKSAT / PAR
3562 VAL = (1. - EXP ( - KDT * DT1))
3563 DDT = DD * VAL
3564 PX = PCPDRP * DT
3565 IF (PX < 0.0) PX = 0.0
3567 INFMAX = (PX * (DDT / (PX + DDT)))/ DT
3568 MXSMC = SMCP (1)
3569 CALL WDFCND (WDF,WCND,MXSMC,SMCMAX,BEXP,DKSAT,DWSAT)
3570 INFMAX = MAX (INFMAX,WCND)
3571 INFMAX = MIN (INFMAX,PX/DT)
3574 IF (PCPDRP > INFMAX) THEN
3575 RUNOFF1 = PCPDRP - INFMAX
3576 PDDUM = INFMAX
3577 END IF
3578 END IF
3579 ! ----------------------------------------------------------------------
3580 ! TOP LAYER
3581 ! ----------------------------------------------------------------------
3582 CALL WDFCND (WDF,WCND,SMCP(1),SMCMAX,BEXP,DKSAT,DWSAT)
3583 DDZ = 1. / ( - .5 * ZSOIL (2) )
3584 AI (1) = 0.0
3585 BI (1) = WDF * DDZ / ( - ZSOIL (1) )
3586 CI (1) = - BI (1)
3587 DSMDZ = (SMCP (1) - SMCP (2) )/( - 0.5 * ZSOIL(2))
3588 RHSTT (1) = (WDF * DSMDZ + WCND- PDDUM + EDIR + ET(1))/ ZSOIL (1)
3589 SSTT = WDF * DSMDZ + WCND+ EDIR + ET(1)
3591 ! ----------------------------------------------------------------------
3592 ! LOOP THRU THE REMAINING SOIL LAYERS, REPEATING THE ABV PROCESS
3593 ! ----------------------------------------------------------------------
3594 DDZ2 = 0.0
3595 DO K = 2,NSOIL-1
3596 DENOM2 = (ZSOIL (K -1) - ZSOIL (K))
3597 IF (K /= NSOIL-1) THEN
3598 MXSMC2 = SMCP (K)
3599 CALL WDFCND (WDF2,WCND2,MXSMC2,SMCMAX,BEXP,DKSAT,DWSAT)
3600 DENOM = (ZSOIL (K -1) - ZSOIL (K +1))
3601 DSMDZ2 = (SMCP (K) - SMCP (K +1)) / (DENOM * 0.5)
3602 DDZ2 = 2.0 / DENOM
3603 CI (K) = - WDF2 * DDZ2 / DENOM2
3604 ELSE
3605 CALL WDFCND (WDF2,WCND2,SMCP(NSOIL-1),SMCMAX,BEXP,DKSAT,DWSAT)
3606 DSMDZ2 = 0.0
3607 CI (K) = 0.0
3608 END IF
3609 NUMER = (WDF2 * DSMDZ2) - (WDF * DSMDZ) &
3610 - WCND+ ET(K)
3611 RHSTT (K) = NUMER / ( - DENOM2)
3612 AI (K) = - WDF * DDZ / DENOM2
3613 BI (K) = - ( AI (K) + CI (K) )
3614 IF (K .eq. NSOIL-1) THEN
3615 RUNOFF2 = 0.0
3616 END IF
3617 IF (K .ne. NSOIL-1) THEN
3618 WDF = WDF2
3619 WCND = WCND2
3620 DSMDZ = DSMDZ2
3621 DDZ = DDZ2
3622 END IF
3623 END DO
3624 ! ----------------------------------------------------------------------
3625 END SUBROUTINE SRT
3626 ! ----------------------------------------------------------------------
3628 SUBROUTINE SSTEP (SMCP,SMC,CMCP,CMC,RHSTT,RHSCT,DT, &
3629 NSOIL,SMCMAX,CMCMAX,RUNOFF3,ZSOIL, &
3630 AI,BI,CI)
3632 ! ----------------------------------------------------------------------
3633 ! SUBROUTINE SSTEP
3634 ! ----------------------------------------------------------------------
3635 ! CALCULATE/UPDATE SOIL MOISTURE CONTENT VALUES AND CANOPY MOISTURE
3636 ! CONTENT VALUES.
3637 ! ----------------------------------------------------------------------
3638 IMPLICIT NONE
3639 INTEGER, INTENT(IN) :: NSOIL
3640 INTEGER :: I, K, KK11
3642 REAL, INTENT(IN) :: CMCMAX, DT, SMCMAX
3643 REAL, INTENT(OUT) :: RUNOFF3
3644 REAL, INTENT(IN) :: CMCP
3645 REAL, INTENT(OUT) :: CMC
3646 REAL, DIMENSION(1:NSOIL), INTENT(IN) :: SMCP, ZSOIL
3647 REAL, DIMENSION(1:NSOIL), INTENT(OUT) :: SMC
3648 REAL, DIMENSION(1:NSOIL), INTENT(INOUT) :: RHSTT
3649 REAL, DIMENSION(1:NSOIL), INTENT(INOUT) :: AI, BI, CI
3650 REAL, DIMENSION(1:NSOIL) :: RHSTTin, SMCOUT,SMCIN
3651 REAL, DIMENSION(1:NSOIL) :: CIin
3652 REAL :: DDZ, RHSCT, WPLUS, STOT
3654 ! ----------------------------------------------------------------------
3655 ! CREATE 'AMOUNT' VALUES OF VARIABLES TO BE INPUT TO THE
3656 ! TRI-DIAGONAL MATRIX ROUTINE.
3657 ! ----------------------------------------------------------------------
3658 DO K = 1,NSOIL-1
3659 RHSTT (K) = RHSTT (K) * DT
3660 AI (K) = AI (K) * DT
3661 BI (K) = 1. + BI (K) * DT
3662 CI (K) = CI (K) * DT
3663 END DO
3664 ! ----------------------------------------------------------------------
3665 ! COPY VALUES FOR INPUT VARIABLES BEFORE CALL TO ROSR12
3666 ! ----------------------------------------------------------------------
3667 DO K = 1,NSOIL-1
3668 RHSTTin (K) = RHSTT (K)
3669 END DO
3670 DO K = 1,NSOIL-1
3671 CIin (K) = CI (K)
3672 END DO
3673 ! ----------------------------------------------------------------------
3674 ! CALL ROSR12 TO SOLVE THE TRI-DIAGONAL MATRIX
3675 ! ----------------------------------------------------------------------
3676 CALL ROSR12 (CI,AI,BI,CIin,RHSTTin,RHSTT,NSOIL-1)
3677 ! ----------------------------------------------------------------------
3678 ! SUM THE PREVIOUS SMC VALUE AND THE MATRIX SOLUTION TO GET A
3679 ! NEW VALUE. MIN ALLOWABLE VALUE OF SMC WILL BE 0.02.
3680 ! RUNOFF3: RUNOFF WITHIN SOIL LAYERS
3681 ! ----------------------------------------------------------------------
3682 WPLUS = 0.0
3683 RUNOFF3 = 0.
3685 DDZ = - ZSOIL (1)
3686 DO K = 1,NSOIL-1
3687 IF (K /= 1) DDZ = ZSOIL (K - 1) - ZSOIL (K)
3688 SMCOUT (K) = SMCP (K) + CI (K) + WPLUS / DDZ
3689 STOT = SMCOUT (K)
3690 IF (STOT > SMCMAX) THEN
3691 IF (K .eq. 1) THEN
3692 DDZ = - ZSOIL (1)
3693 ELSE
3694 KK11 = K - 1
3695 DDZ = - ZSOIL (K) + ZSOIL (KK11)
3696 END IF
3697 WPLUS = (STOT - SMCMAX) * DDZ
3698 ELSE
3699 WPLUS = 0.
3700 END IF
3701 SMC (K) = MAX ( MIN (STOT,SMCMAX),0.066 )
3702 END DO
3704 ! ----------------------------------------------------------------------
3705 ! UPDATE CANOPY WATER CONTENT/INTERCEPTION (CMC). CONVERT RHSCT TO
3706 ! AN 'AMOUNT' VALUE AND ADD TO PREVIOUS CMC VALUE TO GET NEW CMC.
3707 ! ----------------------------------------------------------------------
3708 RUNOFF3 = WPLUS
3709 CMC = CMCP + DT * RHSCT
3710 IF (CMC < 1.E-20) CMC = 0.0
3711 CMC = MIN (CMC,CMCMAX)
3713 ! ----------------------------------------------------------------------
3714 END SUBROUTINE SSTEP
3715 ! ----------------------------------------------------------------------
3717 SUBROUTINE WDFCND (WDF,WCND,SMC,SMCMAX,BEXP,DKSAT,DWSAT)
3719 ! ----------------------------------------------------------------------
3720 ! SUBROUTINE WDFCND
3721 ! ----------------------------------------------------------------------
3722 ! CALCULATE SOIL WATER DIFFUSIVITY AND SOIL HYDRAULIC CONDUCTIVITY.
3723 ! ----------------------------------------------------------------------
3724 IMPLICIT NONE
3725 REAL BEXP
3726 REAL DKSAT
3727 REAL DWSAT
3728 REAL EXPON
3729 REAL FACTR1
3730 REAL FACTR2
3731 REAL SMC
3732 REAL SMCMAX
3733 REAL WCND
3735 ! ----------------------------------------------------------------------
3736 ! CALC THE RATIO OF THE ACTUAL TO THE MAX PSBL SOIL H2O CONTENT
3737 ! ----------------------------------------------------------------------
3738 REAL WDF
3739 FACTR1 = 0.05 / SMCMAX
3741 ! ----------------------------------------------------------------------
3742 ! PREP AN EXPNTL COEF AND CALC THE SOIL WATER DIFFUSIVITY AND CONDUCTIVITY
3743 ! ----------------------------------------------------------------------
3744 FACTR2 = SMC / SMCMAX
3745 FACTR1 = MIN(FACTR1,FACTR2)
3746 EXPON = BEXP + 2.0
3747 WDF = DWSAT * FACTR2 ** EXPON
3748 EXPON = (2.0 * BEXP) + 3.0
3749 WCND = DKSAT * FACTR2 ** EXPON
3751 ! ----------------------------------------------------------------------
3752 END SUBROUTINE WDFCND
3753 ! ----------------------------------------------------------------------
3754 ! SUBROUTINE ROSR12
3755 ! ----------------------------------------------------------------------
3756 ! INVERT (SOLVE) THE TRI-DIAGONAL MATRIX PROBLEM SHOWN BELOW:
3757 ! ### ### ### ### ### ###
3758 ! #B(1), C(1), 0 , 0 , 0 , . . . , 0 # # # # #
3759 ! #A(2), B(2), C(2), 0 , 0 , . . . , 0 # # # # #
3760 ! # 0 , A(3), B(3), C(3), 0 , . . . , 0 # # # # D(3) #
3761 ! # 0 , 0 , A(4), B(4), C(4), . . . , 0 # # P(4) # # D(4) #
3762 ! # 0 , 0 , 0 , A(5), B(5), . . . , 0 # # P(5) # # D(5) #
3763 ! # . . # # . # = # . #
3764 ! # . . # # . # # . #
3765 ! # . . # # . # # . #
3766 ! # 0 , . . . , 0 , A(M-2), B(M-2), C(M-2), 0 # #P(M-2)# #D(M-2)#
3767 ! # 0 , . . . , 0 , 0 , A(M-1), B(M-1), C(M-1)# #P(M-1)# #D(M-1)#
3768 ! # 0 , . . . , 0 , 0 , 0 , A(M) , B(M) # # P(M) # # D(M) #
3769 ! ### ### ### ### ### ###
3770 ! ----------------------------------------------------------------------
3772 SUBROUTINE ROSR12 (P,A,B,C,D,DELTA,NSOIL)
3773 IMPLICIT NONE
3775 INTEGER, INTENT(IN) :: NSOIL
3776 INTEGER :: K, KK
3778 REAL, DIMENSION(1:NSOIL), INTENT(IN):: A, B, D
3779 REAL, DIMENSION(1:NSOIL),INTENT(INOUT):: C,P,DELTA
3781 ! ----------------------------------------------------------------------
3782 ! INITIALIZE EQN COEF C FOR THE LOWEST SOIL LAYER
3783 ! ----------------------------------------------------------------------
3784 C (NSOIL) = 0.0
3785 P (1) = - C (1) / B (1)
3786 DELTA (1) = D (1) / B (1)
3787 DO K = 2,NSOIL
3788 P (K) = - C (K) * ( 1.0 / (B (K) + A (K) * P (K -1)) )
3789 DELTA (K) = (D (K) - A (K)* DELTA (K -1))* (1.0/ (B (K) + A (K)&
3790 * P (K -1)))
3791 END DO
3792 ! ----------------------------------------------------------------------
3793 ! SET P TO DELTA FOR LOWEST SOIL LAYER
3794 ! ----------------------------------------------------------------------
3795 P (NSOIL) = DELTA (NSOIL)
3797 ! ----------------------------------------------------------------------
3798 ! ADJUST P FOR SOIL LAYERS 2 THRU NSOIL
3799 ! ----------------------------------------------------------------------
3800 DO K = 2,NSOIL
3801 KK = NSOIL - K + 1
3802 P (KK) = P (KK) * P (KK +1) + DELTA (KK)
3803 END DO
3804 ! ----------------------------------------------------------------------
3805 END SUBROUTINE ROSR12
3806 !----------------------------------------------------------------------
3808 SUBROUTINE SHFLX (SSOIL,STC,SMC,SMCMAX,NSOIL,T1,DT,YY,ZZ1,ZSOIL, &
3809 TBOT,ZBOT,SMCWLT,DF1,QUARTZ,CSOIL,CAPR)
3811 ! ----------------------------------------------------------------------
3812 ! SUBROUTINE SHFLX
3813 ! ----------------------------------------------------------------------
3814 ! UPDATE THE TEMPERATURE STATE OF THE SOIL COLUMN BASED ON THE THERMAL
3815 ! DIFFUSION EQUATION AND UPDATE THE FROZEN SOIL MOISTURE CONTENT BASED
3816 ! ON THE TEMPERATURE.
3817 ! ----------------------------------------------------------------------
3818 IMPLICIT NONE
3820 INTEGER, INTENT(IN) :: NSOIL
3821 INTEGER :: I
3823 REAL, INTENT(IN) :: DF1,DT,SMCMAX, SMCWLT, TBOT,YY, ZBOT,ZZ1, QUARTZ
3824 REAL, INTENT(IN) :: CSOIL, CAPR
3825 REAL, INTENT(INOUT) :: T1
3826 REAL, INTENT(OUT) :: SSOIL
3827 REAL, DIMENSION(1:NSOIL), INTENT(IN) :: SMC,ZSOIL
3828 REAL, DIMENSION(1:NSOIL), INTENT(INOUT) :: STC
3829 REAL, DIMENSION(1:NSOIL) :: AI, BI, CI, STCF,RHSTS
3831 ! ----------------------------------------------------------------------
3832 ! HRT ROUTINE CALCS THE RIGHT HAND SIDE OF THE SOIL TEMP DIF EQN
3833 ! ----------------------------------------------------------------------
3835 ! Land case
3837 CALL HRT (RHSTS,STC,SMC,SMCMAX,NSOIL,ZSOIL,YY,ZZ1,TBOT, &
3838 ZBOT,DT,DF1,AI,BI,CI,QUARTZ,CSOIL,CAPR)
3840 CALL HSTEP (STCF,STC,RHSTS,DT,NSOIL,AI,BI,CI)
3842 DO I = 1,NSOIL
3843 STC (I) = STCF (I)
3844 ENDDO
3846 ! ----------------------------------------------------------------------
3847 ! CALCULATE SURFACE SOIL HEAT FLUX
3848 ! ----------------------------------------------------------------------
3849 T1 = (YY + (ZZ1- 1.0) * STC (1)) / ZZ1
3850 SSOIL = DF1 * (STC (1) - T1) / (0.5 * ZSOIL (1))
3852 ! ----------------------------------------------------------------------
3853 END SUBROUTINE SHFLX
3854 ! ----------------------------------------------------------------------
3855 ! SUBROUTINE HRT
3856 ! ----------------------------------------------------------------------
3857 ! CALCULATE THE RIGHT HAND SIDE OF THE TIME TENDENCY TERM OF THE SOIL
3858 ! THERMAL DIFFUSION EQUATION. ALSO TO COMPUTE ( PREPARE ) THE MATRIX
3859 ! COEFFICIENTS FOR THE TRI-DIAGONAL MATRIX OF THE IMPLICIT TIME SCHEME.
3860 ! ----------------------------------------------------------------------
3862 SUBROUTINE HRT (RHSTS,STC,SMC,SMCMAX,NSOIL,ZSOIL,YY,ZZ1, &
3863 TBOT,ZBOT,DT,DF1,AI,BI,CI,QUARTZ,CSOIL,CAPR)
3865 IMPLICIT NONE
3866 LOGICAL :: ITAVG
3867 INTEGER, INTENT(IN) :: NSOIL
3868 INTEGER :: I, K
3870 REAL, INTENT(IN) :: DF1, DT,SMCMAX ,TBOT,YY,ZZ1, ZBOT, QUARTZ, CSOIL, CAPR
3871 REAL, DIMENSION(1:NSOIL), INTENT(IN) :: SMC,STC,ZSOIL
3872 REAL, DIMENSION(1:NSOIL), INTENT(OUT) :: RHSTS
3873 REAL, DIMENSION(1:NSOIL), INTENT(OUT) :: AI, BI,CI
3874 REAL :: DDZ, DDZ2, DENOM, DF1K, DTSDZ,DF1N, &
3875 DTSDZ2,HCPCT,QTOT,SSOIL,SICE,TAVG,TBK, &
3876 TBK1,TSNSR,TSURF
3877 REAL, PARAMETER :: CAIR = 1004.0, CH2O = 4.2E6
3880 ! ----------------------------------------------------------------------
3881 ! INITIALIZE LOGICAL FOR SOIL LAYER TEMPERATURE AVERAGING.
3882 ! ----------------------------------------------------------------------
3883 ITAVG = .TRUE.
3885 ! ----------------------------------------------------------------------
3886 ! TOP SOIL LAYER
3887 ! ----------------------------------------------------------------------
3888 HCPCT = SMC (1)* CH2O + (1.0- SMCMAX)* CSOIL + (SMCMAX - SMC (1))&
3889 * CAIR
3890 DDZ = 1.0 / ( -0.5 * ZSOIL (2) )
3891 AI (1) = 0.0
3892 CI (1) = (DF1 * DDZ) / (ZSOIL (1) * HCPCT)
3894 ! ----------------------------------------------------------------------
3895 ! CALCULATE THE VERTICAL SOIL TEMP GRADIENT BTWN THE 1ST AND 2ND SOIL
3896 ! LAYERS. THEN CALCULATE THE SUBSURFACE HEAT FLUX.
3897 ! ----------------------------------------------------------------------
3898 BI (1) = - CI (1) + DF1 / (0.5 * ZSOIL (1) * ZSOIL (1)* HCPCT * &
3899 ZZ1)
3900 DTSDZ = (STC (1) - STC (2)) / (-0.5 * ZSOIL (2))
3901 SSOIL = DF1 * (STC (1) - YY) / (0.5 * ZSOIL (1) * ZZ1)
3902 DENOM = (ZSOIL (1) * HCPCT)
3904 ! ----------------------------------------------------------------------
3905 ! NEXT CAPTURE THE VERTICAL DIFFERENCE OF THE HEAT FLUX AT TOP AND
3906 ! BOTTOM OF FIRST SOIL LAYER FOR USE IN HEAT FLUX CONSTRAINT
3907 ! ----------------------------------------------------------------------
3908 RHSTS (1) = (DF1 * DTSDZ - SSOIL) / DENOM
3909 QTOT = -1.0* RHSTS (1)* DENOM
3910 IF (ITAVG) THEN
3911 TSURF = (YY + (ZZ1-1) * STC (1)) / ZZ1
3912 CALL TBND (STC (1),STC (2),ZSOIL,ZBOT,1,NSOIL,TBK)
3913 ENDIF
3914 DDZ2 = 0.0
3915 DF1N = DF1
3917 ! ----------------------------------------------------------------------
3918 ! LOOP THRU THE REMAINING SOIL LAYERS, REPEATING THE ABOVE PROCESS
3919 ! (EXCEPT SUBSFC OR "GROUND" HEAT FLUX NOT REPEATED IN LOWER LAYERS)
3920 ! ----------------------------------------------------------------------
3921 DO K = 2,NSOIL
3922 ! ----------------------------------------------------------------------
3923 ! THIS SECTION FOR LAYER 2 OR GREATER, BUT NOT LAST LAYER.
3924 ! ----------------------------------------------------------------------
3925 IF (K < NSOIL-1 ) THEN
3926 HCPCT = SMC (K)* CH2O + (1.0- SMCMAX)* CSOIL + (SMCMAX - SMC ( &
3927 K))* CAIR
3928 CALL TDFCND (DF1K, SMC(K), QUARTZ, SMCMAX)
3929 DENOM = 0.5 * ( ZSOIL (K -1) - ZSOIL (K +1) )
3930 DTSDZ2 = (STC (K) - STC (K +1) ) / DENOM
3931 DDZ2 = 2. / (ZSOIL (K -1) - ZSOIL (K +1))
3933 ! ----------------------------------------------------------------------
3934 ! IF TEMPERATURE AVERAGING INVOKED (ITAVG=TRUE; ELSE SKIP): CALCULATE
3935 ! TEMP AT BOTTOM OF LAYER.
3936 ! ----------------------------------------------------------------------
3937 CI (K) = - DF1K * DDZ2 / ( (ZSOIL (K -1) - ZSOIL (K)) * &
3938 HCPCT)
3939 IF (ITAVG) THEN
3940 CALL TBND (STC (K),STC (K +1),ZSOIL,ZBOT,K,NSOIL,TBK1)
3941 END IF
3943 ELSEIF (K == NSOIL-1) THEN
3945 HCPCT = SMC (K)* CH2O + (1.0- SMCMAX)* CSOIL + (SMCMAX- SMC ( &
3946 K))* CAIR
3947 CALL TDFCND (DF1K, SMC(K), QUARTZ, SMCMAX)
3948 DENOM = 0.5 * ( ZSOIL (K -1) - ZSOIL (K +1) )
3949 DTSDZ2 = (STC (K) - STC (K +1) ) / DENOM
3950 DDZ2 = 2. / (ZSOIL (K -1) - ZSOIL (K +1))
3951 !-----------------------------------------------------------------------
3952 ! IF TEMPERATURE AVERAGING INVOKED (ITAVG=TRUE; ELSE SKIP): CALCULATE
3953 ! TEMP AT BOTTOM OF LAST LAYER.
3954 ! ----------------------------------------------------------------------
3955 CI (K) = - DF1K * DDZ2 / ( (ZSOIL (K -1) - ZSOIL (K)) * &
3956 HCPCT)
3957 IF (ITAVG) THEN
3958 CALL TBND (STC (K),TBOT,ZSOIL,ZBOT,K,NSOIL,TBK1)
3959 END IF
3960 ELSE
3961 ! ----------------------------------------------------------------------
3962 ! SPECIAL CASE OF BOTTOM LAYER (CONCRETE ROOF)
3963 ! ----------------------------------------------------------------------
3964 HCPCT = CAPR * 4.1868 * 1.E6
3965 DF1K = 3.24
3966 ! ----------------------------------------------------------------------
3967 ! CALC THE VERTICAL TEMP GRADIENT THRU BOTTOM LAYER.
3968 ! ----------------------------------------------------------------------
3969 DENOM = .5 * (ZSOIL (K -1) + ZSOIL (K)) - ZBOT
3970 DTSDZ2 = (STC (K) - TBOT) / DENOM
3971 ! ----------------------------------------------------------------------
3972 ! IF TEMPERATURE AVERAGING INVOKED (ITAVG=TRUE; ELSE SKIP): CALCULATE
3973 ! TEMP AT BOTTOM OF LAST LAYER.
3974 ! ----------------------------------------------------------------------
3975 CI (K) = 0.
3976 IF (ITAVG) THEN
3977 CALL TBND (STC (K),TBOT,ZSOIL,ZBOT,K,NSOIL,TBK1)
3978 END IF
3979 ! ----------------------------------------------------------------------
3980 ! THIS ENDS SPECIAL LOOP FOR BOTTOM LAYER.
3981 END IF
3982 ! ----------------------------------------------------------------------
3983 ! CALCULATE RHSTS FOR THIS LAYER AFTER CALC'NG A PARTIAL PRODUCT.
3984 ! ----------------------------------------------------------------------
3985 DENOM = ( ZSOIL (K) - ZSOIL (K -1) ) * HCPCT
3986 RHSTS (K) = ( DF1K * DTSDZ2- DF1N * DTSDZ ) / DENOM
3987 QTOT = -1.0* DENOM * RHSTS (K)
3989 ! ----------------------------------------------------------------------
3990 ! CALC MATRIX COEFS, AI, AND BI FOR THIS LAYER.
3991 ! ----------------------------------------------------------------------
3992 AI (K) = - DF1N * DDZ / ( (ZSOIL (K -1) - ZSOIL (K)) * HCPCT)
3994 ! ----------------------------------------------------------------------
3995 ! RESET VALUES OF DF1, DTSDZ, DDZ, AND TBK FOR LOOP TO NEXT SOIL LAYER.
3996 ! ----------------------------------------------------------------------
3997 BI (K) = - (AI (K) + CI (K))
3998 TBK = TBK1
3999 DF1N = DF1K
4000 DTSDZ = DTSDZ2
4001 DDZ = DDZ2
4002 END DO
4003 ! ----------------------------------------------------------------------
4004 END SUBROUTINE HRT
4005 ! ----------------------------------------------------------------------
4007 SUBROUTINE HSTEP (STCOUT,STCIN,RHSTS,DT,NSOIL,AI,BI,CI)
4008 ! CALCULATE/UPDATE THE SOIL TEMPERATURE FIELD.
4009 ! ----------------------------------------------------------------------
4010 IMPLICIT NONE
4011 INTEGER, INTENT(IN) :: NSOIL
4012 INTEGER :: K
4014 REAL, DIMENSION(1:NSOIL), INTENT(IN):: STCIN
4015 REAL, DIMENSION(1:NSOIL), INTENT(OUT):: STCOUT
4016 REAL, DIMENSION(1:NSOIL), INTENT(INOUT):: RHSTS
4017 REAL, DIMENSION(1:NSOIL), INTENT(INOUT):: AI,BI,CI
4018 REAL, DIMENSION(1:NSOIL) :: RHSTSin
4019 REAL, DIMENSION(1:NSOIL) :: CIin
4020 REAL :: DT
4022 ! ----------------------------------------------------------------------
4023 ! CREATE FINITE DIFFERENCE VALUES FOR USE IN ROSR12 ROUTINE
4024 ! ----------------------------------------------------------------------
4025 DO K = 1,NSOIL
4026 RHSTS (K) = RHSTS (K) * DT
4027 AI (K) = AI (K) * DT
4028 BI (K) = 1. + BI (K) * DT
4029 CI (K) = CI (K) * DT
4030 END DO
4031 ! ----------------------------------------------------------------------
4032 ! COPY VALUES FOR INPUT VARIABLES BEFORE CALL TO ROSR12
4033 ! ----------------------------------------------------------------------
4034 DO K = 1,NSOIL
4035 RHSTSin (K) = RHSTS (K)
4036 END DO
4037 DO K = 1,NSOIL
4038 CIin (K) = CI (K)
4039 END DO
4040 ! ----------------------------------------------------------------------
4041 ! SOLVE THE TRI-DIAGONAL MATRIX EQUATION
4042 ! ----------------------------------------------------------------------
4043 CALL ROSR12 (CI,AI,BI,CIin,RHSTSin,RHSTS,NSOIL)
4044 ! ----------------------------------------------------------------------
4045 ! CALC/UPDATE THE SOIL TEMPS USING MATRIX SOLUTION
4046 ! ----------------------------------------------------------------------
4047 DO K = 1,NSOIL
4048 STCOUT (K) = STCIN (K) + CI (K)
4049 END DO
4050 ! ----------------------------------------------------------------------
4051 END SUBROUTINE HSTEP
4052 ! ----------------------------------------------------------------------
4053 ! ----------------------------------------------------------------------
4055 SUBROUTINE TBND (TU,TB,ZSOIL,ZBOT,K,NSOIL,TBND1)
4057 ! ----------------------------------------------------------------------
4058 ! SUBROUTINE TBND
4059 ! ----------------------------------------------------------------------
4060 ! CALCULATE TEMPERATURE ON THE BOUNDARY OF THE LAYER BY INTERPOLATION OF
4061 ! THE MIDDLE LAYER TEMPERATURES
4062 ! ----------------------------------------------------------------------
4063 IMPLICIT NONE
4064 INTEGER, INTENT(IN) :: NSOIL
4065 INTEGER :: K
4066 REAL, INTENT(IN) :: TB, TU, ZBOT
4067 REAL, INTENT(OUT) :: TBND1
4068 REAL, DIMENSION(1:NSOIL), INTENT(IN) :: ZSOIL
4069 REAL :: ZB, ZUP
4071 ! ----------------------------------------------------------------------
4072 ! USE SURFACE TEMPERATURE ON THE TOP OF THE FIRST LAYER
4073 ! ----------------------------------------------------------------------
4074 IF (K == 1) THEN
4075 ZUP = 0.
4076 ELSE
4077 ZUP = ZSOIL (K -1)
4078 END IF
4079 ! ----------------------------------------------------------------------
4080 ! USE DEPTH OF THE CONSTANT BOTTOM TEMPERATURE WHEN INTERPOLATE
4081 ! TEMPERATURE INTO THE LAST LAYER BOUNDARY
4082 ! ----------------------------------------------------------------------
4083 IF (K == NSOIL) THEN
4084 ZB = 2.* ZBOT - ZSOIL (K)
4085 ELSE
4086 ZB = ZSOIL (K +1)
4087 END IF
4088 ! ----------------------------------------------------------------------
4089 ! LINEAR INTERPOLATION BETWEEN THE AVERAGE LAYER TEMPERATURES
4090 ! ----------------------------------------------------------------------
4092 TBND1 = TU + (TB - TU)* (ZUP - ZSOIL (K))/ (ZUP - ZB)
4093 ! ----------------------------------------------------------------------
4094 END SUBROUTINE TBND
4095 ! ----------------------------------------------------------------------
4096 SUBROUTINE TDFCND (DF, SMC, QZ, SMCMAX)
4097 ! ----------------------------------------------------------------------
4098 ! CALCULATE THERMAL CONDUCTIVITY OF THE SOIL
4099 ! ----------------------------------------------------------------------
4100 ! PETERS-LIDARD APPROACH (PETERS-LIDARD et al., 1998)
4101 ! ----------------------------------------------------------------------
4102 IMPLICIT NONE
4103 REAL, INTENT(IN) :: QZ, SMC, SMCMAX
4104 REAL, INTENT(OUT) :: DF
4105 REAL :: AKE, GAMMD, THKDRY, THKO, &
4106 THKQTZ,THKSAT,THKS,THKW,SATRATIO
4108 ! ----------------------------------------------------------------------
4109 ! IF THE SOIL HAS ANY MOISTURE CONTENT COMPUTE A PARTIAL SUM/PRODUCT
4110 ! OTHERWISE USE A CONSTANT VALUE WHICH WORKS WELL WITH MOST SOILS
4111 ! ----------------------------------------------------------------------
4112 ! THKW ......WATER THERMAL CONDUCTIVITY
4113 ! THKQTZ ....THERMAL CONDUCTIVITY FOR QUARTZ
4114 ! THKO ......THERMAL CONDUCTIVITY FOR OTHER SOIL COMPONENTS
4115 ! THKS ......THERMAL CONDUCTIVITY FOR THE SOLIDS COMBINED(QUARTZ+OTHER)
4116 ! SMCMAX ....POROSITY (= SMCMAX)
4117 ! QZ .........QUARTZ CONTENT (SOIL TYPE DEPENDENT)
4118 ! ----------------------------------------------------------------------
4119 ! USE AS IN PETERS-LIDARD, 1998 (MODIF. FROM JOHANSEN, 1975).
4121 ! PABLO GRUNMANN, 08/17/98
4122 ! REFS.:
4123 ! FAROUKI, O.T.,1986: THERMAL PROPERTIES OF SOILS. SERIES ON ROCK
4124 ! AND SOIL MECHANICS, VOL. 11, TRANS TECH, 136 PP.
4125 ! JOHANSEN, O., 1975: THERMAL CONDUCTIVITY OF SOILS. PH.D. THESIS,
4126 ! UNIVERSITY OF TRONDHEIM,
4127 ! PETERS-LIDARD, C. D., ET AL., 1998: THE EFFECT OF SOIL THERMAL
4128 ! CONDUCTIVITY PARAMETERIZATION ON SURFACE ENERGY FLUXES
4129 ! AND TEMPERATURES. JOURNAL OF THE ATMOSPHERIC SCIENCES,
4130 ! VOL. 55, PP. 1209-1224.
4131 ! ----------------------------------------------------------------------
4132 ! NEEDS PARAMETERS
4133 ! POROSITY(SOIL TYPE):
4134 ! POROS = SMCMAX
4135 ! SATURATION RATIO:
4136 ! PARAMETERS W/(M.K)
4137 SATRATIO = SMC / SMCMAX
4138 ! WATER CONDUCTIVITY:
4139 THKW = 0.57
4140 ! THERMAL CONDUCTIVITY OF "OTHER" SOIL COMPONENTS
4141 ! IF (QZ .LE. 0.2) THKO = 3.0
4142 THKO = 2.0
4143 ! QUARTZ' CONDUCTIVITY
4144 THKQTZ = 7.7
4145 ! SOLIDS' CONDUCTIVITY
4146 THKS = (THKQTZ ** QZ)* (THKO ** (1. - QZ))
4148 ! SATURATED THERMAL CONDUCTIVITY
4149 THKSAT = THKS ** (1. - SMCMAX)* THKW ** (SMCMAX)
4151 ! DRY DENSITY IN KG/M3
4152 GAMMD = (1. - SMCMAX)*2700.
4154 ! DRY THERMAL CONDUCTIVITY IN W.M-1.K-1
4155 THKDRY = (0.135* GAMMD+ 64.7)/ (2700. - 0.947* GAMMD)
4157 ! KERSTEN NUMBER (USING "FINE" FORMULA, VALID FOR SOILS CONTAINING AT
4158 ! LEAST 5% OF PARTICLES WITH DIAMETER LESS THAN 2.E-6 METERS.)
4159 ! (FOR "COARSE" FORMULA, SEE PETERS-LIDARD ET AL., 1998).
4161 IF ( SATRATIO > 0.1 ) THEN
4163 AKE = LOG10 (SATRATIO) + 1.0
4165 ! USE K = KDRY
4166 ELSE
4168 AKE = 0.0
4169 END IF
4170 ! THERMAL CONDUCTIVITY
4172 DF = AKE * (THKSAT - THKDRY) + THKDRY
4173 ! ----------------------------------------------------------------------
4174 END SUBROUTINE TDFCND
4175 ! ----------------------------------------------------------------------
4177 FUNCTION kanda_kawai_svf(lp, lf) RESULT (svf)
4178 IMPLICIT NONE
4179 real, intent(in) :: lp, lf
4180 real :: hovl, vloc, vmod, svf
4182 hovl = lf * lp ** (-0.5) / (1. - lp ** 0.5)
4183 vloc = cos(atan(2. * hovl)) * (2. - 4. / piconst * atan(cos(atan(2. * hovl))))
4184 vmod = 0.1120 * lp * vloc - 0.4817 * lp + 0.0246 * vloc + 0.9570
4185 svf = vloc * vmod
4186 END FUNCTION kanda_kawai_svf
4188 !===========================================================================
4189 END MODULE module_sf_urban