SUEWS API Site
Documentation of SUEWS source code
suews_ctrl_input.f95
Go to the documentation of this file.
1 !Reading one line of meteorological forcing data in.
2 !Latest change:
3 ! Feb 2012, LJ: Input fluxes qh and qe changed _obs as well as qn1_obs ending
4 ! Oct 2014, LJ: Variables changed only be used in this part of code and these are passed to calling
5 ! function in MetArray.
6 ! Jan 2015, HCW: Precip_hr, wuh and LAI_hr changed for generic timesteps
7 ! Jan 2016, LJ: Removal of tabs
8 ! Feb 2017, HCW: Added file unit as argument so MetRead can be used for original met forcing file too
9 ! May 2018, TS: imporved the ability in downscaling datetime by introducing a datetime library
10 ! To Do:
11 ! - Check observed SM calculation
12 !---------------------------------------------------------------------------------------------------
13 SUBROUTINE metread(lfn, MetArray, InputmetFormat, ldown_option, NetRadiationMethod, &
14  snowUse, SMDMethod, SoilDepthMeas, SoilRocks, SoilDensity, SmCap)
15 
16  USE defaultnotused
17 
18  IMPLICIT NONE
19 
20  !INPUT
21  REAL(KIND(1d0)), DIMENSION(24)::MetArray !Array leaving the subroutine within
22  !each INTERVAL (defined in RunControl.nml)
23  ! - Met data now provided at a resolution of tstep, HCW Jan 2015
24 
25  REAL(KIND(1d0))::SmCap, &
26  SoilDepthMeas, & !Measured soil depth
27  SoilRocks, & !Rocks on ground
28  SoilDensity !Density of soil
29 
30  INTEGER::InputmetFormat, & !Format of the meteorological forcing file
31  ldown_option, & !Method of calculating Ldown
32  NetRadiationMethod, & !Method of calculating Q*
33  SMDMethod, & !Method of measured soil moisture
34  snowUse
35 
36  ! Variables read in
37  REAL(KIND(1d0))::avkdn, & !Average downwelling shortwave radiation
38  avrh, & !Average relative humidity
39  avu1, & !Average wind speed
40  dectime, & !Decimal time
41  fcld_obs, & !Cloud fraction observed
42  iy, & !Year
43  id, & !Day
44  it, & !Hour
45  imin, & !Minute
46  kdiff, & !Diffuse shortwave radiation
47  kdir, & !Direct shortwave radiation
48  LAI_obs, & !Overall LAI of the study area
49  ldown_obs, & !Downwelling longwave radiation
50  Precip, & !Rainfall [mm]
51  Pres_hPa, & !Station air pressure in hPa
52  Pres_kPa, & !Station air pressure in kPa
53  snowFrac_obs, & !Observed surface fraction of snow (between 0 and 1)
54  qe_obs, & !Observed latent heat flux
55  qf_obs, & !Observed antrhropogeni heat flux
56  qh_obs, & !Observed sensible heat flux
57  qn1_obs, & !Observed net all-wave radiation
58  qs_obs, & !Observed storage heat flux
59  Temp_C, & !Air temperature
60  wdir, & !Wind direction
61  wu_m3, & !Water use provided in met forcing file [m3]
62  xsmd !Measured soil moisture deficit
63 
64  INTEGER::iostat_var, lfn
65 
66  !-----------------------------------------------------------------------------------
67  !-----------------------------------------------------------------------------------
68 
69  IF (inputmetformat == 0) THEN !Default format using LUMPS only
70 
71  READ (lfn, *, iostat=iostat_var) iy, id, it, imin, qn1_obs, avu1, avrh, &
72  temp_c, wdir, pres_kpa, precip, avkdn, snowfrac_obs, ldown_obs, fcld_obs
73 
74  !Set other variables needed while running SUEWS to zero
75  qf_obs = nan
76  qs_obs = nan
77  qh_obs = nan
78  qe_obs = nan
79  xsmd = -99999
80  kdiff = nan
81  kdir = nan
82  wdir = nan
83 
84  ELSEIF (inputmetformat == 10) THEN !SUEWS reading
85  READ (lfn, *, iostat=iostat_var) iy, id, it, imin, qn1_obs, qh_obs, qe_obs, qs_obs, qf_obs, avu1, avrh, &
86  temp_c, pres_kpa, precip, avkdn, snowfrac_obs, ldown_obs, fcld_obs, &
87  wu_m3, xsmd, lai_obs, kdiff, kdir, wdir
88 
89  !write(*,*) 'In LUMPS_MetRead (1)'
90  !write(*,*) 'imin',imin
91  !write(*,*) 'it',it
92  !write(*,*) 'id',id
93  !write(*,*) 'iy',iy
94 
95  !Calculate observed soil moisture deficits from either volumetric or gravimetric SoilStates
96  IF (smdmethod == 1 .AND. xsmd /= -999) THEN !Soil moisture - volumetric
97  xsmd = (smcap - xsmd)*soildepthmeas*soilrocks
98  ELSEIF (smdmethod == 2 .AND. xsmd /= -999) THEN !Soil moisture -gravimetric
99  xsmd = (smcap - xsmd)*soildensity*soildepthmeas*soilrocks
100  ELSE
101  xsmd = -999
102  ENDIF
103 
104  ELSE
105  CALL errorhint(55, 'RunControl.nml, InputMetFormat not usable.', notused, notused, inputmetformat)
106  ENDIF
107 
108  !===============Meteorological variables reading done==========================
109  pres_hpa = pres_kpa*10. ! convert to hPa
110 
111  IF (iostat_var < 0) THEN
112  iostat_var = 0
113  CLOSE (lfn)
114  RETURN
115  ENDIF
116 
117  IF (avkdn < 0) THEN
118  CALL errorhint(27, 'Met Data: avKdn - needed for StoreDrainPrm. resistance, If present, check file not tab delimited', &
119  avkdn, dectime, notusedi)
120  !sg removed this is causing the problems with resistances
121  ! AvKdn=0 !Solar radiation cannot be lower than 1
122  ENDIF
123 
124  IF ((ldown_option == 1) .AND. (ldown_obs < 0)) THEN
125  CALL errorhint(27, 'Met Data: LWdn (ldown_obs) - impact Q* calc', ldown_obs, dectime, notusedi)
126 
127  ELSEIF (ldown_option == 2) THEN
128  IF (fcld_obs == -999.0 .OR. fcld_obs < 0 .OR. fcld_obs > 1) THEN
129  CALL errorhint(27, 'Met Data: flcd_obs - impacts LW & Q* radiation', fcld_obs, dectime, notusedi)
130  ENDIF
131  ENDIF
132 
133  IF (qn1_obs == -999 .AND. netradiationmethod == 0) THEN !If measured Q* is used and it is -999
134  CALL errorhint(27, 'Met Data: Q* - will impact everything', qn1_obs, dectime, notusedi)
135  ENDIF
136 
137  IF (avu1 <= 0) THEN !If wind speed is negative
138  CALL errorhint(27, 'Met Data: avU1 - impacts aeroydnamic resistances', avu1, dectime, notusedi)
139  ENDIF
140 
141  IF (temp_c < -50 .OR. temp_c > 60) THEN !If temperature unrealistic
142  CALL errorhint(27, 'Met Data: Temp_C - beyond what is expected', temp_c, dectime, notusedi)
143  ENDIF
144 
145  IF (avrh > 100 .OR. avrh < 1) THEN !If relative humidity larger than 100%
146  CALL errorhint(27, 'Met Data: avRH - beyond what is expected', avrh, dectime, notusedi)
147  ENDIF
148 
149  IF (pres_kpa < 80) THEN !If pressure too low
150  CALL errorhint(27, 'Met Data: Pres_kPa - too low - this could be fixed in model', pres_kpa, dectime, notusedi)
151  ENDIF
152 
153  IF (precip < 0) THEN !If rain in negative, set it to zero
154  CALL errorhint(27, 'Met Data: Precip - less than 0', precip, dectime, notusedi)
155  ENDIF
156 
157  IF (snowfrac_obs == nan) snowfrac_obs = 0
158 
159  IF (snowuse == 0 .AND. (snowfrac_obs < 0 .OR. snowfrac_obs > 1)) THEN
160  CALL errorhint(27, 'Met Data: snow not between [0 1]', snowfrac_obs, dectime, notusedi)
161  ENDIF
162 
163  IF (xsmd < 0 .AND. smdmethod == 1) THEN !If soil moisture deficit is zero
164  CALL errorhint(27, 'Met Data: xsmd - less than 0', xsmd, dectime, notusedi)
165  ENDIF
166 
167  !Create an array to be printed out.
168  metarray(1:24) = (/iy, id, it, imin, qn1_obs, qh_obs, qe_obs, qs_obs, qf_obs, avu1, &
169  avrh, temp_c, pres_hpa, precip, avkdn, snowfrac_obs, ldown_obs, &
170  fcld_obs, wu_m3, xsmd, lai_obs, kdiff, kdir, wdir/)
171 
172  !write(*,*) 'In LUMPS_MetRead (2)'
173  !write(*,*) 'imin',imin
174  !write(*,*) 'it',it
175  !write(*,*) 'id',id
176  !write(*,*) 'iy',iy
177 
178  RETURN
179 
180 END SUBROUTINE metread
181 
182 ! sg Feb 2012 -- moving all related suborutine together
183 ! subroutine run_control - checks that value is reasonable
184 ! subrotuine skipHeader - jumps over header rows in input files
185 ! Last modified:
186 ! LJ 27 Jan 2016 - Removal of tabs, cleaning of the code
187 !------------------------------------------------------------
188 
189 !Information for the run
190 MODULE run_info
191  IMPLICIT NONE
192  CHARACTER(len=90), DIMENSION(14)::text
193  INTEGER::lim0 = 0, lim1 = 1, lim2 = 2, lim4 = 4, lim3 = 3, lim6 = 6, lim8 = 8, lim12 = 12, lfn_us
194  LOGICAL ::file_qs
195 END MODULE run_info
196 
197 ! run_control
198 ! called from: LUMPS_initial
199 
200 SUBROUTINE run_control(eval, LowerLimit, Upperlimit)
201  ! ver - determines if value to be read is an integer or real and returns the value
202  ! if ver=-9 - then use integer
203  USE run_info
204  IMPLICIT NONE
205  INTEGER::eval, i, lowerlimit, upperlimit
206  CHARACTER(len=4)::check
207 
208  IF (file_qs) THEN
209 101 READ (lfn_us, *) check
210  WRITE (*, *) check
211  DO i = 1, 3
212  IF (check(i:i) == "#") THEN
213  ! write(*,*)check(i:i),i,"check"
214  GOTO 101
215  ELSE
216  backspace(lfn_us)
217  READ (lfn_us, *) eval
218  !write(*,*)eval, TEXT(1)
219  EXIT
220  ENDIF
221  ENDDO
222  ENDIF
223 
224  WRITE (12, 120) eval, text(1)
225 
226  IF (eval < lowerlimit .OR. eval > upperlimit) THEN
227  WRITE (*, *) "Value out of range"
228  WRITE (*, *) eval, text(1)
229  stop
230  ENDIF
231 
232  WRITE (*, 120) eval, text(1)
233 120 FORMAT(i4, 2x, a90)
234 
235  RETURN
236 END SUBROUTINE run_control
237 
238 ! Skip Headers------------------------------
239 SUBROUTINE skipheader(lfn, skip)
241  IMPLICIT NONE
242 
243  INTEGER::skip, lfn, i
244  DO i = 1, skip
245  READ (lfn, *, err=201, iostat=ios_out)
246  END DO
247 
248  RETURN
249 
250 201 reall = REAL(skip)
251  CALL errorhint(20, 'In SkipHeader subroutine.', reall, notused, ios_out)
252 END SUBROUTINE skipheader
253 
254 !-------------------------------------------------------------------------
255 SUBROUTINE inputheadercheck(FileName)
256  ! Checks columns in input files match the columns expected by model code
257  ! Model code columns are defined here
258  ! Latest update:
259  ! MH 21 Jun 2017 - Added parameters to SUEWS_AnthropogenicEmissions.txt
260  ! MH 16 Jun 2017 - Added SUEWS_BiogenCO2.txt
261  ! TS 02 Mar 2016 - AnOHM related variables added
262  ! LJ 27 Jan 2016 - Removal of tabs
263  ! LJ 07 July 2015 - snow albedo removed
264  ! HCW 12 Nov 2014
265  !-------------------------------------------------------------------------
266 
267  USE allocatearray
269  USE defaultnotused
270 
271  IMPLICIT NONE
272 
273  CHARACTER(len=50):: FileName
274 
275  ! ========== Define expected column names here ==========
276  ! =======================================================
277 
278  ! ========== SUEWS_NonVeg.txt =============
279  headernonveg_reqd(ci_code) = "Code"
280  headernonveg_reqd(ci_albmin) = "AlbedoMin"
281  headernonveg_reqd(ci_albmax) = "AlbedoMax"
282  headernonveg_reqd(ci_emis) = "Emissivity"
283  headernonveg_reqd(ci_stormin) = "StorageMin"
284  headernonveg_reqd(ci_stormax) = "StorageMax"
285  headernonveg_reqd(ci_wetthresh) = "WetThreshold"
286  headernonveg_reqd(ci_statelimit) = "StateLimit"
287  headernonveg_reqd(ci_dreq) = "DrainageEq"
288  headernonveg_reqd(ci_drcoef1) = "DrainageCoef1"
289  headernonveg_reqd(ci_drcoef2) = "DrainageCoef2"
290  headernonveg_reqd(ci_soiltcode) = "SoilTypeCode"
291  headernonveg_reqd(ci_snowlimpat) = "SnowLimPatch"
292  headernonveg_reqd(ci_snowlimrem) = "SnowLimRemove"
293  headernonveg_reqd(ci_ohmcode_swet) = "OHMCode_SummerWet"
294  headernonveg_reqd(ci_ohmcode_sdry) = "OHMCode_SummerDry"
295  headernonveg_reqd(ci_ohmcode_wwet) = "OHMCode_WinterWet"
296  headernonveg_reqd(ci_ohmcode_wdry) = "OHMCode_WinterDry"
297  headernonveg_reqd(ci_ohmthresh_sw) = "OHMThresh_SW"
298  headernonveg_reqd(ci_ohmthresh_wd) = "OHMThresh_WD"
299  headernonveg_reqd(ci_estmcode) = "ESTMCode"
300  headernonveg_reqd(ci_cpanohm) = "AnOHM_Cp" ! AnOHM TS
301  headernonveg_reqd(ci_kkanohm) = "AnOHM_Kk" ! AnOHM TS
302  headernonveg_reqd(ci_chanohm) = "AnOHM_Ch" ! AnOHM TS
303 
304  ! ========== SUEWS_Veg.txt ===============
305  headerveg_reqd(cp_code) = "Code"
306  headerveg_reqd(cp_albmin) = "AlbedoMin"
307  headerveg_reqd(cp_albmax) = "AlbedoMax"
308  headerveg_reqd(cp_emis) = "Emissivity"
309  headerveg_reqd(cp_stormin) = "StorageMin"
310  headerveg_reqd(cp_stormax) = "StorageMax"
311  headerveg_reqd(cp_wetthresh) = "WetThreshold"
312  headerveg_reqd(cp_statelimit) = "StateLimit"
313  headerveg_reqd(cp_dreq) = "DrainageEq"
314  headerveg_reqd(cp_drcoef1) = "DrainageCoef1"
315  headerveg_reqd(cp_drcoef2) = "DrainageCoef2"
316  headerveg_reqd(cp_soiltcode) = "SoilTypeCode"
317  headerveg_reqd(cp_snowlimpat) = "SnowLimPatch"
318  headerveg_reqd(cp_baset) = "BaseT"
319  headerveg_reqd(cp_basete) = "BaseTe"
320  headerveg_reqd(cp_gddfull) = "GDDFull"
321  headerveg_reqd(cp_sddfull) = "SDDFull"
322  headerveg_reqd(cp_laimin) = "LAIMin"
323  headerveg_reqd(cp_laimax) = "LAIMax"
324  headerveg_reqd(cp_porositymin) = "PorosityMin"
325  headerveg_reqd(cp_porositymax) = "PorosityMax"
326  headerveg_reqd(cp_gsmax) = "MaxConductance"
327  headerveg_reqd(cp_laieq) = "LAIEq"
328  headerveg_reqd(cp_leafgp1) = "LeafGrowthPower1"
329  headerveg_reqd(cp_leafgp2) = "LeafGrowthPower2"
330  headerveg_reqd(cp_leafop1) = "LeafOffPower1"
331  headerveg_reqd(cp_leafop2) = "LeafOffPower2"
332  headerveg_reqd(cp_ohmcode_swet) = "OHMCode_SummerWet"
333  headerveg_reqd(cp_ohmcode_sdry) = "OHMCode_SummerDry"
334  headerveg_reqd(cp_ohmcode_wwet) = "OHMCode_WinterWet"
335  headerveg_reqd(cp_ohmcode_wdry) = "OHMCode_WinterDry"
336  headerveg_reqd(cp_ohmthresh_sw) = "OHMThresh_SW"
337  headerveg_reqd(cp_ohmthresh_wd) = "OHMThresh_WD"
338  headerveg_reqd(cp_estmcode) = "ESTMCode"
339  headerveg_reqd(cp_cpanohm) = "AnOHM_Cp" ! AnOHM TS
340  headerveg_reqd(cp_kkanohm) = "AnOHM_Kk" ! AnOHM TS
341  headerveg_reqd(cp_chanohm) = "AnOHM_Ch" ! AnOHM TS
342  headerveg_reqd(cp_biogenco2code) = "BiogenCO2Code"
343 
344  ! ========== SUEWS_Water.txt ==================
345  headerwater_reqd(cw_code) = "Code"
346  headerwater_reqd(cw_albmin) = "AlbedoMin"
347  headerwater_reqd(cw_albmax) = "AlbedoMax"
348  headerwater_reqd(cw_emis) = "Emissivity"
349  headerwater_reqd(cw_stormin) = "StorageMin"
350  headerwater_reqd(cw_stormax) = "StorageMax"
351  headerwater_reqd(cw_wetthresh) = "WetThreshold"
352  headerwater_reqd(cw_statelimit) = "StateLimit"
353  headerwater_reqd(cw_waterdepth) = "WaterDepth"
354  headerwater_reqd(cw_dreq) = "DrainageEq"
355  headerwater_reqd(cw_drcoef1) = "DrainageCoef1"
356  headerwater_reqd(cw_drcoef2) = "DrainageCoef2"
357  headerwater_reqd(cw_ohmcode_swet) = "OHMCode_SummerWet"
358  headerwater_reqd(cw_ohmcode_sdry) = "OHMCode_SummerDry"
359  headerwater_reqd(cw_ohmcode_wwet) = "OHMCode_WinterWet"
360  headerwater_reqd(cw_ohmcode_wdry) = "OHMCode_WinterDry"
361  headerwater_reqd(cw_ohmthresh_sw) = "OHMThresh_SW"
362  headerwater_reqd(cw_ohmthresh_wd) = "OHMThresh_WD"
363  headerwater_reqd(cw_estmcode) = "ESTMCode"
364  headerwater_reqd(cw_cpanohm) = "AnOHM_Cp" ! AnOHM TS
365  headerwater_reqd(cw_kkanohm) = "AnOHM_Kk" ! AnOHM TS
366  headerwater_reqd(cw_chanohm) = "AnOHM_Ch" ! AnOHM TS
367 
368  ! ========== SUEWS_Snow.txt ===================
369  headersnow_reqd(cs_code) = "Code"
370  headersnow_reqd(cs_snowrmfactor) = "RadMeltFactor"
371  headersnow_reqd(cs_snowtmfactor) = "TempMeltFactor"
372  headersnow_reqd(cs_snowalbmin) = "AlbedoMin"
373  headersnow_reqd(cs_snowalbmax) = "AlbedoMax"
374  headersnow_reqd(cs_snowemis) = "Emissivity"
375  headersnow_reqd(cs_snowtau_a) = "tau_a"
376  headersnow_reqd(cs_snowtau_f) = "tau_f"
377  headersnow_reqd(cs_snowplimalb) = "PrecipLimAlb"
378  headersnow_reqd(cs_snowsdmin) = "SnowDensMin"
379  headersnow_reqd(cs_snowsdmax) = "SnowDensMax"
380  headersnow_reqd(cs_snowtau_r) = "tau_r"
381  headersnow_reqd(cs_snowcrwmin) = "CRWMin"
382  headersnow_reqd(cs_snowcrwmax) = "CRWMax"
383  headersnow_reqd(cs_snowplimsnow) = "PrecipLimSnow"
384  headersnow_reqd(cs_ohmcode_swet) = "OHMCode_SummerWet"
385  headersnow_reqd(cs_ohmcode_sdry) = "OHMCode_SummerDry"
386  headersnow_reqd(cs_ohmcode_wwet) = "OHMCode_WinterWet"
387  headersnow_reqd(cs_ohmcode_wdry) = "OHMCode_WinterDry"
388  headersnow_reqd(cs_ohmthresh_sw) = "OHMThresh_SW"
389  headersnow_reqd(cs_ohmthresh_wd) = "OHMThresh_WD"
390  headersnow_reqd(cs_estmcode) = "ESTMCode"
391  headersnow_reqd(cs_cpanohm) = "AnOHM_Cp" ! AnOHM TS
392  headersnow_reqd(cs_kkanohm) = "AnOHM_Kk" ! AnOHM TS
393  headersnow_reqd(cs_chanohm) = "AnOHM_Ch" ! AnOHM TS
394 
395  ! ========== SUEWS_Soil.txt ===================
396  headersoil_reqd(cso_code) = "Code"
397  headersoil_reqd(cso_soildepth) = "SoilDepth"
398  headersoil_reqd(cso_soilstcap) = "SoilStoreCap"
399  headersoil_reqd(cso_ksat) = "SatHydraulicCond"
400  headersoil_reqd(cso_soildens) = "SoilDensity"
401  headersoil_reqd(cso_soilinfrate) = "InfiltrationRate"
402  headersoil_reqd(cso_obssmdepth) = "OBS_SMDepth"
403  headersoil_reqd(cso_obssmmax) = "OBS_SMCap"
404  headersoil_reqd(cso_obssnrfrac) = "OBS_SoilNotRocks"
405 
406  ! ========== SUEWS_Conductance.txt ============
407  headercond_reqd(cc_code) = "Code"
408  headercond_reqd(cc_gsg1) = "G1"
409  headercond_reqd(cc_gsg2) = "G2"
410  headercond_reqd(cc_gsg3) = "G3"
411  headercond_reqd(cc_gsg4) = "G4"
412  headercond_reqd(cc_gsg5) = "G5"
413  headercond_reqd(cc_gsg6) = "G6"
414  headercond_reqd(cc_gsth) = "TH"
415  headercond_reqd(cc_gstl) = "TL"
416  headercond_reqd(cc_gss1) = "S1"
417  headercond_reqd(cc_gss2) = "S2"
418  headercond_reqd(cc_gskmax) = "Kmax"
419  headercond_reqd(cc_gsmodel) = "gsModel"
420 
421  ! ========== SUEWS_OHMCoefficients.txt ========
426 
427  ! ========== SUEWS_ESTMCoefficients.txt ========
475  headerestmcoefficients_reqd(ce_alb_ibld) = "Internal_albedo"
476  headerestmcoefficients_reqd(ce_em_ibld) = "Internal_emissivity"
477  headerestmcoefficients_reqd(ce_ch_iwall) = "Internal_CHwall"
478  headerestmcoefficients_reqd(ce_ch_iroof) = "Internal_CHroof"
479  headerestmcoefficients_reqd(ce_ch_ibld) = "Internal_CHbld"
480 
481  ! ========== SUEWS_AnthropogenicEmission.txt ======
492  headeranthropogenic_reqd(ca_ahslopeheating_wd) = "AHSlope_Heating_WD"
493  headeranthropogenic_reqd(ca_ahslopeheating_we) = "AHSlope_Heating_WE"
494  headeranthropogenic_reqd(ca_ahslopecooling_wd) = "AHSlope_Cooling_WD"
495  headeranthropogenic_reqd(ca_ahslopecooling_we) = "AHSlope_Cooling_WE"
496  headeranthropogenic_reqd(ca_tcriticheating_wd) = "TCritic_Heating_WD"
497  headeranthropogenic_reqd(ca_tcriticheating_we) = "TCritic_Heating_WE"
498  headeranthropogenic_reqd(ca_tcriticcooling_wd) = "TCritic_Cooling_WD"
499  headeranthropogenic_reqd(ca_tcriticcooling_we) = "TCritic_Cooling_WE"
500  headeranthropogenic_reqd(ca_enprofwd) = "EnergyUseProfWD"
501  headeranthropogenic_reqd(ca_enprofwe) = "EnergyUseProfWE"
502  headeranthropogenic_reqd(ca_co2mwd) = "ActivityProfWD"
503  headeranthropogenic_reqd(ca_co2mwe) = "ActivityProfWE"
513  headeranthropogenic_reqd(ca_frfossilfuel_heat) = "FrFossilFuel_Heat"
514  headeranthropogenic_reqd(ca_frfossilfuel_nonheat) = "FrFossilFuel_NonHeat"
515  headeranthropogenic_reqd(ca_ef_umolco2perj) = "EF_umolCO2perJ"
517  headeranthropogenic_reqd(ca_fcef_v_kgkmwd) = "FcEF_v_kgkmWD"
518  headeranthropogenic_reqd(ca_fcef_v_kgkmwe) = "FcEF_v_kgkmWE"
519  headeranthropogenic_reqd(ca_co2pointsource) = "CO2PointSource"
520  headeranthropogenic_reqd(ca_trafficunits) = "TrafficUnits"
521 
522  ! ========== SUEWS_Irrigation.txt =============
524  headerirrigation_reqd(cir_iestart) = "Ie_start"
525  headerirrigation_reqd(cir_ieend) = "Ie_end"
526  headerirrigation_reqd(cir_intwu) = "InternalWaterUse"
528  headerirrigation_reqd(cir_h_ponding) = "H_maintain"
535  headerirrigation_reqd(cir_daywat1) = "DayWat(1)"
536  headerirrigation_reqd(cir_daywat2) = "DayWat(2)"
537  headerirrigation_reqd(cir_daywat3) = "DayWat(3)"
538  headerirrigation_reqd(cir_daywat4) = "DayWat(4)"
539  headerirrigation_reqd(cir_daywat5) = "DayWat(5)"
540  headerirrigation_reqd(cir_daywat6) = "DayWat(6)"
541  headerirrigation_reqd(cir_daywat7) = "DayWat(7)"
542  headerirrigation_reqd(cir_daywatper1) = "DayWatPer(1)"
543  headerirrigation_reqd(cir_daywatper2) = "DayWatPer(2)"
544  headerirrigation_reqd(cir_daywatper3) = "DayWatPer(3)"
545  headerirrigation_reqd(cir_daywatper4) = "DayWatPer(4)"
546  headerirrigation_reqd(cir_daywatper5) = "DayWatPer(5)"
547  headerirrigation_reqd(cir_daywatper6) = "DayWatPer(6)"
548  headerirrigation_reqd(cir_daywatper7) = "DayWatPer(7)"
549 
550  ! ========== SUEWS_Profiles.txt ===============
551  headerprofiles_reqd(cpr_code) = "Code"
561  headerprofiles_reqd(cpr_hours(10)) = "9"
562  headerprofiles_reqd(cpr_hours(11)) = "10"
563  headerprofiles_reqd(cpr_hours(12)) = "11"
564  headerprofiles_reqd(cpr_hours(13)) = "12"
565  headerprofiles_reqd(cpr_hours(14)) = "13"
566  headerprofiles_reqd(cpr_hours(15)) = "14"
567  headerprofiles_reqd(cpr_hours(16)) = "15"
568  headerprofiles_reqd(cpr_hours(17)) = "16"
569  headerprofiles_reqd(cpr_hours(18)) = "17"
570  headerprofiles_reqd(cpr_hours(19)) = "18"
571  headerprofiles_reqd(cpr_hours(20)) = "19"
572  headerprofiles_reqd(cpr_hours(21)) = "20"
573  headerprofiles_reqd(cpr_hours(22)) = "21"
574  headerprofiles_reqd(cpr_hours(23)) = "22"
575  headerprofiles_reqd(cpr_hours(24)) = "23"
576 
577  ! ========== SUEWS_WithinGridWaterDist.txt ====
587  headerwgwaterdist_reqd(cwg_tosoilstore) = "ToSoilStore"
588 
589  ! ========== SUEWS_BiogenCO2.txt ======
590  headerbiogen_reqd(cb_code) = "Code"
591  headerbiogen_reqd(cb_alpha) = "alpha"
592  headerbiogen_reqd(cb_beta) = "beta"
593  headerbiogen_reqd(cb_theta) = "theta"
594  headerbiogen_reqd(cb_alpha_enh) = "alpha_enh"
595  headerbiogen_reqd(cb_beta_enh) = "beta_enh"
596  headerbiogen_reqd(cb_resp_a) = "resp_a"
597  headerbiogen_reqd(cb_resp_b) = "resp_b"
598  headerbiogen_reqd(cb_min_r) = "min_respi"
599 
600  ! =======================================================
601 
602  !write(*,*) 'Checking header for ', FileName
603  ! Check columns in input files match model code
604 
605  IF (filename == 'SUEWS_NonVeg.txt') THEN
606  IF (any(headernonveg_file /= headernonveg_reqd)) THEN
607  WRITE (*, *) headernonveg_file == headernonveg_reqd
608  WRITE (*, *) headernonveg_file
609  WRITE (*, *) headernonveg_reqd
610  CALL errorhint(56, 'Names or order of columns in SUEWS_NonVeg.txt does not match model code.', notused, notused, notusedi)
611  ENDIF
612 
613  ELSEIF (filename == 'SUEWS_Veg.txt') THEN
614  IF (any(headerveg_file /= headerveg_reqd)) THEN
615  WRITE (*, *) headerveg_file == headerveg_reqd
616  CALL errorhint(56, 'Names or order of columns in SUEWS_Veg.txt does not match model code.', notused, notused, notusedi)
617  ENDIF
618 
619  ELSEIF (filename == 'SUEWS_Water.txt') THEN
620  IF (any(headerwater_file /= headerwater_reqd)) THEN
621  WRITE (*, *) headerwater_file == headerwater_reqd
622  WRITE (*, *) headerwater_file
623  WRITE (*, *) headerwater_reqd
624  CALL errorhint(56, 'Names or order of columns in SUEWS_Water.txt does not match model code.', notused, notused, notusedi)
625  ENDIF
626 
627  ELSEIF (filename == 'SUEWS_Snow.txt') THEN
628  IF (any(headersnow_file /= headersnow_reqd)) THEN
629  WRITE (*, *) headersnow_file == headersnow_reqd
630  CALL errorhint(56, 'Names or order of columns in SUEWS_Snow.txt does not match model code.', notused, notused, notusedi)
631  ENDIF
632 
633  ELSEIF (filename == 'SUEWS_Soil.txt') THEN
634  IF (any(headersoil_file /= headersoil_reqd)) THEN
635  WRITE (*, *) headersoil_file == headersoil_reqd
636  CALL errorhint(56, 'Names or order of columns in SUEWS_Soil.txt does not match model code.', notused, notused, notusedi)
637  ENDIF
638 
639  ELSEIF (filename == 'SUEWS_Conductance.txt') THEN
640  IF (any(headercond_file /= headercond_reqd)) THEN
641  WRITE (*, *) headercond_file == headercond_reqd
642  CALL errorhint(56, 'Names or order of columns in SUEWS_Cond.txt does not match model code.', notused, notused, notusedi)
643  ENDIF
644 
645  ELSEIF (filename == 'SUEWS_OHMCoefficients.txt') THEN
648  CALL errorhint(56, 'Names or order of columns in SUEWS_OHMCoefficients.txt does not match model code.', &
650  ENDIF
651 
652  ELSEIF (filename == 'SUEWS_ESTMCoefficients.txt') THEN
655  CALL errorhint(56, 'Names or order of columns in SUEWS_ESTMCoefficients.txt does not match model code.', &
657  ENDIF
658 
659  ELSEIF (filename == 'SUEWS_AnthropogenicEmission.txt') THEN
662  CALL errorhint(56, 'Names or order of columns in SUEWS_AnthropogenicEmission.txt does not match model code.', &
664  ENDIF
665 
666  ELSEIF (filename == 'SUEWS_Irrigation.txt') THEN
669  WRITE (*, *) headerirrigation_file
670  WRITE (*, *) headerirrigation_reqd
671  CALL errorhint(56, 'Names or order of columns in SUEWS_Irrigation.txt does not match model code.', &
673  ENDIF
674 
675  ELSEIF (filename == 'SUEWS_Profiles.txt') THEN
676  IF (any(headerprofiles_file /= headerprofiles_reqd)) THEN
678  CALL errorhint(56, 'Names or order of columns in SUEWS_Profiles.txt does not match model code.', &
680  ENDIF
681 
682  ELSEIF (filename == 'SUEWS_WithinGridWaterDist.txt') THEN
685  CALL errorhint(56, 'Names or order of columns in SUEWS_WithinGridWaterDist.txt does not match model code.', &
687  ENDIF
688 
689  ELSEIF (filename == 'SUEWS_BiogenCO2.txt') THEN
690  IF (any(headerbiogen_file /= headerbiogen_reqd)) THEN
691  WRITE (*, *) headerbiogen_file == headerbiogen_reqd
692  CALL errorhint(56, 'Names or order of columns in SUEWS_BiogenCO2.txt does not match model code.', &
694  ENDIF
695 
696  ELSE
697  WRITE (*, *) 'Problem in subroutine InputHeaderCheck. File header not specified in model code for ', filename
699 
700  ENDIF
701 
702 ENDSUBROUTINE inputheadercheck
703 
704 !-------------------------------------------------------------------------
705 !
706 ! TS 05 Jul 2018: No longer needed as interpolation is done through specific subroutines at each required instant
707 !Interpolates hourly profiles provided in SUEWS_Profiles.txt
708 ! to resolution of the model timestep
709 ! HCW 06 Feb 2015
710 !===================================================================================
711 ! SUBROUTINE SUEWS_InterpHourlyProfiles(Gridiv,TstepP_ID,SurfChar_HrProf)
712 !
713 ! USE allocateArray
714 ! USE ColNamesInputFiles
715 ! USE sues_data
716 !
717 ! IMPLICIT NONE
718 !
719 ! INTEGER:: i,j, ii !Used to count over hours and sub-hourly timesteps
720 ! INTEGER:: Gridiv, TstepP_ID
721 ! INTEGER,DIMENSION(24):: SurfChar_HrProf
722 ! REAL(KIND(1d0)):: deltaProf !Change in hourly profiles per model timestep
723 !
724 ! ! Copy value for first hour
725 ! TstepProfiles(Gridiv,TstepP_ID,1) = SurfaceChar(Gridiv,SurfChar_HrProf(1))
726 ! DO i=1,24
727 ! j = (i+1)
728 ! IF(i == 24) j = 1 !If last hour of day, loop round to first hour of day for interpolation
729 ! deltaProf = ((SurfaceChar(Gridiv,SurfChar_HrProf(j)) - SurfaceChar(Gridiv,SurfChar_HrProf(i))))/nsh_real
730 ! DO ii=1,nsh
731 ! IF((nsh*(i-1)+ii+1) < (23*nsh+nsh+1)) THEN
732 ! TstepProfiles(Gridiv,TstepP_ID,(nsh*(i-1)+ii+1)) = SurfaceChar(Gridiv,SurfChar_HrProf(i)) + deltaProf*ii
733 ! ENDIF
734 ! ENDDO
735 ! ENDDO
736 !
737 ! endsubroutine SUEWS_InterpHourlyProfiles
738 !===================================================================================
739 
740 !===================================================================================
741 ! get interpolated profile values at specified time
742 ! NO normalisation performed
743 FUNCTION get_prof_spectime_inst(Hour, Min, Sec, Prof_24h) RESULT(Prof_CurrTime)
745  IMPLICIT NONE
746 
747  INTEGER :: i, j !Used to count over hours and sub-hourly timesteps
748  INTEGER, INTENT(IN) :: Hour, Min, Sec
749  INTEGER :: total_sec, SecPerHour
750  REAL(KIND(1d0)), DIMENSION(0:23), INTENT(IN) :: Prof_24h
751  REAL(KIND(1d0)):: deltaProf !Change in hourly profiles per model timestep
752  REAL(KIND(1d0)) :: Prof_CurrTime
753 
754  total_sec = min*60 + sec
755  secperhour = 3600
756 
757  i = hour
758  j = i + 1
759  IF (j == 24) j = 0
760 
761  deltaprof = (prof_24h(j) - prof_24h(i))/secperhour
762  prof_currtime = prof_24h(hour) + deltaprof*total_sec
763 
764 END FUNCTION get_prof_spectime_inst
765 
766 !===================================================================================
767 ! get interpolated profile values at specified time
768 ! normalise so the AVERAGE of the multipliers is equal to 1
769 FUNCTION get_prof_spectime_mean(Hour, Min, Sec, Prof_24h) RESULT(Prof_CurrTime)
771  IMPLICIT NONE
772 
773  INTEGER :: i, j !Used to count over hours and sub-hourly timesteps
774  INTEGER, INTENT(IN) :: Hour, Min, Sec
775  INTEGER :: total_sec, SecPerHour
776  REAL(KIND(1d0)), DIMENSION(0:23), INTENT(IN) :: Prof_24h
777  REAL(KIND(1d0)), DIMENSION(0:23):: Prof_24h_mean
778  REAL(KIND(1d0)):: deltaProf !Change in hourly profiles per model timestep
779  REAL(KIND(1d0)) :: Prof_CurrTime
780 
781  total_sec = min*60 + sec
782  secperhour = 3600
783 
784  prof_24h_mean = merge(prof_24h/(sum(prof_24h)/size(prof_24h, dim=1)), 0.d0, sum(prof_24h) /= 0) ! prevent zero-division
785  ! print*, Prof_24h_mean
786 
787  i = hour
788  j = i + 1
789  IF (j == 24) j = 0
790 
791  deltaprof = (prof_24h_mean(j) - prof_24h_mean(i))/secperhour
792 
793  ! print*, deltaProf,total_sec
794  prof_currtime = prof_24h_mean(i) + deltaprof*total_sec
795 
796 END FUNCTION get_prof_spectime_mean
797 
798 !===================================================================================
799 ! get interpolated profile values at specified time
800 ! normalise so the SUM of the multipliers is equal to 1
801 FUNCTION get_prof_spectime_sum(Hour, Min, Sec, Prof_24h, dt) RESULT(Prof_CurrTime)
803  IMPLICIT NONE
804 
805  INTEGER :: i, j !Used to count over hours and sub-hourly timesteps
806  INTEGER, INTENT(IN) :: Hour, Min, Sec, dt
807  INTEGER :: total_sec, SecPerHour
808  REAL(KIND(1d0)), DIMENSION(0:23), INTENT(IN) :: Prof_24h
809  REAL(KIND(1d0)), DIMENSION(0:23):: Prof_24h_sum
810  REAL(KIND(1d0)):: deltaProf !Change in hourly profiles per model timestep
811  REAL(KIND(1d0)) :: Prof_CurrTime
812 
813  total_sec = min*60 + sec
814  secperhour = 3600
815 
816  prof_24h_sum = merge(prof_24h/(sum(prof_24h)), 0.d0, sum(prof_24h) /= 0) ! prevent zero-division
817 
818  i = hour
819  j = i + 1
820  IF (j == 24) j = 0
821 
822  deltaprof = (prof_24h_sum(j) - prof_24h_sum(i))/secperhour
823  prof_currtime = prof_24h_sum(hour) + deltaprof*total_sec
824  prof_currtime = prof_currtime*dt/secperhour
825 
826 END FUNCTION get_prof_spectime_sum
827 !===================================================================================
828 
829 ! Subroutines for matching codes in the input files
830 ! could re-write as a generic function later...
831 
832 SUBROUTINE codematchohm(Gridiv, is, SWWD)
833  ! Matches OHM coefficients a1, a2, a3 via OHM codes
834  ! for summer/winter wet/dry conditions
835  ! HCW 03 Nov 2014
836  ! ---------------------------------------------------------
837 
838  USE allocatearray
839  USE initial
841  USE defaultnotused
842 
843  IMPLICIT NONE
844 
845  INTEGER:: gridiv
846  INTEGER:: is
847  CHARACTER(len=4):: SWWD
848 
849  iv5 = 0 ! Reset iv5 to zero
850 
851  IF (swwd == 'SWet') THEN
852 
853  DO iv5 = 1, nlinesohmcoefficients
854  IF (ohmcoefficients_coeff(iv5, co_code) == surfacechar(gridiv, c_ohmcode_swet(is))) THEN
855  EXIT
856  ELSEIF (iv5 == nlinesohmcoefficients) THEN
857  WRITE (*, *) 'Program stopped! OHM code (summer wet)', surfacechar(gridiv, c_ohmcode_swet(is)), &
858  'not found in OHM_Coefficients.txt for surface', is, '.'
859  CALL errorhint(57, 'Cannot find OHM code (summer wet)', surfacechar(gridiv, c_ohmcode_swet(is)), notused, notusedi)
860  ENDIF
861  ENDDO
862 
863  ELSEIF (swwd == 'SDry') THEN
864 
865  DO iv5 = 1, nlinesohmcoefficients
866  IF (ohmcoefficients_coeff(iv5, co_code) == surfacechar(gridiv, c_ohmcode_sdry(is))) THEN
867  EXIT
868  ELSEIF (iv5 == nlinesohmcoefficients) THEN
869  WRITE (*, *) 'Program stopped! OHM code (summer dry)', surfacechar(gridiv, c_ohmcode_sdry(is)), &
870  'not found in OHM_Coefficients.txt for surface', is, '.'
871  CALL errorhint(57, 'Cannot find OHM code (summer dry)', surfacechar(gridiv, c_ohmcode_sdry(is)), notused, notusedi)
872  ENDIF
873  ENDDO
874 
875  ELSEIF (swwd == 'WWet') THEN
876 
877  DO iv5 = 1, nlinesohmcoefficients
878  IF (ohmcoefficients_coeff(iv5, co_code) == surfacechar(gridiv, c_ohmcode_wwet(is))) THEN
879  EXIT
880  ELSEIF (iv5 == nlinesohmcoefficients) THEN
881  WRITE (*, *) 'Program stopped! OHM code (winter wet)', surfacechar(gridiv, c_ohmcode_wwet(is)), &
882  'not found in OHM_Coefficients.txt for surface', is, '.'
883  CALL errorhint(57, 'Cannot find OHM code (winter wet)', surfacechar(gridiv, c_ohmcode_wwet(is)), notused, notusedi)
884  ENDIF
885  ENDDO
886 
887  ELSEIF (swwd == 'WDry') THEN
888 
889  DO iv5 = 1, nlinesohmcoefficients
890  IF (ohmcoefficients_coeff(iv5, co_code) == surfacechar(gridiv, c_ohmcode_wdry(is))) THEN
891  EXIT
892  ELSEIF (iv5 == nlinesohmcoefficients) THEN
893  WRITE (*, *) 'Program stopped! OHM code (winter dry)', surfacechar(gridiv, c_ohmcode_wdry(is)), &
894  'not found in OHM_Coefficients.txt for surface', is, '.'
895  CALL errorhint(57, 'Cannot find OHM code (winter dry)', surfacechar(gridiv, c_ohmcode_wdry(is)), notused, notusedi)
896  ENDIF
897  ENDDO
898 
899  ELSE
900  WRITE (*, *) 'Problem with CodeMatchOHM (in SUEWS_CodeMatch.f95). ', swwd, ' not recognised. Needs to be one of: ', &
901  'SWet = Summer Wet, SDry = Summer Dry, WWet = WinterWet, WDry = Winter Dry. N.B. Case sensitive. '
902  stop
903  ENDIF
904 
905  RETURN
906 ENDSUBROUTINE codematchohm
907 ! ---------------------------------------------------------
908 
909 SUBROUTINE codematchestm(Gridiv, is)
910  ! Matches ESTM coefficients via ESTM code
911  ! Modified HCW 16 Jun 2016 - for SUEWS surface types
912  ! - removed summer/winter wet/dry option
913  ! S.O. 04 Feb 2016
914  ! ---------------------------------------------------------
915 
916  USE allocatearray
917  USE initial
919  USE defaultnotused
920 
921  IMPLICIT NONE
922 
923  INTEGER:: gridiv
924  INTEGER:: is
925 
926  iv5 = 0 ! Reset iv5 to zero
927 
929  IF (estmcoefficients_coeff(iv5, ce_code) == surfacechar(gridiv, c_estmcode(is))) THEN
930  EXIT
931  ELSEIF (iv5 == nlinesestmcoefficients) THEN
932  WRITE (*, *) 'Program stopped! ESTM code', surfacechar(gridiv, c_estmcode(is)), &
933  'not found in ESTM_Coefficients.txt for surface', is, '.'
934  CALL errorhint(57, 'Cannot find ESTM code', surfacechar(gridiv, c_estmcode(is)), notused, notusedi)
935  ENDIF
936  ENDDO
937 
938  RETURN
939 ENDSUBROUTINE codematchestm
940 ! ---------------------------------------------------------
941 
942 SUBROUTINE codematchestm_class(Gridiv, is, ii)
943  ! Matches ESTM coefficients via ESTM codes in SiteSelect for Paved and Bldgs ESTM classes
944  ! HCW 16 Jun 2016
945  ! ---------------------------------------------------------
946 
947  USE allocatearray
948  USE initial
950  USE defaultnotused
951 
952  IMPLICIT NONE
953 
954  INTEGER:: gridiv
955  INTEGER:: is, ii
956 
957  iv5 = 0 ! Reset iv5 to zero
958 
959  IF (is == bldgsurf) THEN
962  EXIT
963  ELSEIF (iv5 == nlinesestmcoefficients) THEN
964  WRITE (*, *) 'Program stopped! ESTM code', surfacechar(gridiv, c_code_estmclass_bldgs(ii)), &
965  'not found in ESTM_Coefficients.txt for surface', is, '.'
966  CALL errorhint(57, 'Cannot find ESTM code', surfacechar(gridiv, c_code_estmclass_bldgs(ii)), notused, notusedi)
967  ENDIF
968  ENDDO
969  ELSEIF (is == pavsurf) THEN
972  EXIT
973  ELSEIF (iv5 == nlinesestmcoefficients) THEN
974  WRITE (*, *) 'Program stopped! ESTM code', surfacechar(gridiv, c_code_estmclass_paved(ii)), &
975  'not found in ESTM_Coefficients.txt for surface', is, '.'
976  CALL errorhint(57, 'Cannot find ESTM code', surfacechar(gridiv, c_code_estmclass_paved(ii)), notused, notusedi)
977  ENDIF
978  ENDDO
979  ELSE
980  WRITE (*, *) 'Problem with CodeMatchESTM_Class (in SUEWS_ctrl_input.f95). ', is, ' not correct. Needs to be either ', &
981  '1 = Paved surfaced, 2 = Bldgs surfaces.'
982  stop
983  ENDIF
984  RETURN
985 ENDSUBROUTINE codematchestm_class
986 ! ---------------------------------------------------------
987 
988 SUBROUTINE codematchprof(Gridiv, SurfaceCharCodeCol)
989  ! Matches Soil characteristics via codes in *SurfaceChar*
990  ! for energy use/water use/snow clearing
991  ! HCW 20 Nov 2014
992  ! ---------------------------------------------------------
993 
994  USE allocatearray
995  USE initial
997  USE defaultnotused
998 
999  IMPLICIT NONE
1000 
1001  INTEGER:: Gridiv
1002  INTEGER:: SurfaceCharCodeCol
1003 
1004  iv5 = 0 ! Reset iv5 to zero
1005 
1006  DO iv5 = 1, nlinesprofiles
1007  IF (profiles_coeff(iv5, cpr_code) == surfacechar(gridiv, surfacecharcodecol)) THEN
1008  EXIT
1009  ELSEIF (iv5 == nlinesprofiles) THEN
1010  WRITE (*, *) 'Program stopped! Profile code ', surfacechar(gridiv, surfacecharcodecol), 'not found in SUEWS_Profiles.txt.'
1011  CALL errorhint(57, 'Cannot find code in SUEWS_Profiles.txt', surfacechar(gridiv, surfacecharcodecol), notused, notusedi)
1012  ENDIF
1013  ENDDO
1014 
1015  RETURN
1016 ENDSUBROUTINE codematchprof
1017 ! ---------------------------------------------------------
1018 
1019 SUBROUTINE codematchdist(rr, CodeCol, codeColSameSurf)
1020  ! Matches within-grid water distribution via codes
1021  ! Checks water cannot flow from one surface to the same surface
1022  ! Checks water distribution fractions sum to 1.
1023  ! HCW 10 Nov 2014
1024  ! ---------------------------------------------------------
1025 
1026  USE allocatearray
1027  USE initial
1028  USE colnamesinputfiles
1029  USE defaultnotused
1030 
1031  IMPLICIT NONE
1032 
1033  INTEGER:: rr
1034  INTEGER:: codeCol, codeColSameSurf
1035 
1036  iv5 = 0 ! Reset iv5 to zero
1037 
1038  DO iv5 = 1, nlineswgwaterdist
1039  IF (wgwaterdist_coeff(iv5, cwg_code) == siteselect(rr, codecol)) THEN
1040  EXIT
1041  ELSEIF (iv5 == nlineswgwaterdist) THEN
1042  WRITE (*, *) 'Program stopped! Within-grid water distribution code ', siteselect(rr, codecol), &
1043  'not found in SUEWS_WaterDistWithinGrid.txt.'
1044  CALL errorhint(57, 'Cannot find code in SUEWS_WaterDistWithinGrid.txt', siteselect(rr, codecol), notused, notusedi)
1045  ENDIF
1046  ENDDO
1047 
1048  ! Check water flow to same surface is zero (previously in RunControlByGridByYear in SUEWS_Initial.f95)
1049  IF (wgwaterdist_coeff(iv5, codecolsamesurf) /= 0) THEN
1050  CALL errorhint(8, 'Diagonal elements should be zero as water cannot move from one surface to the same surface.', &
1051  wgwaterdist_coeff(iv5, codecolsamesurf), notused, notusedi)
1052  ENDIF
1053 
1054  !! QUESTION: MODIFY THIS?
1055  ! Check water either moves to runoff or soilstore, but not to both
1056  ! Model returns an error if both ToRunoff and ToSoilStore are non-zero.
1057  !! - Probably should remove this...
1058  !! - Also look at SUEWS_translate, as the non-zero value goes into WaterDist
1060  CALL errorhint(9, 'One of these (ToRunoff,ToSoilStore) should be zero.', &
1062  ENDIF
1063 
1064  !! Also do for water surface once implemented
1065  IF (codecol /= c_wgwatercode) THEN ! Except for Water surface
1066  ! Check total water distribution from each surface adds up to 1
1067  IF (sum(wgwaterdist_coeff(iv5, cwg_topaved:cwg_tosoilstore)) > 1.0000001 &
1068  .OR. sum(wgwaterdist_coeff(iv5, cwg_topaved:cwg_tosoilstore)) < 0.9999999) THEN
1069  CALL errorhint(8, 'Total water distribution from each surface should add up to 1.', &
1071  ENDIF
1072  ENDIF
1073 
1074  RETURN
1075 ENDSUBROUTINE codematchdist
1076 ! ---------------------------------------------------------
1077 
1078 SUBROUTINE codematchnonveg(rr, CodeCol)
1079  ! Matches Impervious characteristics via codes in SiteSelect
1080  ! HCW 20 Nov 2014
1081  ! ---------------------------------------------------------
1082 
1083  USE allocatearray
1084  USE initial
1085  USE colnamesinputfiles
1086  USE defaultnotused
1087 
1088  IMPLICIT NONE
1089 
1090  INTEGER:: rr
1091  INTEGER:: codeCol
1092 
1093  iv5 = 0 ! Reset iv5 to zero
1094 
1095  DO iv5 = 1, nlinesnonveg
1096  IF (nonveg_coeff(iv5, ci_code) == siteselect(rr, codecol)) THEN
1097  EXIT
1098  ELSEIF (iv5 == nlinesnonveg) THEN
1099  WRITE (*, *) 'Program stopped! NonVeg code ', siteselect(rr, codecol), 'not found in SUEWS_NonVeg.txt.'
1100  CALL errorhint(57, 'Cannot find code in SUEWS_NonVeg.txt', siteselect(rr, codecol), notused, notusedi)
1101  ENDIF
1102  ENDDO
1103 
1104  RETURN
1105 ENDSUBROUTINE codematchnonveg
1106 ! ---------------------------------------------------------
1107 
1108 SUBROUTINE codematchveg(rr, CodeCol)
1109  ! Matches Pervious characteristics via codes in SiteSelect
1110  ! HCW 20 Nov 2014
1111  ! ---------------------------------------------------------
1112 
1113  USE allocatearray
1114  USE initial
1115  USE colnamesinputfiles
1116  USE defaultnotused
1117 
1118  IMPLICIT NONE
1119 
1120  INTEGER:: rr
1121  INTEGER:: codeCol
1122 
1123  iv5 = 0 ! Reset iv5 to zero
1124 
1125  DO iv5 = 1, nlinesveg
1126  IF (veg_coeff(iv5, cp_code) == siteselect(rr, codecol)) THEN
1127  EXIT
1128  ELSEIF (iv5 == nlinesveg) THEN
1129  WRITE (*, *) 'Program stopped! Veg code ', siteselect(rr, codecol), 'not found in SUEWS_Vegs.txt.'
1130  CALL errorhint(57, 'Cannot find code in SUEWS_Veg.txt', siteselect(rr, codecol), notused, notusedi)
1131  ENDIF
1132  ENDDO
1133 
1134  RETURN
1135 ENDSUBROUTINE codematchveg
1136 ! ---------------------------------------------------------
1137 
1138 SUBROUTINE codematchwater(rr, CodeCol)
1139  ! Matches Water characteristics via codes in SiteSelect
1140  ! HCW 20 Nov 2014
1141  ! ---------------------------------------------------------
1142 
1143  USE allocatearray
1144  USE initial
1145  USE colnamesinputfiles
1146  USE defaultnotused
1147 
1148  IMPLICIT NONE
1149 
1150  INTEGER:: rr
1151  INTEGER:: codeCol
1152 
1153  iv5 = 0 ! Reset iv5 to zero
1154 
1155  DO iv5 = 1, nlineswater
1156  IF (water_coeff(iv5, cw_code) == siteselect(rr, codecol)) THEN
1157  EXIT
1158  ELSEIF (iv5 == nlineswater) THEN
1159  WRITE (*, *) 'Program stopped! Water code ', siteselect(rr, codecol), 'not found in SUEWS_Water.txt.'
1160  CALL errorhint(57, 'Cannot find code in SUEWS_Water.txt', siteselect(rr, codecol), notused, notusedi)
1161  ENDIF
1162  ENDDO
1163 
1164  RETURN
1165 ENDSUBROUTINE codematchwater
1166 ! ---------------------------------------------------------
1167 
1168 SUBROUTINE codematchsnow(rr, CodeCol)
1169  ! Matches Snow characteristics via codes in SiteSelect
1170  ! HCW 20 Nov 2014
1171  ! ---------------------------------------------------------
1172 
1173  USE allocatearray
1174  USE initial
1175  USE colnamesinputfiles
1176  USE defaultnotused
1177 
1178  IMPLICIT NONE
1179 
1180  INTEGER:: rr
1181  INTEGER:: codeCol
1182 
1183  iv5 = 0 ! Reset iv5 to zero
1184 
1185  DO iv5 = 1, nlinessnow
1186  IF (snow_coeff(iv5, cs_code) == siteselect(rr, codecol)) THEN
1187  EXIT
1188  ELSEIF (iv5 == nlinessnow) THEN
1189  WRITE (*, *) 'Program stopped! Snow code ', siteselect(rr, codecol), 'not found in SUEWS_Snow.txt.'
1190  CALL errorhint(57, 'Cannot find code in SUEWS_Snow.txt', siteselect(rr, codecol), notused, notusedi)
1191  ENDIF
1192  ENDDO
1193 
1194  RETURN
1195 ENDSUBROUTINE codematchsnow
1196 ! ---------------------------------------------------------
1197 
1198 SUBROUTINE codematchconductance(rr, CodeCol)
1199  ! Matches Conductance characteristics via codes in SiteSelect
1200  ! HCW 20 Nov 2014
1201  ! ---------------------------------------------------------
1202 
1203  USE allocatearray
1204  USE initial
1205  USE colnamesinputfiles
1206  USE defaultnotused
1207 
1208  IMPLICIT NONE
1209 
1210  INTEGER:: rr
1211  INTEGER:: codeCol
1212 
1213  iv5 = 0 ! Reset iv5 to zero
1214 
1215  DO iv5 = 1, nlinesconductance
1216  IF (conductance_coeff(iv5, cc_code) == siteselect(rr, codecol)) THEN
1217  EXIT
1218  ELSEIF (iv5 == nlinesconductance) THEN
1219  WRITE (*, *) 'Program stopped! Conductance code ', siteselect(rr, codecol), 'not found in SUEWS_Conductance.txt.'
1220  CALL errorhint(57, 'Cannot find code in SUEWS_Conductance.txt', siteselect(rr, codecol), notused, notusedi)
1221  ENDIF
1222  ENDDO
1223 
1224  RETURN
1225 ENDSUBROUTINE codematchconductance
1226 ! ---------------------------------------------------------
1227 
1228 SUBROUTINE codematchanthropogenic(rr, CodeCol)
1229  ! Matches AnthropogenicHeat characteristics via codes in SiteSelect
1230  ! HCW 20 Nov 2014
1231  ! MH 21 Jun 2017
1232  ! ---------------------------------------------------------
1233 
1234  USE allocatearray
1235  USE initial
1236  USE colnamesinputfiles
1237  USE defaultnotused
1238 
1239  IMPLICIT NONE
1240 
1241  INTEGER:: rr
1242  INTEGER:: codeCol
1243 
1244  iv5 = 0 ! Reset iv5 to zero
1245 
1246  DO iv5 = 1, nlinesanthropogenic
1247  IF (anthropogenic_coeff(iv5, ca_code) == siteselect(rr, codecol)) THEN
1248  EXIT
1249  ELSEIF (iv5 == nlinesanthropogenic) THEN
1250  WRITE (*, *) 'Program stopped! Anthropogenic code ', siteselect(rr, codecol), &
1251  'not found in SUEWS_AnthropogenicEmission.txt.'
1252  CALL errorhint(57, 'Cannot find code in SUEWS_AnthropogenicEmission.txt', &
1253  siteselect(rr, codecol), notused, notusedi)
1254  ENDIF
1255  ENDDO
1256 
1257  RETURN
1258 ENDSUBROUTINE codematchanthropogenic
1259 ! ---------------------------------------------------------
1260 
1261 SUBROUTINE codematchirrigation(rr, CodeCol)
1262  ! Matches Irrigation characteristics via codes in SiteSelect
1263  ! HCW 20 Nov 2014
1264  ! ---------------------------------------------------------
1265 
1266  USE allocatearray
1267  USE initial
1268  USE colnamesinputfiles
1269  USE defaultnotused
1270 
1271  IMPLICIT NONE
1272 
1273  INTEGER:: rr
1274  INTEGER:: codeCol
1275 
1276  iv5 = 0 ! Reset iv5 to zero
1277 
1278  DO iv5 = 1, nlinesirrigation
1279  IF (irrigation_coeff(iv5, cir_code) == siteselect(rr, codecol)) THEN
1280  EXIT
1281  ELSEIF (iv5 == nlinesirrigation) THEN
1282  WRITE (*, *) 'Program stopped! Irrigation code ', siteselect(rr, codecol), 'not found in SUEWS_Irrigation.txt.'
1283  CALL errorhint(57, 'Cannot find code in SUEWS_Irrigation.txt', siteselect(rr, codecol), notused, notusedi)
1284  ENDIF
1285  ENDDO
1286 
1287  RETURN
1288 ENDSUBROUTINE codematchirrigation
1289 ! ---------------------------------------------------------
1290 
1291 SUBROUTINE codematchsoil(Gridiv, SurfaceCharCodeCol)
1292  ! Matches Soil characteristics via codes in *SurfaceChar*
1293  ! HCW 20 Nov 2014
1294  ! ---------------------------------------------------------
1295 
1296  USE allocatearray
1297  USE initial
1298  USE colnamesinputfiles
1299  USE defaultnotused
1300 
1301  IMPLICIT NONE
1302 
1303  INTEGER:: Gridiv
1304  INTEGER:: SurfaceCharCodeCol
1305 
1306  iv5 = 0 ! Reset iv5 to zero
1307 
1308  DO iv5 = 1, nlinessoil
1309  IF (soil_coeff(iv5, cso_code) == surfacechar(gridiv, surfacecharcodecol)) THEN
1310  EXIT
1311  ELSEIF (iv5 == nlinessoil) THEN
1312  WRITE (*, *) 'Program stopped! Soil code ', surfacechar(gridiv, surfacecharcodecol), 'not found in SUEWS_Soil.txt.'
1313  CALL errorhint(57, 'Cannot find code in SUEWS_Soil.txt', surfacechar(gridiv, surfacecharcodecol), notused, notusedi)
1314  ENDIF
1315  ENDDO
1316 
1317  RETURN
1318 ENDSUBROUTINE codematchsoil
1319 ! ---------------------------------------------------------
1320 
1321 SUBROUTINE codematchbiogen(Gridiv, SurfaceCharCodeCol)
1322  ! Matches Biogen characteristics via codes in *SuraceChar*
1323  ! MH 16 Jun 2017
1324  ! ---------------------------------------------------------
1325 
1326  USE allocatearray
1327  USE initial
1328  USE colnamesinputfiles
1329  USE defaultnotused
1330 
1331  IMPLICIT NONE
1332 
1333  INTEGER:: Gridiv
1334  INTEGER:: SurfaceCharCodeCol
1335 
1336  iv5 = 0 ! Reset iv5 to zero
1337 
1338  DO iv5 = 1, nlinesbiogen
1339  IF (biogen_coeff(iv5, cb_code) == surfacechar(gridiv, surfacecharcodecol)) THEN
1340  EXIT
1341  ELSEIF (iv5 == nlinesbiogen) THEN
1342  WRITE (*, *) 'Program stopped! Biogen code ', surfacechar(gridiv, surfacecharcodecol), 'not found in SUEWS_BiogenCO2.txt.'
1343  CALL errorhint(57, 'Cannot find code in SUEWS_BiogenCO2.txt', surfacechar(gridiv, surfacecharcodecol), notused, notusedi)
1344  ENDIF
1345  ENDDO
1346 
1347  RETURN
1348 ENDSUBROUTINE codematchbiogen
1349 ! ---------------------------------------------------------
1350 
1352  !========================================================================================
1353  ! Disaggregation of meteorological forcing data
1354  ! Code to disaggregate met forcing data from resolution provided to the model time-step
1355  ! Created: HCW 10 Feb 2017
1356  !
1357  ! ---- Key for MetDisaggMethod settings ----
1358  !10 -> linear disaggregation for timestamps representing the end of the averaging period
1359  !20 -> linear disaggregation for instantaneous variables (e.g. air temperature, humidity, pressure, wind speed in WFDEI dataset)
1360  !100 -> evenly distribute rainfall among all subintervals in a rainy interval
1361  !101 -> evenly distribute rainfall among RainAmongN subintervals in a rainy interval
1362  ! - requires RainAmongN to be set in RunControl.nml
1363  ! If KdownZen = 1 -> include additional zenith check in kdown disaggregation
1364  !
1365  ! N.B. wdir downscaling is currently not implemented
1366  !
1367  !========================================================================================
1368 
1369  USE allocatearray
1370  USE colnamesinputfiles
1371  USE data_in
1372  USE initial
1374 
1375  IMPLICIT NONE
1376 
1377 CONTAINS
1378 
1379  !======================================================================================
1380  SUBROUTINE disaggregatemet(iBlock, igrid)
1381  ! Subroutine to disaggregate met forcing data to model time-step
1382  ! HCW 10 Feb 2017
1383  !======================================================================================
1384 
1385  USE sues_data
1386  USE defaultnotused
1387 
1388  IMPLICIT NONE
1389 
1390  INTEGER:: lunit = 100
1391  INTEGER:: tdiff !Time difference (in minutes) between first and second rows of original met forcing file
1392  INTEGER:: i, ii !counter
1393  INTEGER:: iBlock, igrid
1394  INTEGER, DIMENSION(Nper):: seq1Nper
1395  INTEGER, DIMENSION(nsd):: seq1nsd
1396  INTEGER, DIMENSION(nColumnsMetForcingData):: MetDisaggMethod ! Stores method to use for disaggregating met data
1397  REAL(KIND(1d0)), DIMENSION(nColumnsMetForcingData):: MetArrayOrig
1398  REAL(KIND(1d0)), DIMENSION(ReadLinesOrigMetData*Nper, ncolumnsMetForcingData):: Met_tt
1399  REAL(KIND(1d0)), DIMENSION(ReadLinesOrigMetData*Nper):: Met_tt_kdownAdj
1400  CHARACTER(LEN=9), DIMENSION(ncolumnsMetForcingData):: HeaderMet
1401  CHARACTER(LEN=10*ncolumnsMetForcingData):: HeaderMetOut
1402  ! REAL(KIND(1d0)),DIMENSION(ReadLinesOrigMetData):: dectimeOrig
1403  ! REAL(KIND(1d0)),DIMENSION(ReadLinesOrigMetData*Nper):: dectimeDscd, dectimeFast
1404  REAL(KIND(1d0)), DIMENSION(ReadLinesOrigMetData*Nper):: dectimeFast
1405  REAL(KIND(1d0)), DIMENSION(ReadLinesOrigMetData*Nper):: idectime ! sun position at middle of time-step before
1406 
1407  ! INTEGER, DIMENSION(Nper):: temp_iy, temp_id, temp_ih, temp_im, temp_ihm ! temorary varaibles for disaggragation
1408  ! REAL(KIND(1d0)), DIMENSION(Nper):: temp_dectime ! temorary varaibles for disaggragation
1409  ! INTEGER :: Days_of_Year
1410  ! INTEGER, DIMENSION(Nper)::ndays_iy ! number of days in iy
1411 
1412  ! Allocate and initialise arrays to receive original forcing data --------------------
1416  metfordisagg(:, :) = -999
1417  metfordisaggprev(:) = -999
1418  metfordisaggnext(:) = -999
1419  ! Initialise array to receive disaggregated data
1420  met_tt = -999
1421  ! Intialise disaggregation method
1422  metdisaggmethod = -999
1423 
1424  ! Generate useful sequences
1425  seq1nper = (/(i, i=1, nper, 1)/)
1426  seq1nsd = (/(i, i=1, nsd, 1)/)
1427 
1428  ! Get methods to use for disaggregation from RunControl
1429  IF (diagnosedisagg == 1) WRITE (*, *) 'DisaggMethod: ', disaggmethod, 'RainDisaggMethod:', raindisaggmethod
1430  IF (disaggmethod == 1) THEN
1431  metdisaggmethod(:) = 10 !linear disaggregation of averages
1432  ELSEIF (disaggmethod == 2) THEN
1433  metdisaggmethod(:) = 20 !linear disaggregation of instantaneous values
1434  ELSEIF (disaggmethod == 3) THEN !WFDEI set up, where T, Q, pres, U are instantaneous
1435  metdisaggmethod(:) = 10 !linear disaggregation of averages
1436  metdisaggmethod(10:13) = 20 !linear disagg instantaneous values for U, RH, Tair, pres
1437  ELSE
1438  CALL errorhint(2, 'Problem in SUEWS_MetDisagg: DisaggMethod value should be 1, 2, or 3', &
1440  ENDIF
1441  ! Set rainfall
1442  metdisaggmethod(14) = raindisaggmethod
1443 
1444  ! Read data ---------------------------------------------------------------------
1445  IF (diagnosedisagg == 1) WRITE (*, *) 'Reading file: ', trim(fileorigmet)
1446  OPEN (lunit, file=trim(fileorigmet), status='old')
1447  ! CALL skipHeader(lunit,SkipHeaderMet) !Skip header -> read header instead
1448  READ (lunit, *) headermet
1449  !write(*,*) HeaderMet
1450  ! Skip over lines that have already been read and downscaled
1451  IF (skippedlinesorig > 0) THEN
1452  DO i = 1, skippedlinesorig - 1 ! minus 1 here because last line of last block needs to be read again
1453  READ (lunit, *)
1454  ENDDO
1455  ! Read in last line of previous block
1456  CALL metread(lunit, metarrayorig, inputmetformat, ldown_option, netradiationmethod, &
1458  metfordisaggprev(1:ncolumnsmetforcingdata) = metarrayorig
1459  ENDIF
1460  ! print*, 'MetForDisagg',MetForDisagg(1:3,1:4)
1461  ! print*, 'ReadLinesOrigMetDataMax',ReadLinesOrigMetDataMax
1462  ! Read in current block
1463  DO i = 1, readlinesorigmetdatamax
1464  CALL metread(lunit, metarrayorig, inputmetformat, ldown_option, netradiationmethod, &
1466  metfordisagg(i, 1:ncolumnsmetforcingdata) = metarrayorig
1467  ENDDO
1468  ! print*, 'MetForDisagg',MetForDisagg(1:3,1:4)
1469  ! Read in first line of next block (except for last block)
1470  IF (iblock /= readblocksorigmetdata) THEN
1471  CALL metread(lunit, metarrayorig, inputmetformat, ldown_option, netradiationmethod, &
1473  metfordisaggnext(1:ncolumnsmetforcingdata) = metarrayorig
1474  ENDIF
1475  CLOSE (lunit)
1476 
1477  ! Check resolution of original met forcing data -------------------------------------
1478  ! Find time difference (in minutes) between first and second row
1479  tdiff = int(metfordisagg(2, 4) - metfordisagg(1, 4)) !Try using minutes
1480  IF (tdiff == 0) tdiff = int(metfordisagg(2, 3) - metfordisagg(1, 3))*60 !If no difference in minutes, try using hours
1481  IF (tdiff < 0) THEN !If time difference is negative (e.g. change of day), instead use second and third row
1482  tdiff = int(metfordisagg(3, 4) - metfordisagg(2, 4))
1483  IF (tdiff == 0) tdiff = int(metfordisagg(3, 3) - metfordisagg(2, 3))*60 !If no difference in minutes, try using hours
1484  ENDIF
1485  ! Check actual resolution matches specified input resolution
1486  IF (tdiff /= resolutionfilesin/60) THEN
1487  CALL errorhint(2, 'Problem in SUEWS_MetDisagg: timestamps in met forcing file inconsistent with ResolutionFilesIn', &
1488  REAL(ResolutionFilesIn, KIND(1d0)), NotUsed, tdiff*60)
1489  ENDIF
1490 
1491  ! Check file only contains a single year --------------------------------------------
1492  ! Very last data point is allowed to be (should be) timestamped with following year
1493  IF (any(metfordisagg(1:(readlinesorigmetdatamax - 1), 1) /= metfordisagg(1, 1))) THEN
1494  CALL errorhint(3, 'Problem in SUEWS_MetDisagg: multiple years found in original met forcing file.', &
1495  metfordisagg(1, 1), notused, notusedi)
1496  ENDIF
1497 
1498  ! Disaggregate time columns ---------------------------------------------------------
1499  IF (diagnose == 1) WRITE (*, *) 'Disaggregating met forcing data (', trim(fileorigmet), ') to model time-step...'
1500 
1501  CALL disaggregatedatetime(metfordisagg(:, 1:4), tstep, nper, readlinesorigmetdatamax, met_tt(:, 1:4))
1502 
1503  ! Disaggregate other columns --------------------------------------------------------
1504  DO ii = 5, ncolumnsmetforcingdata
1505  IF (ii == 14) THEN !Do something different for rainfall and snowfall (if present)
1506  IF (metdisaggmethod(14) == 100) THEN
1508  IF (all(metfordisagg(:, 16) == -999)) THEN
1509  met_tt(:, 16) = -999
1510  ELSE
1512  ENDIF
1513  ELSEIF (metdisaggmethod(14) == 101) THEN
1514  IF (rainamongn == -999) THEN
1515  CALL errorhint(2, 'Problem in SUEWS_MetDisagg: RainDisaggMethod requires RainAmongN', &
1516  REAL(RainAmongN, KIND(1d0)), NotUsed, RainDisaggMethod)
1517  ELSEIF (rainamongn > nper) THEN
1518  CALL errorhint(2, 'Problem in SUEWS_MetDisagg: RainAmongN > Nper', REAL(Nper, KIND(1d0)), NotUsed, RainAmongN)
1519  ELSE
1520  met_tt(:, 14) = disaggp_amongn(metfordisagg(:, 14), &
1522  IF (all(metfordisagg(:, 16) == -999)) THEN
1523  met_tt(:, 16) = -999
1524  ELSE
1525  met_tt(:, 16) = disaggp_amongn(metfordisagg(:, 16), &
1527  ENDIF
1528  ENDIF
1529  ELSEIF (metdisaggmethod(14) == 102) THEN
1530  IF (all(multrainamongn == -999)) THEN
1531  CALL errorhint(2, 'Problem in SUEWS_MetDisagg: RainDisaggMethod requires MultRainAmongN', &
1532  REAL(MultRainAmongN(1), KIND(1d0)), NotUsed, RainDisaggMethod)
1533  ELSEIF (all(multrainamongnupperi == -999)) THEN
1534  CALL errorhint(2, 'Problem in SUEWS_MetDisagg: RainDisaggMethod requires MultRainAmongNUpperI', &
1535  multrainamongnupperi(1), notused, raindisaggmethod)
1536  ELSEIF (any(multrainamongn > nper)) THEN
1537  CALL errorhint(2, 'Problem in SUEWS_MetDisagg: MultRainAmongN > Nper', REAL(Nper, KIND(1d0)), NotUsed, &
1538  MAXVAL(multrainamongn))
1539  ELSE
1542  IF (all(metfordisagg(:, 16) == -999)) THEN
1543  met_tt(:, 16) = -999
1544  ELSE
1547  ENDIF
1548  ENDIF
1549  ELSE
1550  WRITE (*, *) 'Disaggregation code for rain not recognised'
1551  ENDIF
1552  ELSEIF (ii == 24) THEN !wind direction disaggregation not coded yet...
1553  IF (any(metfordisagg(:, ii) /= -999)) THEN
1554  WRITE (*, *) 'Disaggregation of wind direction not currently implemented!'
1555  ENDIF
1556  ELSE
1557  IF (all(metfordisagg(:, ii) == -999)) THEN
1558  !IF(DiagnoseDisagg==1) write(*,*) 'No data for col.', ii
1559  met_tt(:, ii) = -999
1560  ELSE
1561  met_tt(:, ii) = disagg_lin(metfordisagg(:, ii), metfordisaggprev(ii), metfordisaggnext(ii), metdisaggmethod(ii), &
1563  ENDIF
1564  ENDIF
1565  ENDDO
1566 
1567  ! Adjust kdown disaggregation using zenith angle
1568  IF (kdownzen == 1) THEN
1569  IF (diagnosedisagg == 1) WRITE (*, *) 'Adjusting disaggregated kdown using zenith angle'
1570  met_tt_kdownadj(:) = met_tt(:, 15)
1571  ! Translate location data from SurfaceChar to find solar angles
1572  lat = surfacechar(igrid, c_lat)
1573  lng = surfacechar(igrid, c_lng)
1574  timezone = surfacechar(igrid, c_tz)
1575  alt = surfacechar(igrid, c_alt)
1576  ! Calculate dectime at downscaled time-step
1577  dectimefast(:) = met_tt(:, 2) + met_tt(:, 3)/24.0 + met_tt(:, 4)/(60.0*24.0)
1578  idectime = dectimefast - halftimestep! sun position at middle of timestep before
1579  DO i = 1, (readlinesorigmetdatamax*nper)
1580  CALL narp_cal_sunposition(met_tt(i, 2), idectime(i), timezone, lat, lng, alt, azimuth, zenith_deg)
1581  ! If sun below horizon, set disaggregated kdown to zero
1582  IF (zenith_deg > 90) THEN
1583  !write(*,*) Met_tt(i,1:4)
1584  met_tt_kdownadj(i) = 0.0
1585  ENDIF
1586  ENDDO
1587  ! Redistribute kdown over each day
1588  DO i = 1, (readlinesorigmetdatamax*nper/nsd) ! Loop over each day
1589  met_tt_kdownadj((i - 1)*nsd + seq1nsd) = &
1590  met_tt_kdownadj((i - 1)*nsd + seq1nsd) &
1591  *sum(met_tt((i - 1)*nsd + seq1nsd, 15)) &
1592  /sum(met_tt_kdownadj((i - 1)*nsd + seq1nsd))
1593  ENDDO
1594  ! Copy adjusted kdown back to Met_tt array
1595  met_tt(:, 15) = met_tt_kdownadj(:)
1596  ENDIF
1597 
1598  ! Copy disaggregated data to MetForcingDataArray
1599  metforcingdata(:, 1:24, gridcounter) = met_tt(:, 1:24)
1600 
1601  ! If snow is -999, set to zero (also in LUMPS_metRead.f95)
1602  IF (all(metforcingdata(:, 16, gridcounter) == -999)) metforcingdata(:, 16, gridcounter) = 0
1603 
1604  ! Undo pressure conversion again for writing out
1605  met_tt(:, 13) = met_tt(:, 13)/10.0
1606 
1607  ! Write out disaggregated file ------------------------------------------------------
1608  IF (keeptstepfilesin == 1) THEN
1609  IF (iblock == 1) THEN
1610  ! Prepare header
1611  DO i = 1, ncolumnsmetforcingdata
1612  IF (i == 1) THEN
1613  headermetout = adjustl(headermet(i))
1614  ELSE
1615  headermetout = trim(headermetout)//' '//adjustl(headermet(i))
1616  ENDIF
1617  ENDDO
1618  ! Write out header
1619  OPEN (78, file=trim(filedscdmet), err=112)
1620  WRITE (78, '(a)') headermetout
1621  ELSE
1622  OPEN (78, file=trim(filedscdmet), position='append')!,err=112)
1623  ENDIF
1624  ! Write out data
1625  DO i = 1, (readlinesorigmetdatamax*nper)
1626  WRITE (78, 303) (int(met_tt(i, ii)), ii=1, 4), met_tt(i, 5:ncolumnsmetforcingdata)
1627  ENDDO
1628  !IF(iBlock == ReadBlocksOrigMetData) THEN
1629  ! WRITE(78,'(i2)') -9
1630  ! WRITE(78,'(i2)') -9
1631  !ENDIF
1632  CLOSE (78) !Close output file
1633  ENDIF
1634 
1635 303 FORMAT((i4, 1x), 3(i3, 1x), 9(f12.6, 1x), (f9.4, 1x), 10(f9.4, 1x)) !Allows 4 dp for rainfall
1636 
1637  ! Deallocate arrays -----------------------------------------------------------------
1638  DEALLOCATE (metfordisagg)
1639  DEALLOCATE (metfordisaggprev)
1640  DEALLOCATE (metfordisaggnext)
1641 
1642  RETURN
1643 
1644 112 CALL errorhint(52, trim(filedscdmet), notused, notused, notusedi)
1645 
1646  ENDSUBROUTINE disaggregatemet
1647  !======================================================================================
1648 
1649  !======================================================================================
1650  SUBROUTINE disaggregateestm(iBlock)
1651  ! Subroutine to disaggregate met forcing data to model time-step
1652  ! HCW 10 Feb 2017
1653  !======================================================================================
1654 
1655  USE sues_data
1656  USE defaultnotused
1657 
1658  IMPLICIT NONE
1659 
1660  INTEGER:: lunit = 101
1661  INTEGER:: tdiff !Time difference (in minutes) between first and second rows of original met forcing file
1662  INTEGER:: i, ii !counter
1663  INTEGER:: iBlock
1664  INTEGER, DIMENSION(NperESTM):: seq1NperESTM
1665  INTEGER, DIMENSION(nsd):: seq1nsd
1666  INTEGER, DIMENSION(ncolsESTMdata):: ESTMDisaggMethod ! Stores method to use for disaggregating met data
1667  REAL(KIND(1d0)), DIMENSION(ncolsESTMdata):: ESTMArrayOrig
1668  REAL(KIND(1d0)), DIMENSION(ReadLinesOrigESTMData*NperESTM, ncolsESTMdata):: ESTM_tt
1669  CHARACTER(LEN=9), DIMENSION(ncolsESTMdata):: HeaderESTM
1670  CHARACTER(LEN=10*ncolsESTMdata):: HeaderESTMOut
1671  ! REAL(KIND(1d0)),DIMENSION(ReadLinesOrigESTMData):: dectimeOrig
1672  ! REAL(KIND(1d0)),DIMENSION(ReadLinesOrigESTMData*NperESTM):: dectimeDscd
1673  INTEGER::iostat_var
1674 
1675  ! INTEGER, DIMENSION(NperESTM):: temp_iy, temp_id, temp_ih, temp_im, temp_ihm
1676 
1677  ! Allocate and initialise arrays to receive original forcing data --------------------
1679  ALLOCATE (estmfordisaggprev(ncolsestmdata))
1680  ALLOCATE (estmfordisaggnext(ncolsestmdata))
1681  estmfordisagg(:, :) = -999
1682  estmfordisaggprev(:) = -999
1683  estmfordisaggnext(:) = -999
1684  ! Initialise array to receive disaggregated data
1685  estm_tt = -999
1686  ! Intialise disaggregation method
1687  estmdisaggmethod = -999
1688 
1689  ! Generate useful sequences
1690  seq1nperestm = (/(i, i=1, nperestm, 1)/)
1691  seq1nsd = (/(i, i=1, nsd, 1)/)
1692 
1693  ! Get methods to use for disaggregation from RunControl
1694  ! (N.B.DisaggMethodESTM is set as 1 or 2 in RunControl; ESTMDisaggMethod is array of ncolsESTMdata used here)
1695  IF (disaggmethodestm == 1) THEN
1696  estmdisaggmethod(:) = 10 !linear disaggregation of averages
1697  ELSEIF (disaggmethodestm == 2) THEN
1698  estmdisaggmethod(:) = 20 !linear disaggregation of instantaneous values
1699  ELSE
1700  CALL errorhint(2, 'Problem in SUEWS_ESTMDisagg: DisaggMethodESTM value should be 1 or 2', &
1702  ENDIF
1703 
1704  ! Read data ---------------------------------------------------------------------
1705  IF (diagnosedisaggestm == 1) WRITE (*, *) 'Reading file: ', trim(fileorigestm)
1706  OPEN (lunit, file=trim(fileorigestm), status='old')
1707  ! CALL skipHeader(lunit,SkipHeaderMet) !Skip header -> read header instead
1708  READ (lunit, *) headerestm
1709  !write(*,*) HeaderMet
1710  ! Skip over lines that have already been read and downscaled
1711  IF (skippedlinesorigestm > 0) THEN
1712  DO i = 1, skippedlinesorigestm - 1 ! minus 1 here because last line of last block needs to be read again
1713  READ (lunit, *)
1714  ENDDO
1715  ! Read in last line of previous block
1716  READ (lunit, *, iostat=iostat_var) estmarrayorig
1717  estmfordisaggprev(1:ncolsestmdata) = estmarrayorig
1718  ENDIF
1719  ! Read in current block
1720  DO i = 1, readlinesorigestmdatamax
1721  READ (lunit, *, iostat=iostat_var) estmarrayorig
1722  estmfordisagg(i, 1:ncolsestmdata) = estmarrayorig
1723  ENDDO
1724  ! Read in first line of next block (except for last block)
1725  IF (iblock /= readblocksorigmetdata) THEN
1726  READ (lunit, *, iostat=iostat_var) estmarrayorig
1727  estmfordisaggnext(1:ncolsestmdata) = estmarrayorig
1728  ENDIF
1729  CLOSE (lunit)
1730 
1731  ! Check resolution of original met forcing data -------------------------------------
1732  ! Find time difference (in minutes) between first and second row
1733  tdiff = int(estmfordisagg(2, 4) - estmfordisagg(1, 4)) !Try using minutes
1734  IF (tdiff == 0) tdiff = int(estmfordisagg(2, 3) - estmfordisagg(1, 3))*60 !If no difference in minutes, try using hours
1735  IF (tdiff < 0) THEN !If time difference is negative (e.g. change of day), instead use second and third row
1736  tdiff = int(estmfordisagg(3, 4) - estmfordisagg(2, 4))
1737  IF (tdiff == 0) tdiff = int(estmfordisagg(3, 3) - estmfordisagg(2, 3))*60 !If no difference in minutes, try using hours
1738  ENDIF
1739  ! Check actual resolution matches specified input resolution
1740  IF (tdiff /= resolutionfilesinestm/60) THEN
1741  CALL errorhint(2, 'Problem in SUEWS_ESTMDisagg: timestamps in ESTM forcing file inconsistent with ResolutionFilesInESTM', &
1742  REAL(ResolutionFilesInESTM, KIND(1d0)), NotUsed, tdiff*60)
1743  ENDIF
1744 
1745  ! Disaggregate time columns ---------------------------------------------------------
1746  IF (diagnose == 1) THEN
1747  WRITE (*, *) 'Disaggregating ESTM forcing data (', trim(fileorigestm), ') to model time-step...'
1748  END IF
1749  CALL disaggregatedatetime(estmfordisagg(:, 1:4), tstep, nperestm, readlinesorigmetdatamax, estm_tt(:, 1:4))
1750  ! ! Convert to dectime
1751  ! dectimeOrig = ESTMForDisagg(:,2) + ESTMForDisagg(:,3)/24.0 + ESTMForDisagg(:,4)/(60.0*24.0)
1752  !
1753  ! DO i=1,ReadLinesOrigESTMDataMax
1754  ! ! Downscale dectime using dectimeOrig(i) [becomes timestamp of last subinterval]
1755  ! dectimeDscd(NperESTM*(i-1)+Seq1NperESTM) = dectimeOrig(i) - (tstep/60.0)/(60.0*24.0)*(/(ii, ii=(NperESTM-1),0, -1)/)
1756  ! ! Convert to required formats
1757  ! temp_iy = INT(ESTMForDisagg(i,1)) !Copy year
1758  ! temp_id = FLOOR(dectimeDscd(NperESTM*(i-1)+Seq1NperESTM)) !DOY
1759  ! ! To avoid precision errors, round here
1760  ! ! - this should not be a problem as a difference of 1 = 1 min, so a difference of 0.001 << 1 min
1761  ! temp_ihm = NINT(((dectimeDscd(NperESTM*(i-1)+Seq1NperESTM) - temp_id/1.0)*60.0*24.0)*1000.0)/1000 !Minutes of the day (1440 max)
1762  ! temp_ih = (temp_ihm-MOD(temp_ihm,60))/60 !Hours
1763  ! temp_im = MOD(temp_ihm,60) !Minutes
1764  !
1765  ! IF(dectimeOrig(i) == 1.0000 .AND. i > 1) THEN !If year changes and it is not the beginning of the dataset
1766  ! WRITE(*,*) 'Year change encountered: ', dectimeOrig(i), dectimeOrig(i-1)
1767  ! ! Re-downscale dectime using dectimeOrig(i-1)
1768  ! dectimeDscd(NperESTM*(i-1)+Seq1NperESTM) = dectimeOrig(i-1) + (tstep/60.0)/(60.0*24.0)*Seq1NperESTM
1769  ! ! Convert to required formats
1770  ! temp_iy = INT(ESTMForDisagg(i,1)) !Copy year
1771  ! temp_id = FLOOR(dectimeDscd(NperESTM*(i-1)+Seq1NperESTM)) !DOY
1772  ! temp_ihm = NINT(((dectimeDscd(NperESTM*(i-1)+Seq1NperESTM) - temp_id/1.0)*60.0*24.0)*1000.0)/1000 !Mins of the day (1440 max)
1773  ! temp_ih = (temp_ihm-MOD(temp_ihm,60))/60 !Hours
1774  ! temp_im = MOD(temp_ihm,60) !Minutes
1775  ! ! Adjust year and DOY to account for year change
1776  ! temp_iy(1:(NperESTM-1)) = temp_iy(1:(NperESTM-1)) - 1 !Subtract 1 from year for all except final timestamp
1777  ! temp_id(NperESTM) = 1 !Set day for final timestamp to 1
1778  ! ENDIF
1779  !
1780  ! !IF(i==1 .or. i== ReadlinesOrigESTMDataMax) THEN
1781  ! ! write(*,*) temp_iy
1782  ! ! write(*,*) temp_id
1783  ! ! !write(*,*) temp_ihm
1784  ! ! write(*,*) temp_ih
1785  ! ! write(*,*) temp_im
1786  ! !ENDIF
1787  !
1788  ! ! Copy to ESTM_tt array
1789  ! ESTM_tt(NperESTM*(i-1)+Seq1NperESTM,1) = temp_iy
1790  ! ESTM_tt(NperESTM*(i-1)+Seq1NperESTM,2) = temp_id
1791  ! ESTM_tt(NperESTM*(i-1)+Seq1NperESTM,3) = temp_ih
1792  ! ESTM_tt(NperESTM*(i-1)+Seq1NperESTM,4) = temp_im
1793  !
1794  ! ENDDO
1795 
1796  ! Disaggregate other columns --------------------------------------------------------
1797  ! All other columns are temperatures
1798  DO ii = 5, ncolsestmdata
1799  IF (all(estmfordisagg(:, ii) == -999)) THEN
1800  !IF(DiagnoseDisaggESTM==1) write(*,*) 'No data for col.', ii
1801  estm_tt(:, ii) = -999
1802  ELSE
1803  estm_tt(:, ii) = disagg_lin(estmfordisagg(:, ii), estmfordisaggprev(ii), estmfordisaggnext(ii), estmdisaggmethod(ii), &
1805  ENDIF
1806  ENDDO
1807 
1808  ! Copy disaggregated data to MetForcingDataArray
1810 
1811  ! Write out disaggregated file ------------------------------------------------------
1812  IF (keeptstepfilesin == 1) THEN
1813  IF (iblock == 1) THEN
1814  ! Prepare header
1815  DO i = 1, ncolsestmdata
1816  IF (i == 1) THEN
1817  headerestmout = adjustl(headerestm(i))
1818  ELSE
1819  headerestmout = trim(headerestmout)//' '//adjustl(headerestm(i))
1820  ENDIF
1821  ENDDO
1822  ! Write out header
1823  OPEN (78, file=trim(filedscdestm), err=113)
1824  WRITE (78, '(a)') headerestmout
1825  ELSE
1826  OPEN (78, file=trim(filedscdestm), position='append')!,err=113)
1827  ENDIF
1828  ! Write out data
1830  WRITE (78, 304) (int(estm_tt(i, ii)), ii=1, 4), estm_tt(i, 5:ncolsestmdata)
1831  ENDDO
1832  !IF(iBlock == ReadBlocksOrigMetData) THEN
1833  ! WRITE(78,'(i2)') -9
1834  ! WRITE(78,'(i2)') -9
1835  !ENDIF
1836  CLOSE (78) !Close output file
1837  ENDIF
1838 
1839 304 FORMAT((i4, 1x), 3(i3, 1x), 9(f9.4, 1x))
1840 
1841  ! Deallocate arrays -----------------------------------------------------------------
1842  DEALLOCATE (estmfordisagg)
1843  DEALLOCATE (estmfordisaggprev)
1844  DEALLOCATE (estmfordisaggnext)
1845 
1846  RETURN
1847 
1848 113 CALL errorhint(52, trim(filedscdestm), notused, notused, notusedi)
1849 
1850  ENDSUBROUTINE disaggregateestm
1851  !======================================================================================
1852 
1853  !======================================================================================
1854  SUBROUTINE disaggregatedatetime(DateTimeForDisagg, tstep, Nper, ReadLinesOrigMetDataMax, DateTimeDscd)
1855  USE datetime_module, ONLY: daysinyear
1856  IMPLICIT NONE
1857  INTEGER, INTENT(in) :: tstep, Nper, ReadLinesOrigMetDataMax
1858  REAL(KIND(1d0)), DIMENSION(ReadLinesOrigMetData, 4), INTENT(in):: DateTimeForDisagg
1859  REAL(KIND(1d0)), DIMENSION(ReadLinesOrigMetData*Nper, 4), INTENT(out):: DateTimeDscd
1860 
1861  REAL(KIND(1d0)), DIMENSION(ReadLinesOrigMetData):: dectimeOrig
1862  REAL(KIND(1d0)), DIMENSION(Nper) :: temp_dectime !temorary varaibles for disaggragation
1863  INTEGER, DIMENSION(Nper):: temp_iy, temp_id, temp_ih, temp_im, temp_ihm ! temorary varaibles for disaggragation
1864  INTEGER, DIMENSION(Nper)::ndays_iy ! number of days in iy
1865 
1866  INTEGER :: i, ii
1867  INTEGER, DIMENSION(Nper):: seq1Nper
1868 
1869  ! Generate useful sequences
1870  seq1nper = (/(i, i=1, nper, 1)/)
1871  ! Convert to dectime
1872  ! dectimeOrig = MetForDisagg(:,2) + MetForDisagg(:,3)/24.0 + MetForDisagg(:,4)/(60.0*24.0)
1873  ! correct to dectime(year_start)=0 and dectime(year_end)=days of year (i.e., 365 or 366 if leap year) ! TS 09 May 2018
1874  dectimeorig = (datetimefordisagg(:, 2) - 1) + datetimefordisagg(:, 3)/24.0 + datetimefordisagg(:, 4)/(60.0*24.0)
1875 
1876  DO i = 1, readlinesorigmetdatamax
1877 
1878  ! Downscale dectime using dectimeOrig(i) [becomes timestamp of last subinterval]
1879  temp_dectime = dectimeorig(i) - (tstep/60.0)/(60.0*24.0)*(/(ii, ii=(nper - 1), 0, -1)/)
1880  temp_dectime = nint(temp_dectime*60*60*24)/(60*60*24*1.)
1881 
1882  ! Convert to required formats
1883  ! get year: year-1 if dectime <0, copy `year` from metforcing otherwise
1884  temp_iy = merge(int(datetimefordisagg(i, 1)) - 1, int(datetimefordisagg(i, 1)), temp_dectime < 0)
1885 
1886  ! get day of year:
1887  ndays_iy = daysinyear(temp_iy) ! get days of year for DOY correction when temp_dectime<0 during year-crossing
1888  temp_dectime = merge(temp_dectime + ndays_iy, temp_dectime, temp_dectime < 0) ! correct minus temp_dectime to positive values
1889  temp_id = floor(temp_dectime) + 1 !DOY
1890 
1891  temp_ihm = nint((temp_dectime + 1 - temp_id/1.0)*60.0*24.0) !Minutes of the day (1440 max)
1892  temp_ih = (temp_ihm - mod(temp_ihm, 60))/60 !Hours
1893  temp_ih = merge(temp_ih, 0, mask=(temp_ih < 24))
1894  temp_im = mod(temp_ihm, 60) !Minutes
1895 
1896  ! Copy to Met_tt array
1897  datetimedscd(nper*(i - 1) + seq1nper, 1) = temp_iy
1898  datetimedscd(nper*(i - 1) + seq1nper, 2) = temp_id
1899  datetimedscd(nper*(i - 1) + seq1nper, 3) = temp_ih
1900  datetimedscd(nper*(i - 1) + seq1nper, 4) = temp_im
1901 
1902  ENDDO
1903 
1904  END SUBROUTINE disaggregatedatetime
1905  !======================================================================================
1906 
1907  ! Define functions here:
1908  FUNCTION disagg_lin(Slow, SlowPrev, SlowNext, DisaggType, Nper_loc, ReadLinesOrig_loc, ReadLinesOrigMax_loc, iBlock) RESULT(Fast)
1910  USE defaultnotused
1911  USE sues_data
1912 
1913  IMPLICIT NONE
1914 
1915  INTEGER:: DisaggType !Type of disaggregation: 10 for averaged variables; 20 for instantaneous variables
1916  INTEGER:: Nper_loc !Number of subintervals per interval (local Nper)
1917  INTEGER:: ReadLinesOrig_loc, ReadLinesOrigMax_loc !Number of lines to read in original file (local)
1918  INTEGER:: iBlock
1919  REAL(KIND(1d0)), DIMENSION(ReadLinesOrig_loc*Nper_loc):: Fast !Array to receive disaggregated data
1920  REAL(KIND(1d0)), DIMENSION(ReadLinesOrig_loc):: Slow !Array to disaggregate
1921  REAL(KIND(1d0)):: SlowPrev, SlowNext
1922  INTEGER, DIMENSION(Nper_loc):: FastRows !Group of rows that are filled with each iteration
1923  INTEGER, DIMENSION(FLOOR(Nper_loc/2.0)):: FirstRows10 !Rows at the beginning that are not filled during iteration (for averages)
1924  INTEGER, DIMENSION(Nper_loc - FLOOR(Nper_loc/2.0)):: LastRows10 !Rows at the end that are not filled during iteration
1925  INTEGER, DIMENSION(Nper_loc):: FirstRows20 !Rows at the beginning that are not filled during iteration (for instantaneous)
1926  INTEGER, DIMENSION(Nper_loc):: seq1Nper_loc !1 to Nper_loc
1927  INTEGER:: XNper_loc !XNper_loc = 2 for even Nper_loc; XNper_loc=1 for odd Nper_loc
1928  INTEGER:: i, ii !counters
1929 
1930  ! Calculate XNper_loc (differentiates between disaggregations with odd and even Nper_loc)
1931  IF (mod(nper_loc, 2) == 0) xnper_loc = 2
1932  IF (mod(nper_loc, 2) == 1) xnper_loc = 1
1933 
1934  seq1nper_loc = (/(i, i=1, nper_loc, 1)/)
1935 
1936  ! Setup counters for iteration
1937  IF (disaggtype == 10) THEN
1938  fastrows = floor(nper_loc/2.0) + seq1nper_loc ! Rows to create at model time-step
1939  firstrows10 = (/(i, i=1, (fastrows(1) - 1), 1)/) !For start of dataset
1940  lastrows10 = &
1941  (/(i, i=nper_loc*(readlinesorigmax_loc - 1 - 1) + fastrows(nper_loc) + 1, &
1942  (readlinesorigmax_loc*nper_loc), 1)/) ! For end of dataset
1943  ELSEIF (disaggtype == 20) THEN
1944  fastrows = nper_loc + seq1nper_loc !Rows to create at model time-step
1945  firstrows20 = (/(i, i=1, (fastrows(1) - 1), 1)/) !For start of dataset
1946  ENDIF
1947 
1948  ! Initialise fast array to -999
1949  fast = -999
1950  ! Linearly disaggregate
1951  IF (disaggtype == 10) THEN !Averaged variables
1952  IF (diagnosedisagg == 1) WRITE (*, *) 'Linearly disaggregating averaged variable'
1953  DO i = 1, (readlinesorigmax_loc - 1)
1954  fast(nper_loc*(i - 1) + fastrows) = slow(i) - &
1955  (slow(i + 1) - slow(i))/(xnper_loc*nper_loc) + &
1956  (slow(i + 1) - slow(i))/nper_loc*(/(ii, ii=1, nper_loc, 1)/)
1957  ENDDO
1958 
1959  ! For first few rows, use previous met block
1960  IF (iblock == 1) THEN
1961  fast(firstrows10) = fast(fastrows(1)) !Use repeat values at the start of the year
1962  ELSE
1963  fast(firstrows10) = slowprev - &
1964  (slow(1) - slowprev)/(xnper_loc*nper_loc) + &
1965  (slow(1) - slowprev)/nper_loc* &
1966  (/(ii, ii=(nper_loc - SIZE(firstrows10) + 1), nper_loc, 1)/)
1967  ENDIF
1968  ! For last few rows, use next met block
1969  IF (iblock == readblocksorigmetdata) THEN
1970  fast(lastrows10) = fast(nper_loc*(readlinesorigmax_loc - 1 - 1) + fastrows(nper_loc)) !Use repeat values at the end of the year
1971  ELSE
1972  fast(lastrows10) = slow(readlinesorigmax_loc) - &
1973  (slownext - slow(readlinesorigmax_loc))/(xnper_loc*nper_loc) + &
1974  (slownext - slow(readlinesorigmax_loc))/nper_loc* &
1975  (/(ii, ii=1, SIZE(lastrows10), 1)/)
1976  ENDIF
1977  ELSEIF (disaggtype == 20) THEN !Instantaneous variables
1978  IF (diagnosedisagg == 1) WRITE (*, *) 'Linearly disaggregating instantaneous variable'
1979  DO i = 1, (readlinesorigmax_loc - 1)
1980  fast(nper_loc*(i - 1) + fastrows) = (slow(i) + &
1981  (slow(i + 1) - slow(i))/nper_loc*2*(seq1nper_loc - 1) + &
1982  slow(i))/2
1983  ENDDO
1984  ! For first few rows, use previous met block
1985  IF (iblock == 1) THEN
1986  fast(firstrows20) = fast(fastrows(1)) !Use repeat values at the start of the year
1987  ELSE
1988  fast(firstrows20) = (slowprev + &
1989  (slow(1) - slowprev)/nper_loc*2* &
1990  ((/(ii, ii=(nper_loc - SIZE(firstrows20) + 1), nper_loc, 1)/) - 1) + &
1991  slowprev)/2
1992  ENDIF
1993  !! Last few rows are already filled for the instantaneous value disaggregation
1994  !IF(iBlock==ReadBlocksOrigMetData) THEN
1995  ! Fast(LastRows20) = Fast(Nper_loc*(ReadLinesOrigMax_loc-1-1)+FastRows(Nper_loc)) !Use repeat values at the end of the year
1996  !ELSE
1997  ! Fast(LastRows20) = (Slow(ReadLinesOrigMax_loc) + &
1998  ! (SlowNext-Slow(ReadLinesOrigMax_loc))/Nper_loc*2 * &
1999  ! ((/(ii, ii=1,SIZE(LastRows20), 1)/)-1) + &
2000  ! Slow(ReadLinesOrigMax_loc))/2
2001  !ENDIF
2002  ENDIF
2003 
2004  IF (any(fast(1:readlinesorigmax_loc*nper_loc) == -999)) THEN
2005  WRITE (*, *) 'Problem: -999s (', count(fast(1:readlinesorigmax_loc*nper_loc) == -999), ') in disaggregated data.'
2006  CALL errorhint(13, 'Problem in SUEWS_MetDisagg: -999 values in disaggregated data.', notused, notused, notusedi)
2007  ENDIF
2008 
2009  ENDFUNCTION disagg_lin
2010  !======================================================================================
2011 
2012  !======================================================================================
2013  FUNCTION disaggp_amongn(Slow, amongN, Nper_loc, ReadLinesOrig_loc, ReadLinesOrigMax_loc) RESULT(Fast)
2014  ! Subroutine to disaggregate precipitation by evenly distributing among N subintervals
2015  ! (i.e. equal intensity in N subintervals)
2016  ! See Ward et al. (in review), meanN, 0.5N or 0.25N approach
2017  ! HCW 10 Feb 2017
2018  !======================================================================================
2019 
2020  USE defaultnotused
2021  USE sues_data
2022 
2023  IMPLICIT NONE
2024 
2025  INTEGER:: amongN !Number of subintervals over which rain will be distributed
2026  INTEGER:: Nper_loc !Number of subintervals per interval (local Nper)
2027  INTEGER:: ReadLinesOrig_loc, ReadLinesOrigMax_loc !Number of lines to read in original file (local)
2028  REAL(KIND(1d0)), DIMENSION(ReadLinesOrig_loc*Nper_loc):: Fast !Array to receive disaggregated data
2029  REAL(KIND(1d0)), DIMENSION(ReadLinesOrig_loc):: Slow !Array to disaggregate
2030  INTEGER, DIMENSION(:), ALLOCATABLE:: Subintervals !Array of subintervals that contain rain
2031  INTEGER, DIMENSION(Nper_loc):: seq1Nper_loc !1 to Nper_loc
2032  INTEGER:: i
2033 
2034  ! For each averaging period, get subintervals which will receive rain
2035  ALLOCATE (subintervals(amongn))
2036  subintervals(:) = -999
2037 
2038  seq1nper_loc = (/(i, i=1, nper_loc, 1)/)
2039 
2040  IF (diagnosedisagg == 1) WRITE (*, *) 'Distributing over ', amongn, ' subintervals for variable'
2041  ! If all subintervals are to contain rain, don't need to generate random numbers
2042  IF (amongn == nper_loc) THEN
2043  subintervals(:) = seq1nper_loc
2044  ENDIF
2045  IF (amongn > nper_loc) &
2046  CALL errorhint(2, 'Problem in SUEWS_MetDisagg: no. of rainy periods cannot exceed number of subintervals', &
2047  REAL(Nper_loc, KIND(1d0)), NotUsed, amongN)
2048 
2049  ! Initialise fast array to -999
2050  fast = -999
2051  DO i = 1, readlinesorigmax_loc
2052  fast(nper_loc*(i - 1) + seq1nper_loc) = 0 !Fill all subintervals with zeros initially
2053  IF (slow(i) > 0) THEN !If there is some rainfall during this interval...
2054  IF (amongn < nper_loc) THEN
2055  subintervals(:) = -999
2056  subintervals = randomsamples(amongn, nper_loc)
2057  ENDIF
2058  fast(nper_loc*(i - 1) + subintervals) = slow(i)/amongn
2059  ENDIF
2060  ENDDO
2061 
2062  IF (any(fast(1:readlinesorigmax_loc*nper_loc) == -999)) THEN
2063  WRITE (*, *) 'Problem: -999s (', count(fast(1:readlinesorigmax_loc*nper_loc) == -999), ') in disaggregated data'
2064  CALL errorhint(13, 'Problem in SUEWS_MetDisagg: -999 values in disaggregated data.', notused, notused, notusedi)
2065  ENDIF
2066 
2067  ENDFUNCTION disaggp_amongn
2068  !======================================================================================
2069 
2070  !======================================================================================
2071  FUNCTION disaggp_amongnmult(Slow, multupperI, multamongN, Nper_loc, ReadLinesOrig_loc, ReadLinesOrigMax_loc) RESULT(Fast)
2072  ! Subroutine to disaggregate precipitation by evenly distributing among N subintervals
2073  ! (i.e. equal intensity in N subintervals) for different intensity bins
2074  ! Based on analsysis by Wen Gu
2075  ! HCW 21 Apr 2017
2076  !======================================================================================
2077 
2078  USE defaultnotused
2079  USE sues_data
2080 
2081  IMPLICIT NONE
2082 
2083  REAL(KIND(1d0)), DIMENSION(5):: multupperI !Upper bound of intensity bin
2084  INTEGER, DIMENSION(5):: multamongN !Number of subintervals over which rain will be distributed (array)
2085  INTEGER:: thisamongN !Number of subintervals over which rain will be distributed
2086  INTEGER:: Nper_loc !Number of subintervals per interval (local Nper)
2087  INTEGER:: ReadLinesOrig_loc, ReadLinesOrigMax_loc !Number of lines to read in original file (local)
2088  REAL(KIND(1d0)), DIMENSION(ReadLinesOrig_loc*Nper_loc):: Fast !Array to receive disaggregated data
2089  REAL(KIND(1d0)), DIMENSION(ReadLinesOrig_loc):: Slow !Array to disaggregate
2090  INTEGER, DIMENSION(:), ALLOCATABLE:: Subintervals !Array of subintervals that contain rain
2091  INTEGER, DIMENSION(Nper_loc):: seq1Nper_loc !1 to Nper_loc
2092  INTEGER:: i
2093 
2094  seq1nper_loc = (/(i, i=1, nper_loc, 1)/)
2095 
2096  IF (diagnosedisagg == 1) WRITE (*, *) 'Distributing over variable subintervals depending on intensity for variable'
2097 
2098  ! Initialise fast array to -999
2099  fast = -999
2100  DO i = 1, readlinesorigmax_loc
2101  fast(nper_loc*(i - 1) + seq1nper_loc) = 0 !Fill all subintervals with zeros initially
2102  IF (slow(i) > 0) THEN !If there is some rainfall during this interval...
2103  !Use intensity in this interval to decide number of subintervals to fill with rain
2104  IF (slow(i) <= multupperi(1)) THEN
2105  thisamongn = multamongn(1)
2106  ELSEIF (slow(i) > multupperi(1) .AND. slow(i) <= multupperi(2)) THEN
2107  thisamongn = multamongn(2)
2108  ELSEIF (slow(i) > multupperi(2) .AND. slow(i) <= multupperi(3)) THEN
2109  thisamongn = multamongn(3)
2110  ELSEIF (slow(i) > multupperi(3) .AND. slow(i) <= multupperi(4)) THEN
2111  thisamongn = multamongn(4)
2112  ELSEIF (slow(i) > multupperi(4) .AND. slow(i) <= multupperi(5)) THEN
2113  thisamongn = multamongn(5)
2114  ELSEIF (slow(i) > multupperi(5)) THEN
2115  thisamongn = multamongn(5)
2116  CALL errorhint(4, 'Precip in met forcing file exceeds maxiumum MultRainAmongNUpperI', &
2117  slow(i), multrainamongnupperi(5), notusedi)
2118  ENDIF
2119 
2120  ! For each averaging period, get subintervals which will receive rain
2121  ALLOCATE (subintervals(thisamongn))
2122  subintervals(:) = -999
2123 
2124  IF (thisamongn > nper_loc) CALL errorhint(2, 'Problem in SUEWS_MetDisagg: no. of rainy periods cannot exceed '// &
2125  'number of subintervals', REAL(Nper_loc, KIND(1d0)), NotUsed, thisamongN)
2126 
2127  IF (thisamongn == nper_loc) THEN ! If all subintervals are to contain rain, don't need to generate random numbers
2128  subintervals(:) = seq1nper_loc
2129  ELSEIF (thisamongn < nper_loc) THEN
2130  subintervals = randomsamples(thisamongn, nper_loc)
2131  ENDIF
2132  fast(nper_loc*(i - 1) + subintervals) = slow(i)/thisamongn
2133  !write(*,*) Slow(i), thisamongN
2134  DEALLOCATE (subintervals)
2135  ENDIF
2136  ENDDO
2137 
2138  IF (any(fast(1:readlinesorigmax_loc*nper_loc) == -999)) THEN
2139  WRITE (*, *) 'Problem: -999s (', count(fast(1:readlinesorigmax_loc*nper_loc) == -999), ') in disaggregated data'
2140  CALL errorhint(13, 'Problem in SUEWS_MetDisagg: -999 values in disaggregated data.', notused, notused, notusedi)
2141  ENDIF
2142 
2143  ENDFUNCTION disaggp_amongnmult
2144  !======================================================================================
2145 
2146  !======================================================================================
2147  FUNCTION randomsamples(N, OutOf) RESULT(Samples)
2148  ! Generates N/OutOf random samples without repeats
2149  ! e.g. for N = 3 and OutOf = 12, a possibility for Samples = 7,3,11
2150  ! HCW 10 Feb 2017
2151  !======================================================================================
2152 
2153  IMPLICIT NONE
2154 
2155  INTEGER:: i !counter
2156  INTEGER:: N !number of samples to return
2157  INTEGER:: OutOf !number to sample from
2158  INTEGER:: X !next sample to be added
2159  REAL(KIND(1D0)):: r !random number
2160  INTEGER, DIMENSION(:), ALLOCATABLE:: Samples !Array to receive random samples
2161 
2162  ! Allocate and initialise Samples
2163  ALLOCATE (samples(n))
2164  samples(:) = -999
2165 
2166  ! Generate random sample (no repeats)
2167  i = 0 !Set counter to zero initially
2168  DO WHILE (any(samples == -999))
2169  CALL random_number(r)
2170  x = int(r*outof) + 1
2171  !write(*,*) X
2172  !write(*,*) COUNT(Samples == X)
2173  IF (count(samples == x) == 0) THEN
2174  ! Only keep if this subinterval has not already been selected
2175  i = i + 1
2176  samples(i) = x
2177  ENDIF
2178  !write(*,*) Samples
2179  ENDDO
2180 
2181  ENDFUNCTION randomsamples
2182  !======================================================================================
2183 
2184 ENDMODULE metdisagg
2185 !========================================================================================
integer readlinesorigmetdata
integer inputmetformat
subroutine codematchdist(rr, CodeCol, codeColSameSurf)
character(len=20), dimension(ncolumnsbiogen) headerbiogen_reqd
integer disaggmethodestm
integer, dimension(nsurfincsnow) c_ohmcode_swet
real(kind(1d0)), dimension(:), allocatable metfordisaggprev
character(len=20), dimension(ncolumnssoil) headersoil_reqd
integer diagnosedisaggestm
integer skippedlinesorig
real(kind(1d0)) halftimestep
real(kind(1d0)), dimension(:, :), allocatable water_coeff
subroutine codematchprof(Gridiv, SurfaceCharCodeCol)
subroutine inputheadercheck(FileName)
subroutine codematchnonveg(rr, CodeCol)
integer nlinessoil
integer, parameter ncolumnsmetforcingdata
real(kind(1d0)), dimension(:, :), allocatable irrigation_coeff
integer nlinesestmcoefficients
subroutine narp_cal_sunposition(year, idectime, UTC, locationlatitude, locationlongitude, locationaltitude, sunazimuth, sunzenith)
integer nlinesprofiles
real(kind(1d0)), dimension(:, :), allocatable soil_coeff
real(kind(1d0)) zenith_deg
real(kind(1d0)) function get_prof_spectime_sum(Hour, Min, Sec, Prof_24h, dt)
integer, dimension(5) c_code_estmclass_bldgs
subroutine codematchbiogen(Gridiv, SurfaceCharCodeCol)
integer diagnose
real(kind(1d0)) nan
real(kind(1d0)) alt
subroutine disaggregatemet(iBlock, igrid)
integer keeptstepfilesin
integer readlinesorigestmdatamax
real(kind(1d0)), dimension(:, :), allocatable estmcoefficients_coeff
character(len=150) fileorigmet
real(kind(1d0)), dimension(:, :), allocatable snow_coeff
real(kind(1d0)) notused
real(kind(1d0)), dimension(:), allocatable estmfordisaggnext
character(len=20), dimension(ncolumnswgwaterdist) headerwgwaterdist_file
subroutine codematchconductance(rr, CodeCol)
integer, dimension(nsurfincsnow) c_ohmcode_wwet
integer, dimension(24) cpr_hours
character(len=20), dimension(ncolumnswgwaterdist) headerwgwaterdist_reqd
real(kind(1d0)) soilrocks
integer snowuse
integer nlinesirrigation
integer, dimension(nsurfincsnow) c_estmcode
subroutine skipheader(lfn, skip)
integer, dimension(nsurfincsnow) c_ohmcode_wdry
subroutine metread(lfn, MetArray, InputmetFormat, ldown_option, NetRadiationMethod, snowUse, SMDMethod, SoilDepthMeas, SoilRocks, SoilDensity, SmCap)
real(kind(1d0)), dimension(5) multrainamongnupperi
subroutine codematchestm(Gridiv, is)
character(len=90), dimension(14) text
real(kind(1d0)), dimension(:, :), allocatable surfacechar
integer nlinesveg
character(len=20), dimension(ncolumnsprofiles) headerprofiles_file
character(len=20), dimension(ncolumnsirrigation) headerirrigation_file
subroutine codematchveg(rr, CodeCol)
integer function, dimension(:), allocatable randomsamples(N, OutOf)
character(len=20), dimension(ncolumnsanthropogenic) headeranthropogenic_file
integer nlinesanthropogenic
integer resolutionfilesinestm
subroutine disaggregateestm(iBlock)
real(kind(1d0)) soildensity
real(kind(1d0)), dimension(:), allocatable estmfordisaggprev
integer readlinesorigmetdatamax
integer nlinesconductance
character(len=150) filedscdestm
integer kdownzen
real(kind(1d0)) lng
integer nlinesohmcoefficients
character(len=150) filedscdmet
character(len=20), dimension(ncolumnswater) headerwater_reqd
subroutine codematchohm(Gridiv, is, SWWD)
real(kind(1d0)) soildepthmeas
integer nlineswater
character(len=20), dimension(ncolumnssoil) headersoil_file
subroutine codematchsoil(Gridiv, SurfaceCharCodeCol)
real(kind(1d0)) lat
real(kind(1d0)), dimension(:, :), allocatable siteselect
real(kind(1d0)), dimension(:, :), allocatable ohmcoefficients_coeff
subroutine codematchanthropogenic(rr, CodeCol)
real(kind(1d0)), dimension(:, :), allocatable anthropogenic_coeff
real(kind(1d0)) smcap
real(kind(1d0)) function, dimension(readlinesorig_loc *nper_loc) disagg_lin(Slow, SlowPrev, SlowNext, DisaggType, Nper_loc, ReadLinesOrig_loc, ReadLinesOrigMax_loc, iBlock)
real(kind(1d0)), dimension(:, :), allocatable estmfordisagg
subroutine codematchsnow(rr, CodeCol)
real(kind(1d0)), dimension(:), allocatable metfordisaggnext
character(len=20), dimension(ncolumnsbiogen) headerbiogen_file
integer, parameter ncolsestmdata
real(kind(1d0)), dimension(:, :), allocatable biogen_coeff
integer raindisaggmethod
integer nlinesnonveg
character(len=20), dimension(ncolumnssnow) headersnow_file
integer readblocksorigmetdata
real(kind(1d0)), dimension(:, :), allocatable wgwaterdist_coeff
real(kind(1d0)) function, dimension(readlinesorig_loc *nper_loc) disaggp_amongnmult(Slow, multupperI, multamongN, Nper_loc, ReadLinesOrig_loc, ReadLinesOrigMax_loc)
integer nlinesbiogen
character(len=20), dimension(ncolumnsnonveg) headernonveg_file
integer disaggmethod
real(kind(1d0)) azimuth
real(kind(1d0)), dimension(:, :), allocatable veg_coeff
character(len=20), dimension(ncolumnssnow) headersnow_reqd
integer, dimension(3) c_code_estmclass_paved
character(len=20), dimension(ncolumnsohmcoefficients) headerohmcoefficients_file
character(len=20), dimension(ncolumnswater) headerwater_file
integer ldown_option
real(kind(1d0)) function get_prof_spectime_inst(Hour, Min, Sec, Prof_24h)
character(len=20), dimension(ncolumnsveg) headerveg_file
subroutine run_control(eval, LowerLimit, Upperlimit)
real(kind(1d0)) reall
real(kind(1d0)) function, dimension(readlinesorig_loc *nper_loc) disaggp_amongn(Slow, amongN, Nper_loc, ReadLinesOrig_loc, ReadLinesOrigMax_loc)
character(len=20), dimension(ncolumnsprofiles) headerprofiles_reqd
integer nlinessnow
integer smdmethod
integer nlineswgwaterdist
integer readlinesorigestmdata
integer, dimension(nsurfincsnow) c_ohmcode_sdry
character(len=20), dimension(ncolumnsirrigation) headerirrigation_reqd
character(len=150) fileorigestm
integer resolutionfilesin
subroutine codematchestm_class(Gridiv, is, ii)
integer netradiationmethod
subroutine disaggregatedatetime(DateTimeForDisagg, tstep, Nper, ReadLinesOrigMetDataMax, DateTimeDscd)
real(kind(1d0)), dimension(:, :), allocatable nonveg_coeff
real(kind(1d0)), dimension(:, :, :), allocatable metforcingdata
integer, dimension(5) multrainamongn
character(len=20), dimension(ncolumnsveg) headerveg_reqd
character(len=20), dimension(ncolumnsestmcoefficients) headerestmcoefficients_reqd
integer rainamongn
real(kind(1d0)), dimension(:, :, :), allocatable estmforcingdata
real(kind(1d0)), dimension(:, :), allocatable conductance_coeff
character(len=20), dimension(ncolumnsohmcoefficients) headerohmcoefficients_reqd
integer diagnosedisagg
character(len=20), dimension(ncolumnsanthropogenic) headeranthropogenic_reqd
integer, parameter pavsurf
character(len=20), dimension(ncolumnsestmcoefficients) headerestmcoefficients_file
subroutine errorhint(errh, ProblemFile, VALUE, value2, valueI)
integer, parameter bldgsurf
character(len=20), dimension(ncolumnsconductance) headercond_file
real(kind(1d0)) function get_prof_spectime_mean(Hour, Min, Sec, Prof_24h)
character(len=20), dimension(ncolumnsconductance) headercond_reqd
real(kind(1d0)), dimension(:, :), allocatable profiles_coeff
real(kind(1d0)), dimension(:, :), allocatable metfordisagg
subroutine codematchwater(rr, CodeCol)
integer gridcounter
subroutine codematchirrigation(rr, CodeCol)
integer skippedlinesorigestm
real(kind(1d0)) timezone
character(len=20), dimension(ncolumnsnonveg) headernonveg_reqd