SUEWS API Site
Documentation of SUEWS source code
suews_ctrl_translate.f95
Go to the documentation of this file.
1 !SUEWS_Translate
2 !Translates - new input arrays (v2014b) to existing model variables
3 ! - between arrays for different grids and the model variables
4 !Made by HW&LJ Oct 2014
5 !-----------------------------------------------------------------------------------
6 ! MH 21 Jun 2017 : Added anthropogenic CO2 charasteristic
7 ! MH 16 Jun 2017 : Added biogenic CO2 characteristic
8 ! HCW 13 Dec 2016 : LAIPower and LAIType for all vegetation types now used (previously only DecTr were used)
9 ! HCW 12 Dec 2016 : Switched sign of lng so that input should be -ve for W, +ve for E, as is conventional
10 !Last modified HCW 26 Aug 2016
11 ! NumCapita now uses average of day and night pop density, unless only one is specified
12 !Last modified HCW 06 Jul 2016
13 ! Checks on ESTM fractions
14 ! - default setting to first ESTM Class code if surface not present and ESTM fractions do not sum to 1.
15 !Last modified HCW 29 Jun 2016
16 ! Removed SoilMoistDay and StateDay
17 !Last modified: HCW 16 Jun 2016
18 ! ESTM development for 7 surface types + snow, allowing 3x Paved classes and 5x Bldgs classes
19 ! Currently surface characteristics are averaged here; probably want to average QS instead.
20 !Last modified: TS 13 Apr 2016
21 ! Added AnOHM required variables.
22 !Last modified: LJ 06 Jul 2015
23 ! Changed to read SnowAlb from ModelDailyState instead of SurfaceChar. Location also moved.
24 !Last modified: HCW 03 Jul 2015
25 ! Use PopDensNighttime by default (not PopDensDaytime)
26 !Last modified: HCW 26 Jun 2015
27 ! Translation of DailyState variables from the corresponding '_grids' arrays moved
28 ! earlier in code in order to fix bug in DecidCap, AlbDec, Porosity.
29 !Last modified: HCW 28 Nov 2014
30 !
31 ! To Do:
32 ! - Add AnOHM and ESTM info to FileChoices
33 ! - Check observed soil moisture works correctly!!
34 ! - Adjust model to allow water to runoff and sub-surface soil store for each surface type
35 ! - Adjust model to calculate LAI per surface
36 ! - Adjust model for SM per surface (measured characteristics)
37 !===================================================================================
38 SUBROUTINE suews_translate(Gridiv, ir, iMB)
42  USE data_in
43  USE defaultnotused, ONLY: notused, notusedi
44  USE gis_data, ONLY: &
47  USE mod_z, ONLY: z, z0m, z0m_in, zdm, zdm_in, zzd
48  USE resist, ONLY: g1, g2, g3, g4, g5, g6, th, tl, s1, s2, kmax, gsmodel
49  USE snowmod, ONLY: &
54  USE sues_data, ONLY: &
61  h_maintain, &
65 
66  USE time, ONLY: iy, id, it, imin, isec, dectime, dt_since_start
67  USE estm_data
68  USE wherewhen, ONLY: gridid, gridid_text
69  USE estm_module, ONLY: estm_translate
70 
71  IMPLICIT NONE
72 
73  INTEGER::Gridiv, & !Index of the analysed grid (Gridcounter)
74  ir, & !Meteorological forcing file index (set to zero if SUEWS_Translate called from InitialState)
75  iMB, & !Chunk of met data
76  id_prev
77 
78  INTEGER::iv, j, i
79  !real (Kind(1d0)):: FCskip = -9 !NULL value used for output to FileChoices
80  REAL(KIND(1d0)):: FCskip = -999 !NULL value used for output to FileChoices (changed by HCW 24 May 2016)
81 
82  ! REAL(KIND(1d0)):: z0m_in, zdm_in !Values of z0m and zdm provided in SiteSelect input file (do not get updated unlike z0d and z0m)
83 
84  CHARACTER(len=20):: grid_txt
85  CHARACTER(len=4):: year_txt
86  CHARACTER(len=12)::SsG_YYYY !Site, grid, year string
87 
88  CHARACTER(len=4):: iy_text
89  CHARACTER(len=3):: id_text
90  CHARACTER(len=2):: it_text, imin_text
91 
92  !write(*,*) '---- SUEWS_Translate ----'
93  !write(*,*) 'Year:', SurfaceChar(Gridiv,c_Year)
94  !write(*,*) 'Grid:', SurfaceChar(Gridiv,c_Grid)
95  !write(*,*) 'Gridiv:', Gridiv
96  !write(*,*) 'Met block (iMB or iv):',iMB
97  !write(*,*) 'Met line (ir):',ir
98  !write(*,*) '----'
99 
100  ! =================================================================================
101  ! ======= Translate inputs from SurfaceChar to variable names used in model =======
102  ! =================================================================================
103  ! GridID = GridIDmatrix(Gridiv) ! also in SUEWS_Program so deleted here. TS 10 Jun 2018
104  ! ---- Latitude and longitude
105  lat = surfacechar(gridiv, c_lat)
106  lng = surfacechar(gridiv, c_lng)
107  ! ---- Timezone
108  timezone = surfacechar(gridiv, c_tz)
109  ! ---- Altitude [m]
110  alt = surfacechar(gridiv, c_alt)
111  ! ---- Measurement height [m]
112  z = surfacechar(gridiv, c_z)
113  ! ---- Surface area [ha]
115  ! Change from ha to m2 (was in RunControlByGridByYear)
116  surfacearea = surfacearea_ha*10000 !Change surface area from ha to m^2
117 
118  ! ---- Surface fractions (previously in LUMPS_gis_read)
119  sfr(pavsurf) = surfacechar(gridiv, c_frpaved) ! Paved
120  sfr(bldgsurf) = surfacechar(gridiv, c_frbldgs) ! Bldgs
121  sfr(conifsurf) = surfacechar(gridiv, c_frevetr) ! Everg
122  sfr(decidsurf) = surfacechar(gridiv, c_frdectr) ! Decid
123  sfr(grasssurf) = surfacechar(gridiv, c_frgrass) ! Grass
124  sfr(bsoilsurf) = surfacechar(gridiv, c_frbsoil) ! BSoil
125  sfr(watersurf) = surfacechar(gridiv, c_frwater) ! Water
126 
127  ! Check the surface fractions add up to 1 (or close to 1)
128  IF (sum(sfr) > 1.001 .OR. sum(sfr) < 0.999) &
129  CALL errorhint(10, 'Surface fractions (Fr_) should add up to 1.', sum(sfr), notused, notusedi)
130 
131  ! ---- Irrigated fractions
132  irrfracpaved = surfacechar(gridiv, c_irrpavedfrac) ! Paved
133  irrfracbldgs = surfacechar(gridiv, c_irrbldgsfrac) ! Bldgs
134  irrfracevetr = surfacechar(gridiv, c_irrevetrfrac) ! Everg
135  irrfracdectr = surfacechar(gridiv, c_irrdectrfrac) ! Decid
136  irrfracgrass = surfacechar(gridiv, c_irrgrassfrac) ! Grass
137  irrfracbsoil = surfacechar(gridiv, c_irrbsoilfrac) ! BSoil
138  irrfracwater = surfacechar(gridiv, c_irrwaterfrac) ! Water
139 
140  ! ---------------------------------------------------------------------------------
141  ! --------- Surface cover calculations (previously in LUMPS_gis_read) -------------
142 
143  ! ---- Buildings and trees fraction ----
145 
146  ! ---- Vegetated fraction ----
148  !VegFraction = (sfr(ConifSurf) + sfr(DecidSurf) + sfr(GrassSurf))
149 
150  ! ---- Vegetated fraction (for LUMPS) ----
151  ! For LUMPS, vegetated fraction includes Water and Bare soil surfaces
152  IF (veg_type == 1) THEN ! area vegetated
154  ELSEIF (veg_type == 2) THEN ! area irrigated
156  ENDIF
157 
161  ! ---------------------------------------------------------------------------------
162 
163  ! ---- Heights & frontal areas
164  bldgh = surfacechar(gridiv, c_hbldgs) ! Building height [m]
165  evetreeh = surfacechar(gridiv, c_hevetr) ! Evergreen tree height [m]
166  dectreeh = surfacechar(gridiv, c_hdectr) ! Deciduous tree height [m]
167  IF (sfr(conifsurf) + sfr(decidsurf) > 0.) THEN ! avoid arithmetic error
168  treeh = (evetreeh*sfr(conifsurf) + dectreeh*sfr(decidsurf))/(sfr(conifsurf) + sfr(decidsurf)) ! Average tree height [m]
169  ELSE
170  treeh = 1.
171  END IF
172 
173  faibldg = surfacechar(gridiv, c_faibldgs) ! Frontal area index for buildings
174  faievetree = surfacechar(gridiv, c_faievetr) ! Frontal area index for evergreen trees
175  faidectree = surfacechar(gridiv, c_faidectr) ! Frontal area index for deciduous trees
176  IF (sfr(conifsurf) + sfr(decidsurf) > 0.) THEN ! avoid arithmetic error
177  faitree = (faievetree*sfr(conifsurf) + faidectree*sfr(decidsurf))/(sfr(conifsurf) + sfr(decidsurf)) ! Frontal area index for trees
178  ELSE
179  faitree = 1.
180  END IF
181 
182  z0m = surfacechar(gridiv, c_z0m) ! Roughness length [m]
183  zdm = surfacechar(gridiv, c_zdm) ! Displacement height [m]
184  ! z0m and zdm can vary in time depending on z0method selected. Save the input values here
185  z0m_in = z0m
186  zdm_in = zdm
187 
188  ! ---- Population density [ha-1]
189  ! Weekend fraction added to daytime population density
190  popdensdaytime = surfacechar(gridiv, (/c_popdensday, c_popdensday/)) ! Daytime population density [ha-1]
191  popdensnighttime = surfacechar(gridiv, c_popdensnight) ! Night-time population density [ha-1]
192  IF (popdensdaytime(1) >= 0 .AND. popdensnighttime < 0) popdensnighttime = popdensdaytime(1) !If only daytime data provided, use them
193  IF (popdensdaytime(1) < 0 .AND. popdensnighttime >= 0) popdensdaytime(1) = popdensnighttime !If only night-time data provided, use them
194  popdensdaytime(2) = popdensnighttime + (popdensdaytime(1) - popdensnighttime)*surfacechar(gridiv, c_frpddwe) !Use weekend fraction to daytime population
195  ! the following part has been moved into as can be derived there
196  ! IF (PopDensDaytime(1) >= 0 .AND. PopDensNighttime >= 0) NumCapita(1) = (PopDensDaytime(1) + PopDensNighttime)/2 !If both, use average
197  ! IF (PopDensDaytime(2) >= 0 .AND. PopDensNighttime >= 0) NumCapita(2) = (PopDensDaytime(2) + PopDensNighttime)/2 !If both, use average
198 
199  ! ! IF (PopDensDaytime >= 0 .AND. PopDensNighttime >= 0) NumCapita = (PopDensDaytime + PopDensNighttime)/2 !If both, use average ! moved to , TS 27 Dec 2018
200 
201  ! ---- Traffic rate
202  trafficrate = surfacechar(gridiv, (/c_trafficrate_wd, c_trafficrate_we/)) ! Mean traffic rate within modelled area
203  ! ---- Building energy use
204  qf0_beu = surfacechar(gridiv, (/c_qf0_beu_wd, c_qf0_beu_we/)) ! Building energy use within modelled area
205 
206  ! ---- Albedo [-]
207  alb(1:nsurf) = surfacechar(gridiv, c_albmax) !Use maximum albedos as default value (AlbMin for veg surfaces handled below)
208 
209  ! ---- Set min & max albedo for vegetated surfaces (min albedo not currently used for NonVeg or Water surfaces)
216 
217  ! ---- Emissivity [-]
218  emis(1:nsurf) = surfacechar(gridiv, c_emis)
219  emis_snow = surfacechar(gridiv, c_snowemis)
220 
221  ! ---- Storage capacities [mm]
222  storedrainprm(1, 1:nsurf) = surfacechar(gridiv, c_stormin) ! Minimum
223  storedrainprm(5, 1:nsurf) = surfacechar(gridiv, c_stormax) ! Maximum
224  storedrainprm(6, 1:nsurf) = storedrainprm(1, 1:nsurf) !Set storage capacities for all surface to minimum (DecTr changes with time in Calculations).
225 
226  ! ---- Set min & max storage capacities for DecTr
229  ! ---- Set min & max porosity for DecTr
230  pormin_dec = surfacechar(gridiv, c_porositymin(ivdecid)) ! Minimum
231  pormax_dec = surfacechar(gridiv, c_porositymax(ivdecid)) ! Minimum
232 
233  ! ---- Threshold for wet evaporation [mm]
234  wetthresh(1:nsurf) = surfacechar(gridiv, c_wetthresh)
235 
236  ! ---- Limit for state [mm]
238 
239  ! ---- Water depth [mm]
241 
242  ! ---- Drainage
243  storedrainprm(2, 1:nsurf) = surfacechar(gridiv, c_dreq) ! Drainage equation
244  storedrainprm(3, 1:nsurf) = surfacechar(gridiv, c_drcoef1) ! Drainage coef 1
245  storedrainprm(4, 1:nsurf) = surfacechar(gridiv, c_drcoef2) ! Drainage coef 2
246 
247  ! ---- Limit of SWE (each surface except Water)
248  snowpacklimit(1:(nsurf - 1)) = surfacechar(gridiv, c_snowlimpat(1:(nsurf - 1)))
249 
250  ! ---- Snow limit for removal (only impervious surfaces)
253  !SnowLimBSoil = SurfaceChar(Gridiv,c_SnowLimRem(BSoilSurf)) !Snow clearing not applicable to bare soil surface
254 
255  ! ---- Soil characteristics (each surface except Water)
256  soildepth(1:(nsurf - 1)) = surfacechar(gridiv, c_soildepth(1:(nsurf - 1))) ! Depth of sub-surface soil store [mm]
257  soilstorecap(1:(nsurf - 1)) = surfacechar(gridiv, c_soilstcap(1:(nsurf - 1))) ! Soil store capacity [mm]
258  sathydraulicconduct(1:(nsurf - 1)) = surfacechar(gridiv, c_ksat(1:(nsurf - 1))) ! Hydraulic conductivity of saturated soil [mm s-1]
259  !SoilDensity(1:(nsurf-1)) = SurfaceChar(Gridiv,c_SoilDens(1:(nsurf-1))) ! Soil density [kg m-3]
260  ! Not yet implemented in model
261  !InfiltrationRate (1:(nsurf-1)) = SurfaceChar(Gridiv,c_SoilInfRate(1:(nsurf-1))) ! Infiltration rate [mm h-1]
262 
263  !! Observed soil characteristics
264  !SoilDensity (1:(nsurf-1)) = SurfaceChar(Gridiv,c_SoilDens(1:(nsurf-1))) ! Soil density [kg m-3]
265  !SoilDepthMeas(1:(nsurf-1)) = SurfaceChar(Gridiv,c_ObsSMDepth(1:(nsurf-1)))
266  !SmCap (1:(nsurf-1)) = SurfaceChar(Gridiv,c_ObsSMMax(1:(nsurf-1)))
267  !SoilRocks (1:(nsurf-1)) = SurfaceChar(Gridiv,c_ObsSNRFrac(1:(nsurf-1)))
268  !!Obs soil characteristics now in SUEWS_Soil, i.e. per surface; single value was given previously in FunctionalTypes
269  !!Take first row here for testing !! Need to alter model later...
270  soildensity = surfacechar(gridiv, c_soildens(1)) !!Not sure this works correctly - need to check
272  smcap = surfacechar(gridiv, c_obssmmax(1))
273  soilrocks = surfacechar(gridiv, c_obssnrfrac(1))
274 
275  ! ---- Vegetation characteristics (pervious surfaces)
276  baset(1:nvegsurf) = surfacechar(gridiv, c_baset)
277  basete(1:nvegsurf) = surfacechar(gridiv, c_basete)
278  gddfull(1:nvegsurf) = surfacechar(gridiv, c_gddfull)
279  sddfull(1:nvegsurf) = surfacechar(gridiv, c_sddfull)
280  laimin(1:nvegsurf) = surfacechar(gridiv, c_laimin)
281  laimax(1:nvegsurf) = surfacechar(gridiv, c_laimax)
283 
289  resp_a(1:nvegsurf) = surfacechar(gridiv, c_resp_a)
290  resp_b(1:nvegsurf) = surfacechar(gridiv, c_resp_b)
292 
293  ! ---- LAI characteristics (updated HCW 13 Dec 2016)
294  laitype(1:nvegsurf) = int(surfacechar(gridiv, c_laieq(1:nvegsurf)))
295  laipower(1, 1:nvegsurf) = surfacechar(gridiv, c_leafgp1(1:nvegsurf))
296  laipower(2, 1:nvegsurf) = surfacechar(gridiv, c_leafgp2(1:nvegsurf))
297  laipower(3, 1:nvegsurf) = surfacechar(gridiv, c_leafop1(1:nvegsurf))
298  laipower(4, 1:nvegsurf) = surfacechar(gridiv, c_leafop2(1:nvegsurf))
299 
300  ! ---- LUMPS-related parameters
301  drainrt = surfacechar(gridiv, c_lumpsdr) ! LUMPS Drainage rate [mm h-1]
302  raincover = surfacechar(gridiv, c_lumpscover) ! LUMPS Limit when surface totally wet [mm]
303  rainmaxres = surfacechar(gridiv, c_lumpsmaxres) ! LUMPS Maximum water bucket reservoir [mm]
304 
305  ! ---- NARP-related parameters
306  trans_site = surfacechar(gridiv, c_narptrans) ! NARP atmospheric transmissivity
307 
308  ! ---- Snow-related characteristics
313  tau_a = surfacechar(gridiv, c_snowtau_a)
314  tau_f = surfacechar(gridiv, c_snowtau_f)
318  tau_r = surfacechar(gridiv, c_snowtau_r)
319  crwmin = surfacechar(gridiv, c_snowcrwmin)
320  crwmax = surfacechar(gridiv, c_snowcrwmax)
322 
323  ! ---- Conductance parameters
324  g1 = surfacechar(gridiv, c_gsg1)
325  g2 = surfacechar(gridiv, c_gsg2)
326  g3 = surfacechar(gridiv, c_gsg3)
327  g4 = surfacechar(gridiv, c_gsg4)
328  g5 = surfacechar(gridiv, c_gsg5)
329  g6 = surfacechar(gridiv, c_gsg6)
330  th = surfacechar(gridiv, c_gsth)
331  tl = surfacechar(gridiv, c_gstl)
332  s1 = surfacechar(gridiv, c_gss1)
333  s2 = surfacechar(gridiv, c_gss2)
334  kmax = surfacechar(gridiv, c_gskmax)
335  gsmodel = int(surfacechar(gridiv, c_gsmodel))
336 
337  ! ---- Pipe capacity (was from SiteSpecificParam.txt)
339 
340  ! ---- Water flows (was from SiteSpecificParam.txt)
343 
344  ! ---- Daylight saving (was from ModelledYears.txt)
345  startdls = int(surfacechar(gridiv, c_startdls))
346  endDLS = INT(SurfaceChar(Gridiv, c_EndDLS))
347 
348  ! ---- OHM coeffs (was in SUEWS_OHMnew.f95, subroutine OHMinitialize)
349  ohm_coef = 0 ! Initialise OHM_coef
350  ! Surface types in OHM_coef: Paved, Roof, Conif, Decid, Grass, BareSoil, Water, CANYON, Snow
351  ! No canyon in SurfaceChar, so
352  ! transfer coeffs for surface types 1-7,
353  ! then skip row in OHM_coef (canyon),
354  ! then transfer coeffs for snow surface (8th surface in SurfaceChar; 9th surface in OHM_Coefs)
355  ! Summer wet
356  ohm_coef(1:nsurf, 1, 1) = surfacechar(gridiv, c_a1_swet(1:nsurf)) !1:nsurf a1 Summer wet
357  ohm_coef(nsurf + 1, 1, 1) = surfacechar(gridiv, c_a1_swet(nsurf + 1)) !Snow a1 Summer wet
358  ohm_coef(1:nsurf, 1, 2) = surfacechar(gridiv, c_a2_swet(1:nsurf)) !1:nsurf a2 Summer wet
359  ohm_coef(nsurf + 1, 1, 2) = surfacechar(gridiv, c_a2_swet(nsurf + 1)) !Snow a2 Summer wet
360  ohm_coef(1:nsurf, 1, 3) = surfacechar(gridiv, c_a3_swet(1:nsurf)) !1:nsurf a3 Summer wet
361  ohm_coef(nsurf + 1, 1, 3) = surfacechar(gridiv, c_a3_swet(nsurf + 1)) !Snow a3 Summer wet
362  ! Summer dry
363  ohm_coef(1:nsurf, 2, 1) = surfacechar(gridiv, c_a1_sdry(1:nsurf)) !1:nsurf a1 Summer dry
364  ohm_coef(nsurf + 1, 2, 1) = surfacechar(gridiv, c_a1_sdry(nsurf + 1)) !Snow a1 Summer dry
365  ohm_coef(1:nsurf, 2, 2) = surfacechar(gridiv, c_a2_sdry(1:nsurf)) !1:nsurf a2 Summer dry
366  ohm_coef(nsurf + 1, 2, 2) = surfacechar(gridiv, c_a2_sdry(nsurf + 1)) !Snow a2 Summer dry
367  ohm_coef(1:nsurf, 2, 3) = surfacechar(gridiv, c_a3_sdry(1:nsurf)) !1:nsurf a3 Summer dry
368  ohm_coef(nsurf + 1, 2, 3) = surfacechar(gridiv, c_a3_sdry(nsurf + 1)) !Snow a3 Summer dry
369  ! Winter wet
370  ohm_coef(1:nsurf, 3, 1) = surfacechar(gridiv, c_a1_wwet(1:nsurf)) !1:nsurf a1 Winter wet
371  ohm_coef(nsurf + 1, 3, 1) = surfacechar(gridiv, c_a1_wwet(nsurf + 1)) !Snow a1 Winter wet
372  ohm_coef(1:nsurf, 3, 2) = surfacechar(gridiv, c_a2_wwet(1:nsurf)) !1:nsurf a2 Winter wet
373  ohm_coef(nsurf + 1, 3, 2) = surfacechar(gridiv, c_a2_wwet(nsurf + 1)) !Snow a2 Winter wet
374  ohm_coef(1:nsurf, 3, 3) = surfacechar(gridiv, c_a3_wwet(1:nsurf)) !1:nsurf a3 Winter wet
375  ohm_coef(nsurf + 1, 3, 3) = surfacechar(gridiv, c_a3_wwet(nsurf + 1)) !Snow a3 Winter wet
376  ! Winter dry
377  ohm_coef(1:nsurf, 4, 1) = surfacechar(gridiv, c_a1_wdry(1:nsurf)) !1:nsurf a1 Winter dry
378  ohm_coef(nsurf + 1, 4, 1) = surfacechar(gridiv, c_a1_wdry(nsurf + 1)) !Snow a1 Winter dry
379  ohm_coef(1:nsurf, 4, 2) = surfacechar(gridiv, c_a2_wdry(1:nsurf)) !1:nsurf a2 Winter dry
380  ohm_coef(nsurf + 1, 4, 2) = surfacechar(gridiv, c_a2_wdry(nsurf + 1)) !Snow a2 Winter dry
381  ohm_coef(1:nsurf, 4, 3) = surfacechar(gridiv, c_a3_wdry(1:nsurf)) !1:nsurf a3 Winter dry
382  ohm_coef(nsurf + 1, 4, 3) = surfacechar(gridiv, c_a3_wdry(nsurf + 1)) !Snow a3 Winter dry
383  ! OHM thresholds
384  ohm_threshsw(1:nsurf) = surfacechar(gridiv, c_ohmthresh_sw(1:nsurf)) !1:nsurf
385  ohm_threshsw(nsurf + 1) = surfacechar(gridiv, c_ohmthresh_sw(nsurf + 1)) !Snow
386  ohm_threshwd(1:nsurf) = surfacechar(gridiv, c_ohmthresh_wd(1:nsurf)) !1:nsurf
387  ohm_threshwd(nsurf + 1) = surfacechar(gridiv, c_ohmthresh_wd(nsurf + 1)) !Snow
388 
389  ! ---- ESTM characteristics -------------------------
390  ! HCW 16 Jun 2016
391  ! Wall fraction for ESTM (in SiteSelect.txt)
392  IF (storageheatmethod == 4 .OR. storageheatmethod == 14) THEN
393  areawall = surfacechar(gridiv, c_areawall)
394  fwall = areawall/surfacearea
395 
396  ! Get surface fractions for ESTM classes for Bldgs and Paved surfaces
397  estmsfr_paved = surfacechar(gridiv, c_fr_estmclass_paved) !Dim 3
398  estmsfr_bldgs = surfacechar(gridiv, c_fr_estmclass_bldgs) !Dim 5
399  !Check these sum to 1 and are consistent with sfr of Paved and Bldgs surface types
400  IF (sfr(pavsurf) > 0) THEN !If surface exists, ESTM fractions must be correct
401  IF (sum(estmsfr_paved) > 1.001 .OR. sum(estmsfr_paved) < 0.999) THEN
402  CALL errorhint(10, 'Surface fractions (Fr_ESTMClass_Paved) should sum to 1.', sum(estmsfr_paved), notused, notusedi)
403  ENDIF
404  ELSEIF (sfr(pavsurf) == 0) THEN !If surface does not exist, ESTM fraction does not matter
405  IF (sum(estmsfr_paved) > 1.001 .OR. sum(estmsfr_paved) < 0.999) THEN !If ESTM fractions do not sum to 1, set here
406  estmsfr_paved(1) = 1.000
407  estmsfr_paved(2:3) = 0.000
408  CALL errorhint(67, 'ESTM Paved classes do not sum to 1 (but no Paved surface present).', &
409  sum(estmsfr_paved), notused, notusedi)
410  ENDIF
411  ENDIF
412  IF (sfr(bldgsurf) > 0) THEN
413  IF (sum(estmsfr_bldgs) > 1.001 .OR. sum(estmsfr_bldgs) < 0.999) THEN
414  CALL errorhint(10, 'Surface fractions (Fr_ESTMClass_Bldgs) should sum to 1.', sum(estmsfr_bldgs), notused, notusedi)
415  ENDIF
416  ELSEIF (sfr(bldgsurf) == 0) THEN !If surface does not exist, ESTM fraction does not matter
417  IF (sum(estmsfr_bldgs) > 1.001 .OR. sum(estmsfr_bldgs) < 0.999) THEN !If ESTM fractions do not sum to 1, set here
418  estmsfr_bldgs(1) = 1.000
419  estmsfr_bldgs(2:5) = 0.000
420  CALL errorhint(67, 'ESTM Bldgs classes do not sum to 1 (but no Bldgs surface present).', &
421  sum(estmsfr_bldgs), notused, notusedi)
422  ENDIF
423  ENDIF
424 
425  ! ===== PAVED =====
426  ! First combine characteristics of the 3x Paved classes
427  IF (surfacechar(gridiv, c_estmcode(pavsurf)) == 0) THEN ! If Code = 0, use multiple classes
428  ! Get characteristics of each Paved class
429  DO i = 1, 3
430  zsurf_paved(:, i) = surfacechar(gridiv, (/c_surf_thick1_paved(i), c_surf_thick2_paved(i), c_surf_thick3_paved(i), &
431  c_surf_thick4_paved(i), c_surf_thick5_paved(i)/))
432  ksurf_paved(:, i) = surfacechar(gridiv, (/c_surf_k1_paved(i), c_surf_k2_paved(i), c_surf_k3_paved(i), &
433  c_surf_k4_paved(i), c_surf_k5_paved(i)/))
434  rsurf_paved(:, i) = surfacechar(gridiv, (/c_surf_rhocp1_paved(i), c_surf_rhocp2_paved(i), c_surf_rhocp3_paved(i), &
435  c_surf_rhocp4_paved(i), c_surf_rhocp5_paved(i)/))
436  ENDDO
437  ! Average characteristics of each Paved class according to surface fractions (these sum to 1)
438  zsurf_suewssurfs(:, pavsurf) = zsurf_paved(:, 1)*estmsfr_paved(1) &
439  + zsurf_paved(:, 2)*estmsfr_paved(2) &
440  + zsurf_paved(:, 3)*estmsfr_paved(3)
441  ksurf_suewssurfs(:, pavsurf) = ksurf_paved(:, 1)*estmsfr_paved(1) &
442  + ksurf_paved(:, 2)*estmsfr_paved(2) &
443  + ksurf_paved(:, 3)*estmsfr_paved(3)
444  rsurf_suewssurfs(:, pavsurf) = rsurf_paved(:, 1)*estmsfr_paved(1) &
445  + rsurf_paved(:, 2)*estmsfr_paved(2) &
446  + rsurf_paved(:, 3)*estmsfr_paved(3)
447  ELSEIF (surfacechar(gridiv, c_estmcode(pavsurf)) /= 0) THEN !Otherwise use single values
448  zsurf_suewssurfs(:, pavsurf) = surfacechar(gridiv, &
449  (/c_surf_thick1(pavsurf), c_surf_thick2(pavsurf), c_surf_thick3(pavsurf), &
450  c_surf_thick4(pavsurf), c_surf_thick5(pavsurf)/))
451  ksurf_suewssurfs(:, pavsurf) = surfacechar(gridiv, &
452  (/c_surf_k1(pavsurf), c_surf_k2(pavsurf), c_surf_k3(pavsurf), &
453  c_surf_k4(pavsurf), c_surf_k5(pavsurf)/))
454  rsurf_suewssurfs(:, pavsurf) = surfacechar(gridiv, &
455  (/c_surf_rhocp1(pavsurf), c_surf_rhocp2(pavsurf), c_surf_rhocp3(pavsurf), &
456  c_surf_rhocp4(pavsurf), c_surf_rhocp5(pavsurf)/))
457  ENDIF
458 
459  ! ===== BLDGS =====
460  ! Combine characteristics of 5x Bldgs classes into one
461  IF (surfacechar(gridiv, c_estmcode(bldgsurf)) == 0) THEN ! If Code = 0, use multiple classes
462  ! Get characteristics of each Bldgs class
463  DO i = 1, 5
464  zsurf_bldgs(:, i) = surfacechar(gridiv, (/c_surf_thick1_bldgs(i), c_surf_thick2_bldgs(i), c_surf_thick3_bldgs(i), &
465  c_surf_thick4_bldgs(i), c_surf_thick5_bldgs(i)/))
466  ksurf_bldgs(:, i) = surfacechar(gridiv, (/c_surf_k1_bldgs(i), c_surf_k2_bldgs(i), c_surf_k3_bldgs(i), &
467  c_surf_k4_bldgs(i), c_surf_k5_bldgs(i)/))
468  rsurf_bldgs(:, i) = surfacechar(gridiv, (/c_surf_rhocp1_bldgs(i), c_surf_rhocp2_bldgs(i), c_surf_rhocp3_bldgs(i), &
469  c_surf_rhocp4_bldgs(i), c_surf_rhocp5_bldgs(i)/))
470  zwall_bldgs(:, i) = surfacechar(gridiv, (/c_wall_thick1_bldgs(i), c_wall_thick2_bldgs(i), c_wall_thick3_bldgs(i), &
471  c_wall_thick4_bldgs(i), c_wall_thick5_bldgs(i)/))
472  kwall_bldgs(:, i) = surfacechar(gridiv, (/c_wall_k1_bldgs(i), c_wall_k2_bldgs(i), c_wall_k3_bldgs(i), &
473  c_wall_k4_bldgs(i), c_wall_k5_bldgs(i)/))
474  rwall_bldgs(:, i) = surfacechar(gridiv, (/c_wall_rhocp1_bldgs(i), c_wall_rhocp2_bldgs(i), c_wall_rhocp3_bldgs(i), &
475  c_wall_rhocp4_bldgs(i), c_wall_rhocp5_bldgs(i)/))
476  zibld_bldgs(:, i) = surfacechar(gridiv, (/c_internal_thick1_bldgs(i), c_internal_thick2_bldgs(i), &
477  c_internal_thick3_bldgs(i), &
478  c_internal_thick4_bldgs(i), c_internal_thick5_bldgs(i)/))
479  kibld_bldgs(:, i) = surfacechar(gridiv, (/c_internal_k1_bldgs(i), c_internal_k2_bldgs(i), c_internal_k3_bldgs(i), &
480  c_internal_k4_bldgs(i), c_internal_k5_bldgs(i)/))
481  ribld_bldgs(:, i) = surfacechar(gridiv, (/c_internal_rhocp1_bldgs(i), c_internal_rhocp2_bldgs(i), &
482  c_internal_rhocp3_bldgs(i), &
483  c_internal_rhocp4_bldgs(i), c_internal_rhocp5_bldgs(i)/))
484  nroom_bldgs(i) = surfacechar(gridiv, c_nroom_bldgs(i))
485  alb_ibld_bldgs(i) = surfacechar(gridiv, c_alb_ibld_bldgs(i))
486  em_ibld_bldgs(i) = surfacechar(gridiv, c_em_ibld_bldgs(i))
487  ch_iwall_bldgs(i) = surfacechar(gridiv, c_ch_iwall_bldgs(i))
488  ch_iroof_bldgs(i) = surfacechar(gridiv, c_ch_iroof_bldgs(i))
489  ch_ibld_bldgs(i) = surfacechar(gridiv, c_ch_ibld_bldgs(i))
490  ENDDO
491  ! Average characteristics of each Bldgs class according to surface fractions (these sum to 1)
492  zsurf_suewssurfs(:, bldgsurf) = zsurf_bldgs(:, 1)*estmsfr_bldgs(1) &
493  + zsurf_bldgs(:, 2)*estmsfr_bldgs(2) &
494  + zsurf_bldgs(:, 3)*estmsfr_bldgs(3) &
495  + zsurf_bldgs(:, 4)*estmsfr_bldgs(4) &
496  + zsurf_bldgs(:, 5)*estmsfr_bldgs(5)
497  ksurf_suewssurfs(:, bldgsurf) = ksurf_bldgs(:, 1)*estmsfr_bldgs(1) &
498  + ksurf_bldgs(:, 2)*estmsfr_bldgs(2) &
499  + ksurf_bldgs(:, 3)*estmsfr_bldgs(3) &
500  + ksurf_bldgs(:, 4)*estmsfr_bldgs(4) &
501  + ksurf_bldgs(:, 5)*estmsfr_bldgs(5)
502  rsurf_suewssurfs(:, bldgsurf) = rsurf_bldgs(:, 1)*estmsfr_bldgs(1) &
503  + rsurf_bldgs(:, 2)*estmsfr_bldgs(2) &
504  + rsurf_bldgs(:, 3)*estmsfr_bldgs(3) &
505  + rsurf_bldgs(:, 4)*estmsfr_bldgs(4) &
506  + rsurf_bldgs(:, 5)*estmsfr_bldgs(5)
507  !Wall
508  zwall = zwall_bldgs(:, 1)*estmsfr_bldgs(1) &
509  + zwall_bldgs(:, 2)*estmsfr_bldgs(2) &
510  + zwall_bldgs(:, 3)*estmsfr_bldgs(3) &
511  + zwall_bldgs(:, 4)*estmsfr_bldgs(4) &
512  + zwall_bldgs(:, 5)*estmsfr_bldgs(5)
513  kwall = kwall_bldgs(:, 1)*estmsfr_bldgs(1) &
514  + kwall_bldgs(:, 2)*estmsfr_bldgs(2) &
515  + kwall_bldgs(:, 3)*estmsfr_bldgs(3) &
516  + kwall_bldgs(:, 4)*estmsfr_bldgs(4) &
517  + kwall_bldgs(:, 5)*estmsfr_bldgs(5)
518  rwall = rwall_bldgs(:, 1)*estmsfr_bldgs(1) &
519  + rwall_bldgs(:, 2)*estmsfr_bldgs(2) &
520  + rwall_bldgs(:, 3)*estmsfr_bldgs(3) &
521  + rwall_bldgs(:, 4)*estmsfr_bldgs(4) &
522  + rwall_bldgs(:, 5)*estmsfr_bldgs(5)
523  !Internal
524  zibld = zibld_bldgs(:, 1)*estmsfr_bldgs(1) &
525  + zibld_bldgs(:, 2)*estmsfr_bldgs(2) &
526  + zibld_bldgs(:, 3)*estmsfr_bldgs(3) &
527  + zibld_bldgs(:, 4)*estmsfr_bldgs(4) &
528  + zibld_bldgs(:, 5)*estmsfr_bldgs(5)
529  kibld = kibld_bldgs(:, 1)*estmsfr_bldgs(1) &
530  + kibld_bldgs(:, 2)*estmsfr_bldgs(2) &
531  + kibld_bldgs(:, 3)*estmsfr_bldgs(3) &
532  + kibld_bldgs(:, 4)*estmsfr_bldgs(4) &
533  + kibld_bldgs(:, 5)*estmsfr_bldgs(5)
534  ribld = ribld_bldgs(:, 1)*estmsfr_bldgs(1) &
535  + ribld_bldgs(:, 2)*estmsfr_bldgs(2) &
536  + ribld_bldgs(:, 3)*estmsfr_bldgs(3) &
537  + ribld_bldgs(:, 4)*estmsfr_bldgs(4) &
538  + ribld_bldgs(:, 5)*estmsfr_bldgs(5)
539 
540  nroom = nroom_bldgs(1)*estmsfr_bldgs(1) &
541  + nroom_bldgs(2)*estmsfr_bldgs(2) &
542  + nroom_bldgs(3)*estmsfr_bldgs(3) &
543  + nroom_bldgs(4)*estmsfr_bldgs(4) &
544  + nroom_bldgs(5)*estmsfr_bldgs(5)
545  alb_ibld = alb_ibld_bldgs(1)*estmsfr_bldgs(1) &
546  + alb_ibld_bldgs(2)*estmsfr_bldgs(2) &
547  + alb_ibld_bldgs(3)*estmsfr_bldgs(3) &
548  + alb_ibld_bldgs(4)*estmsfr_bldgs(4) &
549  + alb_ibld_bldgs(5)*estmsfr_bldgs(5)
550  em_ibld = em_ibld_bldgs(1)*estmsfr_bldgs(1) &
551  + em_ibld_bldgs(2)*estmsfr_bldgs(2) &
552  + em_ibld_bldgs(3)*estmsfr_bldgs(3) &
553  + em_ibld_bldgs(4)*estmsfr_bldgs(4) &
554  + em_ibld_bldgs(5)*estmsfr_bldgs(5)
555  ch_iwall = ch_iwall_bldgs(1)*estmsfr_bldgs(1) &
556  + ch_iwall_bldgs(2)*estmsfr_bldgs(2) &
557  + ch_iwall_bldgs(3)*estmsfr_bldgs(3) &
558  + ch_iwall_bldgs(4)*estmsfr_bldgs(4) &
559  + ch_iwall_bldgs(5)*estmsfr_bldgs(5)
560  ch_iroof = ch_iroof_bldgs(1)*estmsfr_bldgs(1) &
561  + ch_iroof_bldgs(2)*estmsfr_bldgs(2) &
562  + ch_iroof_bldgs(3)*estmsfr_bldgs(3) &
563  + ch_iroof_bldgs(4)*estmsfr_bldgs(4) &
564  + ch_iroof_bldgs(5)*estmsfr_bldgs(5)
565  ch_ibld = ch_ibld_bldgs(1)*estmsfr_bldgs(1) &
566  + ch_ibld_bldgs(2)*estmsfr_bldgs(2) &
567  + ch_ibld_bldgs(3)*estmsfr_bldgs(3) &
568  + ch_ibld_bldgs(4)*estmsfr_bldgs(4) &
569  + ch_ibld_bldgs(5)*estmsfr_bldgs(5)
570 
571  ELSEIF (surfacechar(gridiv, c_estmcode(bldgsurf)) /= 0) THEN !Otherwise use single values
572  zsurf_suewssurfs(:, bldgsurf) = surfacechar(gridiv, (/c_surf_thick1(bldgsurf), c_surf_thick2(bldgsurf), &
573  c_surf_thick3(bldgsurf), &
574  c_surf_thick4(bldgsurf), c_surf_thick5(bldgsurf)/))
575  ksurf_suewssurfs(:, bldgsurf) = surfacechar(gridiv, (/c_surf_k1(bldgsurf), c_surf_k2(bldgsurf), c_surf_k3(bldgsurf), &
576  c_surf_k4(bldgsurf), c_surf_k5(bldgsurf)/))
577  rsurf_suewssurfs(:, bldgsurf) = surfacechar(gridiv, (/c_surf_rhocp1(bldgsurf), c_surf_rhocp2(bldgsurf), &
578  c_surf_rhocp3(bldgsurf), &
579  c_surf_rhocp4(bldgsurf), c_surf_rhocp5(bldgsurf)/))
580  zwall = surfacechar(gridiv, (/c_wall_thick1, c_wall_thick2, c_wall_thick3, c_wall_thick4, c_wall_thick5/))
581  kwall = surfacechar(gridiv, (/c_wall_k1, c_wall_k2, c_wall_k3, c_wall_k4, c_wall_k5/))
582  rwall = surfacechar(gridiv, (/c_wall_rhocp1, c_wall_rhocp2, c_wall_rhocp3, c_wall_rhocp4, c_wall_rhocp5/))
583  zibld = surfacechar(gridiv, &
584  (/c_internal_thick1, c_internal_thick2, c_internal_thick3, c_internal_thick4, c_internal_thick5/))
585  kibld = surfacechar(gridiv, (/c_internal_k1, c_internal_k2, c_internal_k3, c_internal_k4, c_internal_k5/))
586  ribld = surfacechar(gridiv, &
587  (/c_internal_rhocp1, c_internal_rhocp2, c_internal_rhocp3, c_internal_rhocp4, c_internal_rhocp5/))
588 
589  nroom = surfacechar(gridiv, c_nroom)
590  alb_ibld = surfacechar(gridiv, c_alb_ibld)
591  em_ibld = surfacechar(gridiv, c_em_ibld)
592  ch_iwall = surfacechar(gridiv, c_ch_iwall)
593  ch_iroof = surfacechar(gridiv, c_ch_iroof)
594  ch_ibld = surfacechar(gridiv, c_ch_ibld)
595  ENDIF
596 
597  !For other surfaces, only one ESTM class
598  DO iv = conifsurf, nsurfincsnow
599  zsurf_suewssurfs(:, iv) = surfacechar(gridiv, (/c_surf_thick1(iv), c_surf_thick2(iv), c_surf_thick3(iv), &
600  c_surf_thick4(iv), c_surf_thick5(iv)/))
601  ksurf_suewssurfs(:, iv) = surfacechar(gridiv, (/c_surf_k1(iv), c_surf_k2(iv), c_surf_k3(iv), &
602  c_surf_k4(iv), c_surf_k5(iv)/))
603  rsurf_suewssurfs(:, iv) = surfacechar(gridiv, (/c_surf_rhocp1(iv), c_surf_rhocp2(iv), c_surf_rhocp3(iv), &
604  c_surf_rhocp4(iv), c_surf_rhocp5(iv)/))
605  ENDDO
606 
607  ! Now combine SUEWS surfaces into ESTM facets
608  !Surface fractions for ESTM facets (moved from SUEWS_ESTM_initials HCW 16 Jun 2016)
609  !roof = Bldgs
610  froof = sfr(bldgsurf)
611  !ground = all except Bldgs
612  fground = sfr(pavsurf) + sfr(conifsurf) + sfr(decidsurf) + sfr(grasssurf) + sfr(bsoilsurf) + sfr(watersurf)
613  !veg = EveTr, DecTr, Grass
614  fveg = sfr(conifsurf) + sfr(decidsurf) + sfr(grasssurf)
615 
616  ! Ground = all except buildings (exclude snow at the moment)
617  zground = 0
618  kground = 0
619  rground = 0
620  DO iv = 1, nsurf
621  IF (iv /= bldgsurf .AND. fground /= 0) THEN !Bldgs surface excluded from ground facet
622  zground = zground + zsurf_suewssurfs(:, iv)*sfr(iv)/fground !Normalised by ground fraction
623  kground = kground + ksurf_suewssurfs(:, iv)*sfr(iv)/fground !Normalised by ground fraction
624  rground = rground + rsurf_suewssurfs(:, iv)*sfr(iv)/fground !Normalised by ground fraction
625  ELSEIF (fground == 0.) THEN !check fground==0 (or HW==0) scenario to avoid division-by-zero error, TS 21 Jul 2016
626  zground = zground + 0.01
627  kground = kground + 0.01
628  rground = rground + 0.01
629  ! PRINT*, zground
630  ! PRINT*, kground
631  ! PRINT*, rground
632  ENDIF
633  ENDDO
634  ! Roof = buildings
635  zroof = zsurf_suewssurfs(:, bldgsurf)
636  kroof = ksurf_suewssurfs(:, bldgsurf)
637  rroof = rsurf_suewssurfs(:, bldgsurf)
638 
639  ! the following initialisation is problematic: TS 01 Mar 2019
640  ! what would happen if zground(5)>0? Nground is initialised NOWHERE!
641  ! initialise these variables as 5 so if z_sfc(5)>0 happens, the numbers are still correct, TS 06 Aug 2019
642  nground = 5
643  nroof = 5
644  nwall = 5
645  nibld = 5
646  DO i = 1, 5
647  IF (zground(i) <= 0) THEN
648  nground = i - 1
649  EXIT
650  ENDIF
651  ENDDO
652  DO i = 1, 5
653  IF (zroof(i) <= 0) THEN
654  nroof = i - 1
655  EXIT
656  ENDIF
657  ENDDO
658  DO i = 1, 5
659  IF (zwall(i) <= 0) THEN
660  nwall = i - 1
661  EXIT
662  ENDIF
663  ENDDO
664  DO i = 1, 5
665  IF (zibld(i) <= 0) THEN
666  nibld = i - 1
667  EXIT
668  ENDIF
669  ENDDO
670  ENDIF ! ESTM related translation finished here.
671 
672  ! ---- AnOHM related ------------------------------
673  IF (storageheatmethod == 3) THEN
674  cpanohm(1:nsurf) = surfacechar(gridiv, c_cpanohm) ! AnOHM TS
675  kkanohm(1:nsurf) = surfacechar(gridiv, c_kkanohm) ! AnOHM TS
676  chanohm(1:nsurf) = surfacechar(gridiv, c_chanohm) ! AnOHM TS
677 
678  ! cp and k are estimated from ESTM coefficients:
679  ! cpAnOHM(1:nsurf)=rSurf_SUEWSsurfs(1,1:nsurf)
680  ! kkAnOHM(1:nsurf)=kSurf_SUEWSsurfs(1,1:nsurf)
681  ! IF ( ir ==1 .AND. iMb ==1) THEN
682  ! PRINT*, 'StoreDrainPrm',PavSurf,':'
683  ! PRINT'(a10,x,5f10.2)', 'Depth',zSurf_SUEWSsurfs(:,i)
684  ! PRINT'(a10,x,5es10.2)', 'RhoCp',rSurf_SUEWSsurfs(:,i)
685  ! PRINT'(a10,x,5es10.2)', 'avg_RhoCp',cpAnOHM(i)
686  ! PRINT'(a10,x,5es10.2)', 'k',kSurf_SUEWSsurfs(:,i)
687  ! PRINT'(a10,x,5es10.2)', 'avg_k',kkAnOHM(i)
688  ! PRINT'(a10,x,5es10.2)', 'avg_Ch',chAnOHM(i)
689  !
690  ! END IF
691  ! DO i = 1, nsurf, 1
692  ! ! filter out invalid z values
693  ! WHERE ( zSurf_SUEWSsurfs(:,i) == -999. ) zSurf_SUEWSsurfs(:,i)=0
694  !
695  ! ! cp: weight-averaged by depth
696  ! cpAnOHM(i)=DOT_PRODUCT(rSurf_SUEWSsurfs(:,i),zSurf_SUEWSsurfs(:,i))/SUM(zSurf_SUEWSsurfs(:,i))
697  ! ! IF ( i==PavSurf .AND. ir ==1 .AND. iMb ==1) THEN
698  ! ! PRINT*, 'StoreDrainPrm',i,':'
699  ! ! PRINT'(a10,x,5f10.2)', 'Depth',zSurf_SUEWSsurfs(:,i)
700  ! ! PRINT'(a10,x,5es10.2)', 'RhoCp',rSurf_SUEWSsurfs(:,i)
701  ! ! PRINT'(a10,x,5es10.2)', 'avg_RhoCp',cpAnOHM(i)
702  ! !
703  ! ! END IF
704  !
705  ! ! 1/k: weight-averaged by depth
706  ! kkAnOHM(i)=DOT_PRODUCT(1/kSurf_SUEWSsurfs(:,i),zSurf_SUEWSsurfs(:,i))/SUM(zSurf_SUEWSsurfs(:,i))
707  ! kkAnOHM(i)=1/kkAnOHM(i)
708  ! ! IF ( i==PavSurf .AND. ir ==1 .AND. iMb ==1) THEN
709  ! ! PRINT'(a10,x,5es10.2)', 'k',kSurf_SUEWSsurfs(:,i)
710  ! ! PRINT'(a10,x,5es10.2)', 'avg_k',kkAnOHM(i)
711  ! ! PRINT'(a10,x,5es10.2)', 'avg_Ch',chAnOHM(i)
712  ! !
713  ! ! PRINT'(a10,x,7f10.2)', 'fractions:',sfr
714  ! !
715  ! ! END IF
716  !
717  !
718  ! ! restore invalid z values
719  ! WHERE ( zSurf_SUEWSsurfs(:,i) == 0 ) zSurf_SUEWSsurfs(:,i)=nan
720  ! ! IF ( i==PavSurf .AND. ir ==1 .AND. iMb ==1) THEN
721  ! ! PRINT'(a10,x,5f10.2)', 'Depth',zSurf_SUEWSsurfs(:,i)
722  ! ! PRINT*, '*****************'
723  ! ! END IF
724  !
725  ! END DO
726  END IF
727 
728  ! ---- QF coeffs (was in SUEWS_SAHP.f95, subroutine SAHP_Coefs)
729  baset_hc = -999 ! Initialise QF coeffs
730  qf_a = 0
731  qf_b = 0
732  qf_c = 0
733  ah_min = 0
734  baset_heating = 0
735  baset_cooling = 0
736  ah_slope_heating = 0
737  ah_slope_cooling = 0
738 
739  baset_hc = surfacechar(gridiv, c_baset_hc)
740  qf_a = surfacechar(gridiv, (/c_qf_a1, c_qf_a2/))
741  qf_b = surfacechar(gridiv, (/c_qf_b1, c_qf_b2/))
742  qf_c = surfacechar(gridiv, (/c_qf_c1, c_qf_c2/))
743  ah_min = surfacechar(gridiv, (/c_ahmin_wd, c_ahmin_we/))
744  ah_slope_heating = surfacechar(gridiv, (/c_ahslopeheating_wd, c_ahslopeheating_we/))
745  ah_slope_cooling = surfacechar(gridiv, (/c_ahslopecooling_wd, c_ahslopecooling_we/))
746  baset_heating = surfacechar(gridiv, (/c_tcriticheating_wd, c_tcriticheating_we/))
747  baset_cooling = surfacechar(gridiv, (/c_tcriticcooling_wd, c_tcriticcooling_we/))
748  enprofwd = surfacechar(gridiv, c_enprofwd)
749  enprofwe = surfacechar(gridiv, c_enprofwe)
750  co2mwd = surfacechar(gridiv, c_co2mwd)
751  co2mwe = surfacechar(gridiv, c_co2mwe)
752  traffprofwd = surfacechar(gridiv, c_traffprofwd)
753  traffprofwe = surfacechar(gridiv, c_traffprofwe)
754  popprofwd = surfacechar(gridiv, c_popprofwd)
755  popprofwe = surfacechar(gridiv, c_popprofwe)
756  minqfmetab = surfacechar(gridiv, c_minqfmetab)
757  maxqfmetab = surfacechar(gridiv, c_maxqfmetab)
758  minfcmetab = surfacechar(gridiv, c_minfcmetab)
759  maxfcmetab = surfacechar(gridiv, c_maxfcmetab)
760  frfossilfuel_heat = surfacechar(gridiv, c_frfossilfuel_heat)
761  frfossilfuel_nonheat = surfacechar(gridiv, c_frfossilfuel_nonheat)
762  ef_umolco2perj = surfacechar(gridiv, c_ef_umolco2perj)
763  enef_v_jkm = surfacechar(gridiv, c_enef_v_jkm)
764  fcef_v_kgkm = surfacechar(gridiv, (/c_fcef_v_kgkmwd, c_fcef_v_kgkmwe/))
765  co2pointsource = surfacechar(gridiv, c_co2pointsource)
766  trafficunits = surfacechar(gridiv, c_trafficunits)
767 
768  ! ---- Irrigation
769  ie_start = int(surfacechar(gridiv, c_iestart))
770  ie_end = int(surfacechar(gridiv, c_ieend))
771  internalwateruse_h = surfacechar(gridiv, c_intwu)
772  faut = surfacechar(gridiv, c_faut)
773  h_maintain = surfacechar(gridiv, c_h_maintain) !Automatic irrigation model coefficients [mm d-1]; [mm d-1 degC-1]; [mm d-2]
774  ie_a = surfacechar(gridiv, c_ie_a) !Automatic irrigation model coefficients [mm d-1]; [mm d-1 degC-1]; [mm d-2]
775  ie_m = surfacechar(gridiv, c_ie_m) !Manual irrigation model coefficients [mm d-1]; [mm d-1 degC-1]; [mm d-2]
776  daywat = surfacechar(gridiv, c_daywat)
777  daywatper = surfacechar(gridiv, c_daywatper)
778 
779  ! ---- Hourly profiles
780  ahprof_24hr(0:23, 1) = surfacechar(gridiv, c_hrprofenusewd) ! Anthropogenic heat, weekdays
781  ahprof_24hr(0:23, 2) = surfacechar(gridiv, c_hrprofenusewe) ! Anthropogenic heat, weekends
782  wuprofm_24hr(0:23, 1) = surfacechar(gridiv, c_hrprofwumanuwd) ! Water use, manual, weekdays
783  wuprofm_24hr(0:23, 2) = surfacechar(gridiv, c_hrprofwumanuwe) ! Water use, manual, weekends
784  wuprofa_24hr(0:23, 1) = surfacechar(gridiv, c_hrprofwuautowd) ! Water use, automatic, weekdays
785  wuprofa_24hr(0:23, 2) = surfacechar(gridiv, c_hrprofwuautowe) ! Water use, automatic, weekends
786  snowprof_24hr(0:23, 1) = surfacechar(gridiv, c_hrprofsnowcwd) ! Snow clearing, weekdays
787  snowprof_24hr(0:23, 2) = surfacechar(gridiv, c_hrprofsnowcwe) ! Snow clearing, weekends
788  humactivity_24hr(0:23, 1) = surfacechar(gridiv, c_hrprofhumactivitywd) ! Human activity, weekdays
789  humactivity_24hr(0:23, 2) = surfacechar(gridiv, c_hrprofhumactivitywe) ! Human activity, weekends
790  traffprof_24hr(0:23, 1) = surfacechar(gridiv, c_hrproftraffwd) ! Traffic, weekdays
791  traffprof_24hr(0:23, 2) = surfacechar(gridiv, c_hrproftraffwe) ! Traffic, weekends
792  popprof_24hr(0:23, 1) = surfacechar(gridiv, c_hrprofpopwd) ! Population, weekdays
793  popprof_24hr(0:23, 2) = surfacechar(gridiv, c_hrprofpopwe) ! Population, weekends
794 
795  ! ---- Profiles at the resolution of model time step
796  ! AHProf_tstep(:,1) = TstepProfiles(Gridiv,cTP_EnUseWD,:) ! Anthropogenic heat, weekdays
797  ! AHProf_tstep(:,2) = TstepProfiles(Gridiv,cTP_EnUseWE,:) ! Anthropogenic heat, weekends
798  ! WUProfM_tstep(:,1) = TstepProfiles(Gridiv,cTP_WUManuWD,:) ! Water use, manual, weekdays
799  ! WUProfM_tstep(:,2) = TstepProfiles(Gridiv,cTP_WUManuWE,:) ! Water use, manual, weekends
800  ! WUProfA_tstep(:,1) = TstepProfiles(Gridiv,cTP_WUAutoWD,:) ! Water use, automatic, weekdays
801  ! WUProfA_tstep(:,2) = TstepProfiles(Gridiv,cTP_WUAutoWE,:) ! Water use, automatic, weekends
802  ! HumActivity_tstep(:,1) = TstepProfiles(Gridiv,cTP_HumActivityWD,:) ! Human activity, weekdays
803  ! HumActivity_tstep(:,2) = TstepProfiles(Gridiv,cTP_HumActivityWE,:) ! Human activity, weekends
804  ! TraffProf_tstep(:,1) = TstepProfiles(Gridiv,cTP_TraffProfWD,:) !Traffic, weekdays
805  ! TraffProf_tstep(:,2) = TstepProfiles(Gridiv,cTP_TraffProfWE,:) !Traffic, weekends
806  ! PopProf_tstep(:,1) = TstepProfiles(Gridiv,cTP_PopProfWD,:) !Population, weekdays
807  ! PopProf_tstep(:,2) = TstepProfiles(Gridiv,cTP_PopProfWE,:) !Population, weekends
808 
809  ! ---- Within-grid water distribution
810  ! N.B. Rows and columns of WaterDist are the other way round to the input info
811  !! Model currently does not include above-ground flow from the Water surface
812  !! - Probably should adjust WaterDist to have nsurf columns so that Water can behave like the other surfaces.
813  ! Model returns an error if both ToRunoff and ToSoilStore are non-zero (in CodeMatchDist)
814  ! For impervious surfaces, water goes to runoff; for pervious surfaces, water goes to soilstore
815  waterdist(pavsurf, 1:(nsurf - 1)) = surfacechar(gridiv, c_wgtopaved(1:(nsurf - 1)))
816  waterdist(bldgsurf, 1:(nsurf - 1)) = surfacechar(gridiv, c_wgtobldgs(1:(nsurf - 1)))
817  waterdist(conifsurf, 1:(nsurf - 1)) = surfacechar(gridiv, c_wgtoevetr(1:(nsurf - 1)))
818  waterdist(decidsurf, 1:(nsurf - 1)) = surfacechar(gridiv, c_wgtodectr(1:(nsurf - 1)))
819  waterdist(grasssurf, 1:(nsurf - 1)) = surfacechar(gridiv, c_wgtograss(1:(nsurf - 1)))
820  waterdist(bsoilsurf, 1:(nsurf - 1)) = surfacechar(gridiv, c_wgtobsoil(1:(nsurf - 1)))
821  waterdist(watersurf, 1:(nsurf - 1)) = surfacechar(gridiv, c_wgtowater(1:(nsurf - 1)))
822  ! Runoff or SoilStore row !!Change later to allow both Runoff and SoilStore
823  DO iv = 1, (nsurf - 1)
824  IF (surfacechar(gridiv, c_wgtorunoff(iv)) /= 0) THEN
825  waterdist((nsurf + 1), iv) = surfacechar(gridiv, c_wgtorunoff(iv))
826  ELSE
827  waterdist((nsurf + 1), iv) = surfacechar(gridiv, c_wgtosoilstore(iv))
828  ENDIF
829  ENDDO
830 
831  ! Access required DailyState variables for the current grid (moved HCW 26 Jun 2015)
832  ! HDD(:,:) = HDD_grids(:,:,Gridiv)
833  ! GDD(:,:) = GDD_grids(:,:,Gridiv)
834  ! LAI(:,:) = LAI_grids(:,:,Gridiv)
835  ! WUDay(:,:) = WUDay_grids(:,:,Gridiv)
836  ! AlbDecTr(:) = AlbDecTr_grids(:,Gridiv)
837  ! DecidCap(:) = DecidCap_grids(:,Gridiv)
838  ! Porosity(:) = Porosity_grids(:,Gridiv)
839  ! AlbEveTr(:) = AlbEveTr_grids(:,Gridiv)
840  ! AlbGrass(:) = AlbGrass_grids(:,Gridiv)
841  ! SnowAlb = ModelDailyState(Gridiv, cMDS_SnowAlb)
842 
843  !! ---- Between-grid water distribution
844 !!! Need to make these larger than MaxNumberOfGrids (and recode), as each grid can have 8 connections
845  !!GridConnections(1,) = SurfaceChar(Gridiv,c_Grid)
846  !!GridConnectionsFrac() = SurfaceChar(Gridiv,55)
847  !!GridConnections(2,) = SurfaceChar(Gridiv,54)
848  !
849  !! Fraction of water from each grid
850 !!! N.B. will need to check input files are correctly set up
851  !GridToFrac(1:NConns) = (SurfaceChar(Gridiv,55:69:2))
852  !! Grid where water goes to
853  !GridTo(1:NConns) = (SurfaceChar(Gridiv,54:68:2))
854  !! Come back to this later
855 
856  ! =================================================================================
857 
858  !-----------------------------------------------------
859  !-----------------------------------------------------
860  ! load snow related properties for NARP
861  if (snowuse == 1) narp_emis_snow = surfacechar(gridiv, c_snowemis)
862  !NARP_CONFIGURATION if net radiation is to be modelled
863  IF (netradiationmethod > 0) THEN
864  narp_lat = surfacechar(gridiv, c_lat)
865  narp_long = surfacechar(gridiv, c_lng) ! New sun_position_v2 use degrees FL
866  narp_year = int(surfacechar(gridiv, c_year))
867  narp_tz = timezone !not every 5-min
868  narp_trans_site = trans_site
869  !INTERVAL IS ONLY RELEVANT TO LUPCORR
870  !ALL OTHER CALCULATIONS ARE INTERVAL INDEPENDENT
871  !NB FOR INTERVALS LONGER THAN 15 MINUTES ERRORS IN KCLEAR WILL BE GREATER
872 
873  ! Commented out HCW 04 Mar 2015
874  !NARP_NPERHOUR=MAX(3600/t_INTERVAL,1) !!Check this
875  !IF(ALLOCATED(NARP_KDOWN_HR)) DEALLOCATE(NARP_KDOWN_HR)
876  !ALLOCATE(NARP_KDOWN_HR(NARP_NPERHOUR))
877  !NARP_KDOWN_HR=0.
878 
879  !IF (ldown_option==4.or.ldown_option==5) then !Added by LJ
880  ! INIITIALIZE SMITH DAY OF YEAR GRID G
881  ! NARP_G=SMITHLAMBDA(NINT(LAT))
882  !ENDIF
883  ENDIF
884 
885  !=================================================================================
886  ! When SUEWS_Translate is called from InitialState (ir=0), inputs need translating
887  IF (ir == 0) THEN
888  !write(*,*) 'This should be seen only when called from InitialState and ir is 0. ir:',ir
889 
890  ! =============================================================================
891  ! === Translate inputs from ModelDailyState to variable names used in model ===
892  ! =============================================================================
893 
894  ! Get id_prev from ModelDailyState
895  id_prev = int(modeldailystate(gridiv, cmds_id_prev))
896 
897  snowfallcum = modeldailystate(gridiv, cmds_snowfallcum)
898  snowalb = modeldailystate(gridiv, cmds_snowalb)
899 
900  porosity_id = modeldailystate(gridiv, cmds_porosity)
901  albdectr_id = modeldailystate(gridiv, cmds_albdectr)
902  albevetr_id = modeldailystate(gridiv, cmds_albevetr)
903  albgrass_id = modeldailystate(gridiv, cmds_albgrass)
904  decidcap_id = modeldailystate(gridiv, cmds_decidcap)
905 
906  decidcap_id_grids(gridiv) = decidcap_id
907  albdectr_id_grids(gridiv) = albdectr_id
908  albevetr_id_grids(gridiv) = albevetr_id
909  albgrass_id_grids(gridiv) = albgrass_id
910  porosity_id_grids(gridiv) = porosity_id
911 
912  ! ---- Phenology
913  ! ---- LAI
914  lai_id = 0
915  lai_id(ivconif) = modeldailystate(gridiv, cmds_laiinitialevetr)
916  lai_id(ivdecid) = modeldailystate(gridiv, cmds_laiinitialdectr)
917  lai_id(ivgrass) = modeldailystate(gridiv, cmds_laiinitialgrass)
918 
919  ! GDD_id: GDD Values for one day
920  gdd_id = modeldailystate(gridiv, cmds_gdd1_0)
921  ! SDD_id: SDD Values for one day
922  sdd_id = modeldailystate(gridiv, cmds_gdd2_0)
923  ! Tmin, Tmax: daily minimum and maximum temperatures
924  tmin_id = modeldailystate(gridiv, cmds_gddmin)
925  tmax_id = modeldailystate(gridiv, cmds_gddmax)
926  ! length of daylight
927  lenday_id = 0
928 
929  ! ---- Heating degree days, HDD
930  ! HDD = 0
931  ! HDD(id_prev,1) = ModelDailyState(Gridiv,cMDS_HDD1) ! 1 = Heating
932  ! HDD(id_prev,2) = ModelDailyState(Gridiv,cMDS_HDD2) ! 2 = Cooling
933  ! HDD(id_prev-3,3) = ModelDailyState(Gridiv,cMDS_TempCOld3) ! 3 will become average
934  ! HDD(id_prev-2,3) = ModelDailyState(Gridiv,cMDS_TempCOld2)
935  ! HDD(id_prev-1,3) = ModelDailyState(Gridiv,cMDS_TempCOld1)
936  ! HDD(id_prev,3) = ModelDailyState(Gridiv,cMDS_TempC)
937  ! 4 = 5 day running mean
938  ! 5 = daily precip total
939  ! HDD(id_prev,6) = ModelDailyState(Gridiv,cMDS_DaysSinceRain) ! 6 = days since rain
940 
941  ! ---- Heating degree days, HDD_id: HDD Values for one day
942  hdd_id(1:6) = 0
943 
944  ! HDD_id(1)=ModelDailyState(Gridiv,cMDS_HDD1)
945  ! HDD_id(2)=ModelDailyState(Gridiv,cMDS_HDD2)
946  ! HDD_id(3)=ModelDailyState(Gridiv,cMDS_TempC)
947  ! ! 4 = 5 day running mean
948  ! ! 5 = daily precip total
949  ! HDD_id(6) = ModelDailyState(Gridiv,cMDS_DaysSinceRain)
950 
951  ! Save required DailyState variables for the current grid (HCW 27 Nov 2014)
952  hdd_id_grids(:, gridiv) = hdd_id(:)
953  gdd_id_grids(:, gridiv) = gdd_id(:)
954  sdd_id_grids(:, gridiv) = sdd_id(:)
955  tmin_id_grids(gridiv) = tmin_id
956  tmax_id_grids(gridiv) = tmax_id
957  lenday_id_grids(gridiv) = lenday_id
958  lai_id_grids(:, gridiv) = lai_id(:)
959 
960  ! daily water use
961  wuday_id = 0
962  wuday_id_grids(:, gridiv) = wuday_id(:)
963 
964  ! AlbDecTr_grids(:,Gridiv) = AlbDecTr(:)
965  ! AlbEveTr_grids(:,Gridiv) = AlbEveTr(:)
966  ! AlbGrass_grids(:,Gridiv) = AlbGrass(:)
967  ! DecidCap_grids(:,Gridiv) = DecidCap(:)
968  ! Porosity_grids(:,Gridiv) = Porosity(:)
969 
970  ! ---- Snow density of each surface
971  snowdens(1:nsurf) = modeldailystate(gridiv, cmds_snowdens(1:nsurf))
972 
973  ! =============================================================================
974  ! === Translate inputs from ModelOutputData to variable names used in model ===
975  ! =============================================================================
976  ! ---- Above-ground state
977  state_id(1:nsurf) = modeloutputdata(0, cmod_state(1:nsurf), gridiv)
978  ! stateDay(0,Gridiv,1:nsurf) = ModelOutputData(0,cMOD_State(1:nsurf),Gridiv)
979  ! ---- Below-ground
980  soilstore_id(1:nsurf) = modeloutputdata(0, cmod_soilstate(1:nsurf), gridiv)
981  ! soilmoistDay(0,Gridiv,1:nsurf) = ModelOutputData(0,cMOD_SoilState(1:nsurf),Gridiv)
982  ! ---- Snow fraction
983  snowfrac(1:nsurf) = modeloutputdata(0, cmod_snowfrac(1:nsurf), gridiv)
984  ! ---- Snow water equivalent in SnowPack
985  snowpack(1:nsurf) = modeloutputdata(0, cmod_snowpack(1:nsurf), gridiv)
986  ! ---- Liquid (melted) water in SnowPack
987  snowwater(1:nsurf) = modeloutputdata(0, cmod_snowwaterstate(1:nsurf), gridiv)
988 
989  ENDIF !ir = 0
990  !=================================================================================
991 
992  !=========================== Write FileChoices.txt ===============================
993  !=================================================================================
994  ! Do once per grid per year (was in SUEWS_Initial.f95)
995  IF (ir == 1 .AND. imb == 1) THEN !For first row of first block only
996  !write(*,*) 'Writing to FileChoices for first chunk of met data per year per grid'
997  filechoices = trim(fileoutputpath)//trim(filecode)//'_FileChoices.txt'
998  OPEN (12, file=filechoices, position='append')
999 
1000  WRITE (grid_txt, '(I5)') int(surfacechar(gridiv, c_grid))
1001  WRITE (year_txt, '(I4)') int(surfacechar(gridiv, c_year))
1002  WRITE (ssg_yyyy, '(A12)') trim(filecode)//trim(adjustl(grid_txt))//'_'//trim(adjustl(year_txt))
1003 
1004  !write(12,*) '--------------------------------------------------------------------------------'
1005  WRITE (12, *) '----- '//trim(adjustl(ssg_yyyy))//' Surface characteristics'//' -----'
1006  ! Characteristics that apply to some or all surface types
1007  WRITE (12, '(8a10,a16)') 'Paved', 'Bldgs', 'EveTr', 'DecTr', 'Grass', 'BSoil', 'Water', 'Snow', ' SurfType'
1008  WRITE (12, 120) (sfr(iv), iv=1, nsurf), fcskip, ' SurfFr'
1009  WRITE (12, 120) fcskip, fcskip, irrfracevetr, irrfracdectr, irrfracgrass, fcskip, fcskip, fcskip, ' IrrFr'
1010  WRITE (12, 120) fcskip, fcskip, wuareaevetr_m2, wuareadectr_m2, wuareagrass_m2, fcskip, fcskip, fcskip, ' WaterUseArea'
1011  WRITE (12, 120) fcskip, bldgh, evetreeh, dectreeh, fcskip, fcskip, fcskip, fcskip, ' H'
1012  WRITE (12, 120) fcskip, faibldg, faievetree, faidectree, fcskip, fcskip, fcskip, fcskip, ' FAI'
1013  WRITE (12, 120) fcskip, fcskip, albmin_evetr, albmin_dectr, albmin_grass, fcskip, fcskip, snowalbmin, ' AlbedoMin'
1014  WRITE (12, 120) fcskip, fcskip, albmax_evetr, albmax_dectr, albmax_grass, fcskip, fcskip, snowalbmax, ' AlbedoMax'
1015  !write(12,120) (alb(iv),iv=1,nsurf),SnowAlb, ' Albedo' ! This is instantaneous value (not provided as input)
1016  WRITE (12, 120) (emis(iv), iv=1, nsurf), emis_snow, ' Emissivity'
1017  WRITE (12, 120) fcskip, fcskip, (baset(iv), iv=1, nvegsurf), fcskip, fcskip, fcskip, ' BaseT'
1018  WRITE (12, 120) fcskip, fcskip, (basete(iv), iv=1, nvegsurf), fcskip, fcskip, fcskip, ' BaseTe'
1019  WRITE (12, 120) (storedrainprm(1, iv), iv=1, nsurf), fcskip, ' StorageMin'
1020  WRITE (12, 120) (storedrainprm(5, iv), iv=1, nsurf), fcskip, ' StorageMax'
1021  WRITE (12, 120) (wetthresh(iv), iv=1, nsurf), fcskip, ' WetThreshold'
1022  WRITE (12, 120) (statelimit(iv), iv=1, nsurf), fcskip, ' StateLimit'
1023  WRITE (12, 120) (storedrainprm(2, iv), iv=1, nsurf), fcskip, ' DrainageEq' !real
1024  WRITE (12, 120) (storedrainprm(3, iv), iv=1, nsurf), fcskip, ' DrainageCoef1'
1025  WRITE (12, 120) (storedrainprm(4, iv), iv=1, nsurf), fcskip, ' DrainageCoef2'
1026  WRITE (12, 120) fcskip, fcskip, (gddfull(iv), iv=1, nvegsurf), fcskip, fcskip, fcskip, ' GDDFull'
1027  WRITE (12, 120) fcskip, fcskip, (sddfull(iv), iv=1, nvegsurf), fcskip, fcskip, fcskip, ' SDDFull'
1028  WRITE (12, 120) fcskip, fcskip, (laimin(iv), iv=1, nvegsurf), fcskip, fcskip, fcskip, ' LAIMin'
1029  WRITE (12, 120) fcskip, fcskip, (laimax(iv), iv=1, nvegsurf), fcskip, fcskip, fcskip, ' LAIMax'
1030  WRITE (12, 120) fcskip, fcskip, fcskip, pormin_dec, fcskip, fcskip, fcskip, fcskip, ' PorosityMin'
1031  WRITE (12, 120) fcskip, fcskip, fcskip, pormax_dec, fcskip, fcskip, fcskip, fcskip, ' PorosityMax'
1032  WRITE (12, '(2f10.3,3i10, 3f10.3,a16)') fcskip, fcskip, laitype(1:nvegsurf), fcskip, fcskip, fcskip, ' LAIEq' !integer
1033  WRITE (12, '(2f10.3,3f10.5,3f10.3,a16)') fcskip, fcskip, laipower(1, 1:nvegsurf), fcskip, fcskip, fcskip, ' LAI_LeafGP1'
1034  WRITE (12, '(2f10.3,3f10.5,3f10.3,a16)') fcskip, fcskip, laipower(2, 1:nvegsurf), fcskip, fcskip, fcskip, ' LAI_LeafGP2'
1035  WRITE (12, '(2f10.3,3f10.5,3f10.3,a16)') fcskip, fcskip, laipower(3, 1:nvegsurf), fcskip, fcskip, fcskip, ' LAI_LeafOP1'
1036  WRITE (12, '(2f10.3,3f10.5,3f10.3,a16)') fcskip, fcskip, laipower(4, 1:nvegsurf), fcskip, fcskip, fcskip, ' LAI_LeafOP2'
1037  WRITE (12, 120) fcskip, fcskip, (maxconductance(iv), iv=1, nvegsurf), fcskip, fcskip, fcskip, ' MaxCond'
1038  WRITE (12, 120) (soildepth(iv), iv=1, (nsurf - 1)), fcskip, fcskip, ' SoilDepth'
1039  WRITE (12, 120) (soilstorecap(iv), iv=1, (nsurf - 1)), fcskip, fcskip, ' SoilStoreCap'
1040  WRITE (12, '(6f10.5,2f10.3,a16)') (sathydraulicconduct(iv), iv=1, (nsurf - 1)), fcskip, fcskip, ' SatHydraulicConduct'
1041  ! Not currently coded, but add these later: SoilDensity, InfiltrationRate, OBS_SMDept, OBS_SMCap, OBS_SoilNotRocks
1042  WRITE (12, 120) (snowpacklimit(iv), iv=1, (nsurf - 1)), fcskip, fcskip, ' SnowLimPatch'
1043  WRITE (12, 120) snowlimpaved, snowlimbldg, fcskip, fcskip, fcskip, fcskip, fcskip, fcskip, ' SnowLimRemove'
1044  WRITE (12, 120) (ohm_coef(1:nsurf, 1, 1)), ohm_coef(nsurf + 1, 1, 1), ' OHM_a1_Sum_Wet'
1045  WRITE (12, 120) (ohm_coef(1:nsurf, 2, 1)), ohm_coef(nsurf + 1, 2, 1), ' OHM_a1_Sum_Dry'
1046  WRITE (12, 120) (ohm_coef(1:nsurf, 3, 1)), ohm_coef(nsurf + 1, 3, 1), ' OHM_a1_Win_Wet'
1047  WRITE (12, 120) (ohm_coef(1:nsurf, 4, 1)), ohm_coef(nsurf + 1, 4, 1), ' OHM_a1_Win_Dry'
1048  WRITE (12, 120) (ohm_coef(1:nsurf, 1, 2)), ohm_coef(nsurf + 1, 1, 2), ' OHM_a2_Sum_Wet'
1049  WRITE (12, 120) (ohm_coef(1:nsurf, 2, 2)), ohm_coef(nsurf + 1, 2, 2), ' OHM_a2_Sum_Dry'
1050  WRITE (12, 120) (ohm_coef(1:nsurf, 3, 2)), ohm_coef(nsurf + 1, 3, 2), ' OHM_a2_Win_Wet'
1051  WRITE (12, 120) (ohm_coef(1:nsurf, 4, 2)), ohm_coef(nsurf + 1, 4, 2), ' OHM_a2_Win_Dry'
1052  WRITE (12, 120) (ohm_coef(1:nsurf, 1, 3)), ohm_coef(nsurf + 1, 1, 3), ' OHM_a3_Sum_Wet'
1053  WRITE (12, 120) (ohm_coef(1:nsurf, 2, 3)), ohm_coef(nsurf + 1, 2, 3), ' OHM_a3_Sum_Dry'
1054  WRITE (12, 120) (ohm_coef(1:nsurf, 3, 3)), ohm_coef(nsurf + 1, 3, 3), ' OHM_a3_Win_Wet'
1055  WRITE (12, 120) (ohm_coef(1:nsurf, 4, 3)), ohm_coef(nsurf + 1, 4, 3), ' OHM_a3_Win_Dry'
1056  WRITE (12, 120) (ohm_threshsw(1:nsurf)), ohm_threshsw(nsurf + 1), ' OHMthreshold_SW'
1057  WRITE (12, 120) (ohm_threshwd(1:nsurf)), ohm_threshwd(nsurf + 1), ' OHMthreshold_WD'
1058 
1059  WRITE (12, *) '----- '//trim(adjustl(ssg_yyyy))//' Snow parameters'//' -----'
1060  WRITE (12, '(a12,11a10)') 'Grid', 'RadMeltF', 'TempMeltF', 'tau_a', 'tau_f', 'PLimAlb', 'SDensMin', 'SDensMax', &
1061  'tau_r', 'CRWMin', 'CRWMax', 'PLimSnow'
1062  WRITE (12, '(a12,11f10.4)') ssg_yyyy, radmeltfact, tempmeltfact, tau_a, tau_f, preciplimitalb, snowdensmin, snowdensmax, &
1063  tau_r, crwmin, crwmax, preciplimit
1064 
1065  WRITE (12, *) '----- '//trim(adjustl(ssg_yyyy))//' Conductance parameters'//' -----'
1066  WRITE (12, '(a12,12a10)') 'Grid', 'G1', 'G2', 'G3', 'G4', 'G5', 'G6', 'TH', 'TL', 'S1', 'S2', 'Kmax', 'gsModel'
1067  WRITE (12, '(a12,11f10.3,i3)') ssg_yyyy, g1, g2, g3, g4, g5, g6, th, tl, s1, s2, kmax, gsmodel
1068 
1069  WRITE (12, *) '----- '//trim(adjustl(ssg_yyyy))//' Energy-use parameters'//' -----'
1070  WRITE (12, '(a12,11a10)') 'Grid', 'PopDensDaytime', 'BaseT_HC', 'QF_A_WD', 'QF_A_WE', 'QF_B_WD', 'QF_B_WE', 'QF_C_WD', &
1071  'QF_C_WE', 'AH_Min', 'AH_Slope', 'BaseT_Heating'
1072  WRITE (12, '(a12,11f10.3)') ssg_yyyy, popdensdaytime, baset_hc, qf_a(1:2), qf_b(1:2), qf_c(1:2), &
1073  ah_min, ah_slope_heating, baset_heating
1074 
1075  WRITE (12, *) '----- '//trim(adjustl(ssg_yyyy))//' Water-use parameters'//' -----'
1076  WRITE (12, '(a12,10a10)') 'Grid', 'IeStart', 'IeEnd', 'IntWatUse', 'Faut', &
1077  'Ie_a1', 'Ie_a2', 'Ie_a3', 'Ie_m1', 'Ie_m2', 'Ie_m3'
1078  WRITE (12, '(a12,2i10,8f10.3)') ssg_yyyy, ie_start, ie_end, internalwateruse_h, faut, &
1079  ie_a(1:3), ie_m(1:3)
1080 
1081  WRITE (12, *) '----- '//trim(adjustl(ssg_yyyy))//' Weekly profiles'//' -----'
1082  WRITE (12, '(a12,7a10, a16)') 'Grid', '1_Sun', '2_Mon', '3_Tue', '4_Wed', '5_Thu', '6_Fri', '7_Sat', ' DayOfWeek'
1083  WRITE (12, '(a12,7f10.3,a16)') ssg_yyyy, daywat(1:7), ' Irr allowed'
1084  WRITE (12, '(a12,7f10.3,a16)') ssg_yyyy, daywatper(1:7), ' Frac properties'
1085 
1086  WRITE (12, *) '----- '//trim(adjustl(ssg_yyyy))//' Hourly profiles'//' -----'
1087  WRITE (12, '(a12,24i10,a20)') 'Grid', (iv, iv=0, 23), 'HourOfDay'
1088  WRITE (12, 121) ssg_yyyy, ahprof_24hr(0:23, 1), ' Anthrop heat WD'
1089  WRITE (12, 121) ssg_yyyy, ahprof_24hr(0:23, 2), ' Anthrop heat WE'
1090  WRITE (12, 121) ssg_yyyy, wuprofm_24hr(0:23, 1), ' Manual water use WD'
1091  WRITE (12, 121) ssg_yyyy, wuprofm_24hr(0:23, 2), ' Manual water use WE'
1092  WRITE (12, 121) ssg_yyyy, wuprofa_24hr(0:23, 1), ' Auto. water use WD'
1093  WRITE (12, 121) ssg_yyyy, wuprofa_24hr(0:23, 2), ' Auto. water use WE'
1094  WRITE (12, 121) ssg_yyyy, snowprof_24hr(0:23, 1), ' Snow clearing WD'
1095  WRITE (12, 121) ssg_yyyy, snowprof_24hr(0:23, 2), ' Snow clearing WE'
1096 
1097  WRITE (12, *) '----- '//trim(adjustl(ssg_yyyy))//' Within-grid water distribution'//' -----'
1098  WRITE (12, '(9a10)') 'ToPaved', 'ToBldgs', 'ToEveTr', 'ToDecTr', 'ToGrass', 'ToBSoil', 'ToWater', 'ToROorSS'
1099 
1100  DO iv = 1, (nsurf - 1)
1101  WRITE (12, '(8f10.4)') (waterdist(j, iv), j=1, nsurf + 1)
1102  ENDDO
1103 
1104  WRITE (12, *) '----- '//trim(adjustl(ssg_yyyy))//' Other parameters'//' -----'
1105  WRITE (12, '(a12,7a10)') 'Grid', 'FlowChange', 'ROToWater', 'PipeCap', & ! Water-related
1106  'DrRate', 'Cover', 'MaxRes', & ! LUMPS-related
1107  'Trans' ! NARP-related
1108  WRITE (12, '(a12,7f10.3)') ssg_yyyy, flowchange, runofftowater, pipecapacity, &
1109  drainrt, raincover, rainmaxres, &
1110  trans_site
1111 
1112  WRITE (12, *) '----- '//trim(adjustl(ssg_yyyy))//' Site parameters'//' -----'
1113  WRITE (12, '(a12,9a10)') &
1114  'Grid', 'lat', 'lon', 'tz', 'alt', 'SurfA_ha', 'z', 'PopDensNighttime', 'z0_input', 'zd_input', 'StartDLS', 'EndDLS'
1115  WRITE (12, '(a12,4f10.4,f10.2,4f10.4,2i10)') &
1116  ssg_yyyy, lat, lng*(-1.0), timezone, alt, surfacearea_ha, z, popdensnighttime, z0m_in, zdm_in, &
1117  startdls, enddls! DayLightSavingDay(1:2)
1118 
1119  WRITE (12, *) ''
1120 
1121  CLOSE (12)
1122 
1123  !==============================================================================
1124  ! Check input values are reasonable ===========================================
1125 
1126  ! Coefficients for anthropogenic heat models ----------------------------------
1127  IF (emissionsmethod == 1) THEN !Loridan et al. (2011) calculation
1128  IF (ah_min(1) == 0 .AND. ah_slope_heating(1) == 0 .AND. baset_heating(1) == 0) THEN
1129  CALL errorhint(53, 'Check QF calculation coefficients.', notused, notused, emissionsmethod)
1130  ENDIF
1131 
1132  ELSEIF (emissionsmethod == 2) THEN !Jarvi et al. (2011) calculation
1133  IF (sum(qf_a) == 0 .AND. sum(qf_b) == 0 .AND. sum(qf_c) == 0) THEN
1134  CALL errorhint(54, 'Check QF calculation coefficients.', notused, notused, emissionsmethod)
1135  ENDIF
1136  ENDIF
1137 
1138  ! Morphometric parameters -----------------------------------------------------
1139  IF (roughlenmommethod == 1) THEN !z0, zd values provided in input file
1140  ! Check z0m and zd are reasonable
1141  IF (z0m < 0.00001) CALL errorhint(1, 'z0 value provided is very small (RoughLenMomMethod=1).', z0m, notused, gridid)
1142  IF (zdm < 0.00001) CALL errorhint(1, 'zd value provided is very small (RoughLenMomMethod=1).', zdm, notused, gridid)
1143  zzd = z - zdm
1144  ELSEIF (roughlenmommethod == 3) THEN !z0, zd calculated using FAI provided in input file
1145  ! Check FAIs reasonable
1146  IF (faibldg < 0) CALL errorhint(1, &
1147  'FAI_Bldgs value provided is very small (RoughLenMomMethod=3)', &
1148  faibldg, notused, gridid)
1149  IF (faitree < 0) CALL errorhint(1, &
1150  'FAI_EveTr/DecTr value provided is very small (RoughLenMomMethod=3)', &
1151  faitree, notused, gridid)
1152  ENDIF
1153 
1154  ENDIF !End for first row of first block only ===================================
1155 
1156  !=================================================================================
1157  !For each row of the met forcing file (ir), translate correct info for each grid
1158  ! into model variables
1159  IF (ir > 0) THEN
1160  ! =============================================================================
1161  ! === Translate met data from MetForcingData to variable names used in model ==
1162  ! =============================================================================
1163  iy = int(metforcingdata(ir, 1, gridiv)) !Integer variables
1164  id = int(metforcingdata(ir, 2, gridiv))
1165  it = int(metforcingdata(ir, 3, gridiv))
1166  imin = int(metforcingdata(ir, 4, gridiv))
1167  isec = 0 ! NOT used by SUEWS but by WRF-SUEWS via the cal_main interface
1168  qn1_obs = metforcingdata(ir, 5, gridiv) !Real values (kind(1d0))
1169  qh_obs = metforcingdata(ir, 6, gridiv)
1170  qe_obs = metforcingdata(ir, 7, gridiv)
1171  qs_obs = metforcingdata(ir, 8, gridiv)
1172  qf_obs = metforcingdata(ir, 9, gridiv)
1173  avu1 = metforcingdata(ir, 10, gridiv)
1174  avrh = metforcingdata(ir, 11, gridiv)
1175  temp_c = metforcingdata(ir, 12, gridiv)
1176  press_hpa = metforcingdata(ir, 13, gridiv)
1177  precip = metforcingdata(ir, 14, gridiv)
1178  avkdn = metforcingdata(ir, 15, gridiv)
1179  snowfrac_obs = metforcingdata(ir, 16, gridiv)
1180  ldown_obs = metforcingdata(ir, 17, gridiv)
1181  fcld_obs = metforcingdata(ir, 18, gridiv)
1182  wu_m3 = metforcingdata(ir, 19, gridiv)
1183  xsmd = metforcingdata(ir, 20, gridiv)
1184  lai_obs = metforcingdata(ir, 21, gridiv)
1185  kdiff = metforcingdata(ir, 22, gridiv)
1186  kdir = metforcingdata(ir, 23, gridiv)
1187  wdir = metforcingdata(ir, 24, gridiv)
1188 
1189  ! get qn1 memory for previous time steps
1190  dqndt = dqndt_grids(gridiv)
1191  qn1_av = qn1_av_grids(gridiv)
1192  tair_av = tair_av_grids(gridiv)
1193 
1194  IF (snowuse == 1) THEN
1195  dqnsdt = dqnsdt_grids(gridiv)
1196  qn1_s_av = qn1_s_av_grids(gridiv)
1197  ENDIF
1198 
1199  ! added by TS 29 Jun 2018 to remove annual loops in main calculation
1200  gdd_id = gdd_id_grids(:, gridiv)
1201  sdd_id = sdd_id_grids(:, gridiv)
1202  tmin_id = tmin_id_grids(gridiv)
1203  tmax_id = tmax_id_grids(gridiv)
1204  lenday_id = lenday_id_grids(gridiv)
1205  hdd_id = hdd_id_grids(:, gridiv)
1206  lai_id = lai_id_grids(:, gridiv)
1207  wuday_id = wuday_id_grids(:, gridiv)
1208 
1209  decidcap_id = decidcap_id_grids(gridiv)
1210  albdectr_id = albdectr_id_grids(gridiv)
1211  albevetr_id = albevetr_id_grids(gridiv)
1212  albgrass_id = albgrass_id_grids(gridiv)
1213  porosity_id = porosity_id_grids(gridiv)
1214 
1215  ! get met array for one grid used in AnOHM
1216  metforcingdata_grid = metforcingdata(:, :, gridiv)
1217 
1218  ! Calculate dectime
1219  dectime = REAL(id - 1, KIND(1d0)) + REAL(it, kind(1d0))/24 + REAL(imin, kind(1d0))/(60*24)
1220  ! Create datetime stamp for error/warnings file
1221  WRITE (iy_text, '(i4)') iy
1222  WRITE (id_text, '(i3)') id
1223  WRITE (it_text, '(i2)') it
1224  WRITE (imin_text, '(i2)') imin
1225  ! datetime = TRIM(ADJUSTL(iy_text))//' '//TRIM(ADJUSTL(id_text))//' '//TRIM(ADJUSTL(it_text))//' '//TRIM(ADJUSTL(imin_text))
1226  WRITE (gridid_text, '(i10)') gridid
1227 
1228  ! =============================================================================
1229  ! === Translate values from ModelDailyState to variable names used in model ===
1230  ! =============================================================================
1231  ! porosity(id) = ModelDailyState(Gridiv,cMDS_porosity)
1232  ! albDecTr(id) = ModelDailyState(Gridiv,cMDS_albDecTr)
1233  ! albEveTr(id) = ModelDailyState(Gridiv,cMDS_albEveTr)
1234  ! albGrass(id) = ModelDailyState(Gridiv,cMDS_albGrass)
1235  ! DecidCap(id) = ModelDailyState(Gridiv,cMDS_DecidCap)
1236 
1237  ! SnowfallCum is instantaneous values and should be translated at each tstep, TS 17 Sep 2019
1238  snowfallcum = modeldailystate(gridiv, cmds_snowfallcum)
1239  ! ---- Snow density of each surface
1240  snowdens(1:nsurf) = modeldailystate(gridiv, cmds_snowdens(1:nsurf))
1241  ! ---- Snow albedo
1242  snowalb = modeldailystate(gridiv, cmds_snowalb)
1243 
1244  ! =============================================================================
1245  ! === Translate values from ModelOutputData to variable names used in model ===
1246  ! =============================================================================
1247  ! ---- Above-ground state
1248  state_id(1:nsurf) = modeloutputdata(ir - 1, cmod_state(1:nsurf), gridiv)
1249  ! ---- Below-ground state
1250  soilstore_id(1:nsurf) = modeloutputdata(ir - 1, cmod_soilstate(1:nsurf), gridiv)
1251  ! ---- Snow fraction
1252  snowfrac(1:nsurf) = modeloutputdata(ir - 1, cmod_snowfrac(1:nsurf), gridiv)
1253  ! ---- Snow water equivalent in SnowPack
1254  snowpack(1:nsurf) = modeloutputdata(ir - 1, cmod_snowpack(1:nsurf), gridiv)
1255  ! ---- Liquid (melted) water in SnowPack
1256  snowwater(1:nsurf) = modeloutputdata(ir - 1, cmod_snowwaterstate(1:nsurf), gridiv)
1257 
1258  ! ---- ice fraction
1259  icefrac = icefrac_grids(:, gridiv)
1260 
1261  !Also translate ESTM forcing data
1262  IF (storageheatmethod == 4 .OR. storageheatmethod == 14) THEN
1263  ! write(*,*) 'Translating ESTM forcing data'
1264  ts5mindata(ir, 1:ncolsestmdata) = estmforcingdata(ir, 1:ncolsestmdata, gridiv)
1265  ts5mindata_ir(1:ncolsestmdata) = estmforcingdata(ir, 1:ncolsestmdata, gridiv)
1266  CALL estm_translate(gridiv)
1267  ENDIF
1268 
1269  ENDIF !ir>0 !===================================================================
1270 
1271  ! --------------------------------------------------------------------------------
1272  ! Check Initial Conditions are reasonable ----------------------------------------
1273  IF (ir == 1 .AND. imb == 1) THEN !For first row of first block only
1274  CALL checkinitial
1275  ENDIF
1276  ! --------------------------------------------------------------------------------
1277 
1278  ! ======================================================================
1279  ! write out initial conditions for debugging supy
1280  if (ir == 1 .and. imb == 1) then
1281  filestateinit = trim(fileoutputpath)//trim(filecode)//trim(adjustl(grid_txt))//'_'//trim(adjustl(year_txt))//'_state_init.txt'
1282  OPEN (12, file=filestateinit, position='rewind')
1283 
1284  write (12, *) '&state_init'
1285  write (12, *) 'aerodynamicresistancemethod=', aerodynamicresistancemethod
1286  write (12, *) 'ah_min=', ah_min
1287  write (12, *) 'ahprof_24hr=', ahprof_24hr
1288  write (12, *) 'ah_slope_cooling=', ah_slope_cooling
1289  write (12, *) 'ah_slope_heating=', ah_slope_heating
1290  write (12, *) 'alb=', alb
1291  write (12, *) 'albmax_dectr=', albmax_dectr
1292  write (12, *) 'albmax_evetr=', albmax_evetr
1293  write (12, *) 'albmax_grass=', albmax_grass
1294  write (12, *) 'albmin_dectr=', albmin_dectr
1295  write (12, *) 'albmin_evetr=', albmin_evetr
1296  write (12, *) 'albmin_grass=', albmin_grass
1297  write (12, *) 'alpha_bioco2=', alpha_bioco2
1298  write (12, *) 'alpha_enh_bioco2=', alpha_enh_bioco2
1299  write (12, *) 'alt=', alt
1300  write (12, *) 'avkdn=', avkdn
1301  write (12, *) 'avrh=', avrh
1302  write (12, *) 'avu1=', avu1
1303  write (12, *) 'baset=', baset
1304  write (12, *) 'basete=', basete
1305  write (12, *) 'BaseT_HC=', baset_hc
1306  write (12, *) 'beta_bioco2=', beta_bioco2
1307  write (12, *) 'beta_enh_bioco2=', beta_enh_bioco2
1308  write (12, *) 'bldgh=', bldgh
1309  write (12, *) 'capmax_dec=', capmax_dec
1310  write (12, *) 'capmin_dec=', capmin_dec
1311  write (12, *) 'chanohm=', chanohm
1312  write (12, *) 'co2pointsource=', co2pointsource
1313  write (12, *) 'cpanohm=', cpanohm
1314  write (12, *) 'crwmax=', crwmax
1315  write (12, *) 'crwmin=', crwmin
1316  write (12, *) 'daywat=', daywat
1317  write (12, *) 'daywatper=', daywatper
1318  write (12, *) 'dectreeh=', dectreeh
1319  write (12, *) 'diagnose=', diagnose
1320  write (12, *) 'diagqn=', diagqn
1321  write (12, *) 'diagqs=', diagqs
1322  write (12, *) 'drainrt=', drainrt
1323  write (12, *) 'dt_since_start=', dt_since_start
1324  write (12, *) 'dqndt=', dqndt
1325  write (12, *) 'qn1_av=', qn1_av
1326  write (12, *) 'dqnsdt=', dqnsdt
1327  write (12, *) 'qn1_s_av=', qn1_s_av
1328  write (12, *) 'ef_umolco2perj=', ef_umolco2perj
1329  write (12, *) 'emis=', emis
1330  write (12, *) 'emissionsmethod=', emissionsmethod
1331  write (12, *) 'enef_v_jkm=', enef_v_jkm
1332  write (12, *) 'enddls=', enddls
1333  write (12, *) 'evetreeh=', evetreeh
1334  write (12, *) 'faibldg=', faibldg
1335  write (12, *) 'faidectree=', faidectree
1336  write (12, *) 'faievetree=', faievetree
1337  write (12, *) 'faut=', faut
1338  write (12, *) 'fcef_v_kgkm=', fcef_v_kgkm
1339  write (12, *) 'fcld_obs=', fcld_obs
1340  write (12, *) 'flowchange=', flowchange
1341  write (12, *) 'frfossilfuel_heat=', frfossilfuel_heat
1342  write (12, *) 'frfossilfuel_nonheat=', frfossilfuel_nonheat
1343  write (12, *) 'g1=', g1
1344  write (12, *) 'g2=', g2
1345  write (12, *) 'g3=', g3
1346  write (12, *) 'g4=', g4
1347  write (12, *) 'g5=', g5
1348  write (12, *) 'g6=', g6
1349  write (12, *) 'gdd_id=', gdd_id
1350  write (12, *) 'gddfull=', gddfull
1351  write (12, *) 'gridiv=', gridiv
1352  write (12, *) 'gsmodel=', gsmodel
1353  write (12, *) 'hdd_id=', hdd_id
1354  write (12, *) 'humactivity_24hr=', humactivity_24hr
1355  write (12, *) 'icefrac=', icefrac
1356  write (12, *) 'id=', id
1357  write (12, *) 'ie_a=', ie_a
1358  write (12, *) 'ie_end=', ie_end
1359  write (12, *) 'ie_m=', ie_m
1360  write (12, *) 'ie_start=', ie_start
1361  write (12, *) 'imin=', imin
1362  write (12, *) 'internalwateruse_h=', internalwateruse_h
1363  write (12, *) 'IrrFracEveTr=', irrfracevetr
1364  write (12, *) 'IrrFracDecTr=', irrfracdectr
1365  write (12, *) 'irrfracgrass=', irrfracgrass
1366  write (12, *) 'isec=', isec
1367  write (12, *) 'it=', it
1368  write (12, *) 'evapmethod=', evapmethod
1369  write (12, *) 'iy=', iy
1370  write (12, *) 'kkanohm=', kkanohm
1371  write (12, *) 'kmax=', kmax
1372  write (12, *) 'lai_id=', lai_id
1373  write (12, *) 'laicalcyes=', laicalcyes
1374  write (12, *) 'laimax=', laimax
1375  write (12, *) 'laimin=', laimin
1376  write (12, *) 'lai_obs=', lai_obs
1377  write (12, *) 'laipower=', laipower
1378  write (12, *) 'laitype=', laitype
1379  write (12, *) 'lat=', lat
1380  write (12, *) 'lenday_id=', lenday_id
1381  write (12, *) 'ldown_obs=', ldown_obs
1382  write (12, *) 'lng=', lng
1383  write (12, *) 'maxconductance=', maxconductance
1384  write (12, *) 'maxfcmetab=', maxfcmetab
1385  write (12, *) 'maxqfmetab=', maxqfmetab
1386  write (12, *) 'snowwater=', snowwater
1387  ! write (12, *) 'metforcingdata_grid=', metforcingdata_grid
1388  write (12, *) 'minfcmetab=', minfcmetab
1389  write (12, *) 'minqfmetab=', minqfmetab
1390  write (12, *) 'min_res_bioco2=', min_res_bioco2
1391  write (12, *) 'narp_emis_snow=', narp_emis_snow
1392  write (12, *) 'narp_trans_site=', narp_trans_site
1393  write (12, *) 'netradiationmethod=', netradiationmethod
1394  write (12, *) 'ohm_coef=', ohm_coef
1395  write (12, *) 'ohmincqf=', ohmincqf
1396  write (12, *) 'ohm_threshsw=', ohm_threshsw
1397  write (12, *) 'ohm_threshwd=', ohm_threshwd
1398  write (12, *) 'pipecapacity=', pipecapacity
1399  write (12, *) 'popdensdaytime=', popdensdaytime
1400  write (12, *) 'popdensnighttime=', popdensnighttime
1401  write (12, *) 'popprof_24hr=', popprof_24hr
1402  write (12, *) 'pormax_dec=', pormax_dec
1403  write (12, *) 'pormin_dec=', pormin_dec
1404  write (12, *) 'precip=', precip
1405  write (12, *) 'preciplimit=', preciplimit
1406  write (12, *) 'preciplimitalb=', preciplimitalb
1407  write (12, *) 'press_hpa=', press_hpa
1408  write (12, *) 'qf0_beu=', qf0_beu
1409  write (12, *) 'qf_a=', qf_a
1410  write (12, *) 'qf_b=', qf_b
1411  write (12, *) 'qf_c=', qf_c
1412  write (12, *) 'qn1_obs=', qn1_obs
1413  write (12, *) 'qh_obs=', qh_obs
1414  write (12, *) 'qs_obs=', qs_obs
1415  write (12, *) 'qf_obs=', qf_obs
1416  write (12, *) 'radmeltfact=', radmeltfact
1417  write (12, *) 'raincover=', raincover
1418  write (12, *) 'rainmaxres=', rainmaxres
1419  write (12, *) 'resp_a=', resp_a
1420  write (12, *) 'resp_b=', resp_b
1421  write (12, *) 'roughlenheatmethod=', roughlenheatmethod
1422  write (12, *) 'roughlenmommethod=', roughlenmommethod
1423  write (12, *) 'runofftowater=', runofftowater
1424  write (12, *) 's1=', s1
1425  write (12, *) 's2=', s2
1426  write (12, *) 'sathydraulicconduct=', sathydraulicconduct
1427  write (12, *) 'sddfull=', sddfull
1428  write (12, *) 'sdd_id=', sdd_id
1429  write (12, *) 'sfr=', sfr
1430  write (12, *) 'smdmethod=', smdmethod
1431  write (12, *) 'snowalb=', snowalb
1432  write (12, *) 'snowalbmax=', snowalbmax
1433  write (12, *) 'snowalbmin=', snowalbmin
1434  write (12, *) 'snowpacklimit=', snowpacklimit
1435  write (12, *) 'snowdens=', snowdens
1436  write (12, *) 'snowdensmax=', snowdensmax
1437  write (12, *) 'snowdensmin=', snowdensmin
1438  write (12, *) 'snowfallcum=', snowfallcum
1439  write (12, *) 'snowfrac=', snowfrac
1440  write (12, *) 'snowlimbldg=', snowlimbldg
1441  write (12, *) 'snowlimpaved=', snowlimpaved
1442  write (12, *) 'snowfrac_obs=', snowfrac_obs
1443  write (12, *) 'snowpack=', snowpack
1444  write (12, *) 'snowprof_24hr=', snowprof_24hr
1445  write (12, *) 'snowuse=', snowuse
1446  write (12, *) 'soildepth=', soildepth
1447  write (12, *) 'soilstore_id=', soilstore_id
1448  write (12, *) 'soilstorecap=', soilstorecap
1449  write (12, *) 'stabilitymethod=', stabilitymethod
1450  write (12, *) 'startdls=', startdls
1451  write (12, *) 'state_id=', state_id
1452  write (12, *) 'statelimit=', statelimit
1453  write (12, *) 'storageheatmethod=', storageheatmethod
1454  write (12, *) 'storedrainprm=', storedrainprm
1455  write (12, *) 'surfacearea=', surfacearea
1456  write (12, *) 'tair_av=', tair_av
1457  write (12, *) 'tau_a=', tau_a
1458  write (12, *) 'tau_f=', tau_f
1459  write (12, *) 'tau_r=', tau_r
1460  write (12, *) 'tmax_id=', tmax_id
1461  write (12, *) 'tmin_id=', tmin_id
1462  write (12, *) 'BaseT_Cooling=', baset_cooling
1463  write (12, *) 'BaseT_Heating=', baset_heating
1464  write (12, *) 'temp_c=', temp_c
1465  write (12, *) 'tempmeltfact=', tempmeltfact
1466  write (12, *) 'th=', th
1467  write (12, *) 'theta_bioco2=', theta_bioco2
1468  write (12, *) 'timezone=', timezone
1469  write (12, *) 'tl=', tl
1470  write (12, *) 'trafficrate=', trafficrate
1471  write (12, *) 'trafficunits=', trafficunits
1472  write (12, *) 'traffprof_24hr=', traffprof_24hr
1473  ! write (12, *) 'ts5mindata_ir=', ts5mindata_ir
1474  write (12, *) 'tstep=', tstep
1475  write (12, *) 'tstep_prev=', tstep_prev
1476  write (12, *) 'veg_type=', veg_type
1477  write (12, *) 'waterdist=', waterdist
1478  write (12, *) 'waterusemethod=', waterusemethod
1479  write (12, *) 'wetthresh=', wetthresh
1480  write (12, *) 'wu_m3=', wu_m3
1481  write (12, *) 'wuday_id=', wuday_id
1482  write (12, *) 'decidcap_id=', decidcap_id
1483  write (12, *) 'albdectr_id=', albdectr_id
1484  write (12, *) 'albevetr_id=', albevetr_id
1485  write (12, *) 'albgrass_id=', albgrass_id
1486  write (12, *) 'porosity_id=', porosity_id
1487  write (12, *) 'wuprofa_24hr=', wuprofa_24hr
1488  write (12, *) 'wuprofm_24hr=', wuprofm_24hr
1489  write (12, *) 'xsmd=', xsmd
1490  write (12, *) 'z=', z
1491  write (12, *) 'z0m_in=', z0m_in
1492  write (12, *) 'zdm_in=', zdm_in
1493  write (12, *) '/'
1494 
1495  WRITE (12, *) ''
1496 
1497  CLOSE (12)
1498 
1499  end if
1500 
1501  ! ======================================================================
1502 
1503  RETURN
1504 
1505 120 FORMAT(8f10.3, a16) !format (10g10.2)
1506 121 FORMAT(a12, 24f10.4, a20)
1507 
1508 END SUBROUTINE suews_translate
1509 !===================================================================================
1510 
1511 !SUEWS_TranslateBack
1512 !Translates model variables to arrays for each grid
1513 !Runs at the end of SUEWS_Calculations to store correct info for each grid
1514 !Made by HW Nov 2014
1515 !-----------------------------------------------------------------------------------
1516 !Last modified:LJ 14 Sep 2015
1517 ! HCW 28 Nov 2014
1518 !===================================================================================
1519 SUBROUTINE suews_translateback(Gridiv, ir, irMax)
1521  USE allocatearray
1522  USE colnamesinputfiles
1524  USE data_in
1525  USE defaultnotused
1526  USE gis_data
1527  USE initial
1528  USE mod_z
1529  USE resist
1530  USE snowmod
1531  USE sues_data
1532  USE time
1533 
1534  IMPLICIT NONE
1535 
1536  INTEGER::Gridiv, & ! Index of the analysed grid (Gridcounter)
1537  ir, & ! Meteorological forcing file index (set to zero if SUEWS_Translate called from InitialState)
1538  irMax ! Last row in current chunk of met data
1539 
1540  ! =============================================================================
1541  ! === Translate values from variable names used in model to ModelDailyState ===
1542  ! =============================================================================
1543 
1544  ! ModelDailyState(Gridiv,cMDS_porosity) = porosity(id)
1545  ! ModelDailyState(Gridiv,cMDS_albDecTr) = albDecTr(id)
1546  ! ModelDailyState(Gridiv,cMDS_albEveTr) = albEveTr(id)
1547  ! ModelDailyState(Gridiv,cMDS_albGrass) = albGrass(id)
1548  ! ModelDailyState(Gridiv,cMDS_DecidCap) = DecidCap(id)
1550 
1556 
1557  ! Save required DailyState variables for the current grid (HCW 27 Nov 2014)
1558  ! HDD_grids(:,:,Gridiv) = HDD(:,:)
1559  ! GDD_grids(:,:,Gridiv) = GDD(:,:)
1560  ! LAI_grids(:,:,Gridiv) = LAI(:,:)
1561  ! WUDay_grids(:,:,Gridiv) = WUDay(:,:)
1562  ! AlbDecTr_grids(:,Gridiv) = AlbDecTr(:)
1563  ! AlbEveTr_grids(:,Gridiv) = AlbEveTr(:)
1564  ! AlbGrass_grids(:,Gridiv) = AlbGrass(:)
1565  ! DecidCap_grids(:,Gridiv) = DecidCap(:)
1566  ! Porosity_grids(:,Gridiv) = Porosity(:)
1567 
1568  ! ! update qn1 memory with values of current time step
1569  ! qn1_store(:,Gridiv) = qn1_store_grid(:)
1570  ! qn1_av_store(:,Gridiv) = qn1_av_store_grid(:)
1571  ! IF (SnowUse == 1) THEN
1572  ! qn1_S_store(:,Gridiv) = qn1_S_store_grid(:)
1573  ! qn1_S_av_store(:,Gridiv) = qn1_S_av_store_grid(:)
1574  ! ENDIF
1575 
1576  ! update averaged qn1 memory
1577  dqndt_grids(gridiv) = dqndt
1578  qn1_av_grids(gridiv) = qn1_av
1579  tair_av_grids(gridiv) = tair_av
1580  IF (snowuse == 1) THEN
1581  dqnsdt_grids(gridiv) = dqnsdt
1582  qn1_s_av_grids(gridiv) = qn1_s_av
1583  ENDIF
1584 
1585  ! added by TS 29 Jun 2018 to remove annual loops in main calculation
1586  gdd_id_grids(:, gridiv) = gdd_id
1587  sdd_id_grids(:, gridiv) = sdd_id
1588  tmin_id_grids(gridiv) = tmin_id
1589  tmax_id_grids(gridiv) = tmax_id
1590  lenday_id_grids(gridiv) = lenday_id
1591  hdd_id_grids(:, gridiv) = hdd_id
1592  lai_id_grids(:, gridiv) = lai_id
1593  wuday_id_grids(:, gridiv) = wuday_id
1594 
1595  decidcap_id_grids(gridiv) = decidcap_id
1596  albdectr_id_grids(gridiv) = albdectr_id
1597  albevetr_id_grids(gridiv) = albevetr_id
1598  albgrass_id_grids(gridiv) = albgrass_id
1599  porosity_id_grids(gridiv) = porosity_id
1600 
1601  ! ---- Snow density of each surface
1604 
1605  ! =============================================================================
1606  ! === Translate values from variable names used in model to ModelOutputData ===
1607  ! =============================================================================
1608 
1609  modeloutputdata(ir, cmod_state(1:nsurf), gridiv) = state_id(1:nsurf)
1611  modeloutputdata(ir, cmod_snowfrac(1:nsurf), gridiv) = snowfrac(1:nsurf)
1612  modeloutputdata(ir, cmod_snowpack(1:nsurf), gridiv) = snowpack(1:nsurf)
1614 
1615  ! ---- ice fraction
1616  icefrac_grids(:, gridiv) = icefrac
1617 
1618  IF (ir == irmax) THEN !Store variables ready for next chunk of met data
1619  modeloutputdata(0, cmod_state(1:nsurf), gridiv) = state_id(1:nsurf)
1621  modeloutputdata(0, cmod_snowfrac(1:nsurf), gridiv) = snowfrac(1:nsurf)
1622  modeloutputdata(0, cmod_snowpack(1:nsurf), gridiv) = snowpack(1:nsurf)
1624  ENDIF
1625 
1626  RETURN
1627 endsubroutine suews_translateback
1628 !===================================================================================
real(kind(1d0)) areazh
integer, dimension(nvegsurf) c_theta_bioco2
integer, dimension(nvegsurf) c_resp_b
real(kind(1d0)) evetreeh
real(kind(1d0)) qn1_s_av
real(kind(1d0)) impervfraction
real(kind(1d0)) tau_r
real(kind(1d0)), dimension(nsurf) snowdens
integer, dimension(nvegsurf) c_leafgp2
real(kind(1d0)) g1
integer, dimension(nsurf) c_dreq
real(kind(1d0)), dimension(maxnumberofgrids) tmin_id_grids
real(kind(1d0)) trans_site
real(kind(1d0)) s2
real(kind(1d0)), dimension(nvegsurf) theta_bioco2
integer gsmodel
real(kind(1d0)), dimension(3) ie_m
integer, dimension(nsurf) c_snowlimpat
real(kind(1d0)), dimension(nsurf) icefrac
integer, dimension(nvegsurf) c_alpha_bioco2
real(kind(1d0)), dimension(nsurf) snowwater
integer, dimension(nvegsurf) c_leafop2
real(kind(1d0)) albmin_grass
real(kind(1d0)) g5
real(kind(1d0)) snowfallcum
real(kind(1d0)) albmin_evetr
real(kind(1d0)), dimension(nsurf) state_id
real(kind(1d0)) vegfraction
real(kind(1d0)) pipecapacity
integer, dimension(nsurf) c_wetthresh
real(kind(1d0)), dimension(:), allocatable qn1_s_av_grids
real(kind(1d0)) qn1_av
real(kind(1d0)) kmax
integer, dimension(nsurf) c_soildepth
integer, dimension(nsurf) c_snowlimrem
integer, dimension(nsurf) c_obssmdepth
real(kind(1d0)), dimension(nvegsurf) laimax
real(kind(1d0)), dimension(:), allocatable qn1_av_grids
real(kind(1d0)), dimension(1) waterdepth
integer, dimension(nvegsurf) c_beta_enh_bioco2
real(kind(1d0)) snowlimpaved
real(kind(1d0)), dimension(2) popdensdaytime
real(kind(1d0)) z
integer, parameter bsoilsurf
real(kind(1d0)) th
real(kind(1d0)), dimension(nvegsurf) maxconductance
real(kind(1d0)) alt
real(kind(1d0)) snowalb
real(kind(1d0)) tempmeltfact
real(kind(1d0)), dimension(nsurf) soilstore_id
real(kind(1d0)) preciplimit
real(kind(1d0)) faibldg
real(kind(1d0)), dimension(nvegsurf) alpha_enh_bioco2
real(kind(1d0)) notused
integer, dimension(nvegsurf) c_sddfull
real(kind(1d0)), dimension(4, nvegsurf) laipower
subroutine suews_translateback(Gridiv, ir, irMax)
integer, dimension(nsurf) c_albmin
real(kind(1d0)) veg_fr
real(kind(1d0)) emis_snow
real(kind(1d0)), dimension(maxnumberofgrids) albevetr_id_grids
real(kind(1d0)), dimension(nvegsurf, maxnumberofgrids) lai_id_grids
real(kind(1d0)), dimension(maxnumberofgrids) lenday_id_grids
subroutine estm_translate(Gridiv)
real(kind(1d0)), dimension(nsurf) snowpacklimit
real(kind(1d0)) internalwateruse_h
real(kind(1d0)) surfacearea_ha
real(kind(1d0)), dimension(12, maxnumberofgrids) hdd_id_grids
real(kind(1d0)) soilrocks
real(kind(1d0)) albmax_evetr
real(kind(1d0)), dimension(nvegsurf) alpha_bioco2
real(kind(1d0)), dimension(maxnumberofgrids) albdectr_id_grids
integer snowuse
real(kind(1d0)) zdm
real(kind(1d0)) albevetr_id
integer, parameter nsurf
integer, dimension(nvegsurf) c_porositymin
real(kind(1d0)), dimension(nvegsurf) resp_a
real(kind(1d0)) g2
integer isec
real(kind(1d0)) wuareagrass_m2
real(kind(1d0)) tmax_id
real(kind(1d0)), dimension(:, :), allocatable surfacechar
real(kind(1d0)) h_maintain
real(kind(1d0)) irrfracbsoil
integer, dimension(nvegsurf) c_porositymax
integer, dimension(nvegsurf) c_laimax
integer, dimension(nsurf) c_drcoef2
integer, dimension(nsurf) c_obssmmax
real(kind(1d0)), dimension(nsurf) soilstorecap
real(kind(1d0)), dimension(nvegsurf) sdd_id
integer startdls
real(kind(1d0)), dimension(nsurf, maxnumberofgrids) icefrac_grids
real(kind(1d0)) zdm_in
real(kind(1d0)) nonwaterfraction
real(kind(1d0)) snowdensmin
integer aerodynamicresistancemethod
subroutine suews_translate(Gridiv, ir, iMB)
real(kind(1d0)) tair_av
real(kind(1d0)) faievetree
real(kind(1d0)) zzd
real(kind(1d0)) soildensity
real(kind(1d0)) g3
real(kind(1d0)) wuareadectr_m2
real(kind(1d0)) albmax_grass
integer, parameter conifsurf
real(kind(1d0)) irrfracbldgs
real(kind(1d0)), dimension(nsurf) sathydraulicconduct
real(kind(1d0)), dimension(nsurf) statelimit
integer, dimension(nvegsurf) c_leafgp1
integer id
real(kind(1d0)), dimension(nvegsurf) gddfull
real(kind(1d0)) rainmaxres
real(kind(1d0)), dimension(12) hdd_id
real(kind(1d0)) tau_f
real(kind(1d0)) lng
real(kind(1d0)) capmax_dec
real(kind(1d0)) irrfracdectr
real(kind(1d0)) irrfracpaved
real(kind(1d0)) dqnsdt
integer, dimension(nvegsurf) laitype
real(kind(1d0)) raincover
integer imin
real(kind(1d0)) snowalbmax
real(kind(1d0)), dimension(maxnumberofgrids) porosity_id_grids
real(kind(1d0)), dimension(nvegsurf) resp_b
real(kind(1d0)) pormax_dec
integer, dimension(nsurf) c_ksat
real(kind(1d0)) soildepthmeas
integer, dimension(nsurf) c_soildens
real(kind(1d0)) dectreeh
real(kind(1d0)) faut
real(kind(1d0)) tmin_id
real(kind(1d0)) lat
real(kind(1d0)), dimension(nsurf) emis
real(kind(1d0)), dimension(nsurf) soildepth
real(kind(1d0)) faitree
integer, parameter grasssurf
real(kind(1d0)) tau_a
real(kind(1d0)) radmeltfact
integer, dimension(nvegsurf) c_laimin
real(kind(1d0)) irrfracgrass
real(kind(1d0)) irrfracevetr
real(kind(1d0)) smcap
real(kind(1d0)), dimension(0:23, 2) snowprof_24hr
real(kind(1d0)), dimension(:, :), allocatable modeldailystate
integer iy
integer, dimension(nvegsurf) c_laieq
real(kind(1d0)) wuareaevetr_m2
integer, dimension(nsurf) cmod_state
real(kind(1d0)), dimension(nsurf) snowfrac
integer roughlenheatmethod
integer, dimension(nvegsurf) c_min_res_bioco2
integer, dimension(nsurf) c_stormin
integer, dimension(nsurf) c_soilstcap
real(kind(1d0)) g4
integer dt_since_start
real(kind(1d0)) faidectree
real(kind(1d0)), dimension(nvegsurf) laimin
integer, dimension(nsurf) c_obssnrfrac
real(kind(1d0)) surfacearea
integer, parameter nvegsurf
subroutine checkinitial
real(kind(1d0)), dimension(nsurf) wetthresh
integer, dimension(nsurf) c_albmax
real(kind(1d0)) pervfraction
integer, dimension(nsurf) cmod_snowfrac
character(len=10) gridid_text
real(kind(1d0)), dimension(3) ie_a
real(kind(1d0)) z0m_in
real(kind(1d0)) treeh
real(kind(1d0)), dimension(nvegsurf) lai_id
integer, dimension(nsurf) c_stormax
integer, dimension(nvegsurf) c_gsmax
real(kind(1d0)) irrfracwater
real(kind(1d0)), dimension(2) trafficrate
integer, dimension(nsurf) c_statelimit
real(kind(1d0)), dimension(9, maxnumberofgrids) wuday_id_grids
integer, dimension(nsurf) cmod_snowwaterstate
real(kind(1d0)) lenday_id
real(kind(1d0)) popdensnighttime
real(kind(1d0)) porosity_id
real(kind(1d0)), dimension(maxnumberofgrids) decidcap_id_grids
integer, dimension(nvegsurf) c_leafop1
integer, dimension(nvegsurf) c_resp_a
real(kind(1d0)) albmin_dectr
integer, dimension(nsurf) cmod_snowpack
real(kind(1d0)), dimension(nvegsurf) basete
real(kind(1d0)) albgrass_id
real(kind(1d0)), dimension(:), allocatable dqnsdt_grids
real(kind(1d0)) drainrt
real(kind(1d0)), dimension(nsurf) snowpack
integer, dimension(nvegsurf) c_basete
real(kind(1d0)), dimension(6, nsurf) storedrainprm
integer, dimension(nvegsurf) c_alpha_enh_bioco2
real(kind(1d0)), dimension(nvegsurf, maxnumberofgrids) sdd_id_grids
real(kind(1d0)), dimension(maxnumberofgrids) albgrass_id_grids
real(kind(1d0)), dimension(7) daywatper
real(kind(1d0)) dqndt
real(kind(1d0)), dimension(nsurf) sfr
real(kind(1d0)) snowlimbldg
real(kind(1d0)) g6
integer, dimension(nvegsurf) c_beta_bioco2
integer, dimension(nsurf) cmod_soilstate
real(kind(1d0)), dimension(9) wuday_id
real(kind(1d0)), dimension(nvegsurf) min_res_bioco2
real(kind(1d0)) crwmin
real(kind(1d0)) bldgh
integer it
real(kind(1d0)) snowdensmax
integer, parameter decidsurf
real(kind(1d0)) z0m
real(kind(1d0)), dimension(:), allocatable tair_av_grids
real(kind(1d0)) dectime
integer, parameter pavsurf
integer, dimension(nsurf) cmds_snowdens
real(kind(1d0)), dimension(nsurf) alb
real(kind(1d0)) preciplimitalb
real(kind(1d0)) pormin_dec
real(kind(1d0)) tl
real(kind(1d0)), dimension(nvegsurf) beta_enh_bioco2
integer, dimension(nvegsurf) c_baset
real(kind(1d0)), dimension(:), allocatable dqndt_grids
subroutine errorhint(errh, ProblemFile, VALUE, value2, valueI)
integer stabilitymethod
real(kind(1d0)), dimension(:, :, :), allocatable modeloutputdata
real(kind(1d0)) capmin_dec
integer, parameter bldgsurf
integer, parameter watersurf
integer, dimension(nsurf) c_drcoef1
real(kind(1d0)), dimension(nvegsurf) gdd_id
real(kind(1d0)), dimension(nvegsurf) sddfull
integer, dimension(nvegsurf) c_gddfull
real(kind(1d0)) crwmax
real(kind(1d0)), dimension(maxnumberofgrids) tmax_id_grids
real(kind(1d0)), dimension(2) qf0_beu
real(kind(1d0)), dimension(nvegsurf) baset
integer, dimension(nsurf) c_emis
real(kind(1d0)) snowalbmin
real(kind(1d0)), dimension(7) daywat
real(kind(1d0)) runofftowater
real(kind(1d0)) flowchange
real(kind(1d0)), dimension(nvegsurf, maxnumberofgrids) gdd_id_grids
real(kind(1d0)), dimension(nvegsurf) beta_bioco2
real(kind(1d0)) albdectr_id
integer, parameter ivdecid
real(kind(1d0)) decidcap_id
real(kind(1d0)) timezone
real(kind(1d0)) albmax_dectr
real(kind(1d0)) s1