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