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"
534  headerirrigation_reqd(cir_daywat1) = "DayWat(1)"
535  headerirrigation_reqd(cir_daywat2) = "DayWat(2)"
536  headerirrigation_reqd(cir_daywat3) = "DayWat(3)"
537  headerirrigation_reqd(cir_daywat4) = "DayWat(4)"
538  headerirrigation_reqd(cir_daywat5) = "DayWat(5)"
539  headerirrigation_reqd(cir_daywat6) = "DayWat(6)"
540  headerirrigation_reqd(cir_daywat7) = "DayWat(7)"
541  headerirrigation_reqd(cir_daywatper1) = "DayWatPer(1)"
542  headerirrigation_reqd(cir_daywatper2) = "DayWatPer(2)"
543  headerirrigation_reqd(cir_daywatper3) = "DayWatPer(3)"
544  headerirrigation_reqd(cir_daywatper4) = "DayWatPer(4)"
545  headerirrigation_reqd(cir_daywatper5) = "DayWatPer(5)"
546  headerirrigation_reqd(cir_daywatper6) = "DayWatPer(6)"
547  headerirrigation_reqd(cir_daywatper7) = "DayWatPer(7)"
548 
549  ! ========== SUEWS_Profiles.txt ===============
550  headerprofiles_reqd(cpr_code) = "Code"
560  headerprofiles_reqd(cpr_hours(10)) = "9"
561  headerprofiles_reqd(cpr_hours(11)) = "10"
562  headerprofiles_reqd(cpr_hours(12)) = "11"
563  headerprofiles_reqd(cpr_hours(13)) = "12"
564  headerprofiles_reqd(cpr_hours(14)) = "13"
565  headerprofiles_reqd(cpr_hours(15)) = "14"
566  headerprofiles_reqd(cpr_hours(16)) = "15"
567  headerprofiles_reqd(cpr_hours(17)) = "16"
568  headerprofiles_reqd(cpr_hours(18)) = "17"
569  headerprofiles_reqd(cpr_hours(19)) = "18"
570  headerprofiles_reqd(cpr_hours(20)) = "19"
571  headerprofiles_reqd(cpr_hours(21)) = "20"
572  headerprofiles_reqd(cpr_hours(22)) = "21"
573  headerprofiles_reqd(cpr_hours(23)) = "22"
574  headerprofiles_reqd(cpr_hours(24)) = "23"
575 
576  ! ========== SUEWS_WithinGridWaterDist.txt ====
586  headerwgwaterdist_reqd(cwg_tosoilstore) = "ToSoilStore"
587 
588  ! ========== SUEWS_BiogenCO2.txt ======
589  headerbiogen_reqd(cb_code) = "Code"
590  headerbiogen_reqd(cb_alpha) = "alpha"
591  headerbiogen_reqd(cb_beta) = "beta"
592  headerbiogen_reqd(cb_theta) = "theta"
593  headerbiogen_reqd(cb_alpha_enh) = "alpha_enh"
594  headerbiogen_reqd(cb_beta_enh) = "beta_enh"
595  headerbiogen_reqd(cb_resp_a) = "resp_a"
596  headerbiogen_reqd(cb_resp_b) = "resp_b"
597  headerbiogen_reqd(cb_min_r) = "min_respi"
598 
599  ! =======================================================
600 
601  !write(*,*) 'Checking header for ', FileName
602  ! Check columns in input files match model code
603 
604  IF (filename == 'SUEWS_NonVeg.txt') THEN
605  IF (any(headernonveg_file /= headernonveg_reqd)) THEN
606  WRITE (*, *) headernonveg_file == headernonveg_reqd
607  WRITE (*, *) headernonveg_file
608  WRITE (*, *) headernonveg_reqd
609  CALL errorhint(56, 'Names or order of columns in SUEWS_NonVeg.txt does not match model code.', notused, notused, notusedi)
610  ENDIF
611 
612  ELSEIF (filename == 'SUEWS_Veg.txt') THEN
613  IF (any(headerveg_file /= headerveg_reqd)) THEN
614  WRITE (*, *) headerveg_file == headerveg_reqd
615  CALL errorhint(56, 'Names or order of columns in SUEWS_Veg.txt does not match model code.', notused, notused, notusedi)
616  ENDIF
617 
618  ELSEIF (filename == 'SUEWS_Water.txt') THEN
619  IF (any(headerwater_file /= headerwater_reqd)) THEN
620  WRITE (*, *) headerwater_file == headerwater_reqd
621  WRITE (*, *) headerwater_file
622  WRITE (*, *) headerwater_reqd
623  CALL errorhint(56, 'Names or order of columns in SUEWS_Water.txt does not match model code.', notused, notused, notusedi)
624  ENDIF
625 
626  ELSEIF (filename == 'SUEWS_Snow.txt') THEN
627  IF (any(headersnow_file /= headersnow_reqd)) THEN
628  WRITE (*, *) headersnow_file == headersnow_reqd
629  CALL errorhint(56, 'Names or order of columns in SUEWS_Snow.txt does not match model code.', notused, notused, notusedi)
630  ENDIF
631 
632  ELSEIF (filename == 'SUEWS_Soil.txt') THEN
633  IF (any(headersoil_file /= headersoil_reqd)) THEN
634  WRITE (*, *) headersoil_file == headersoil_reqd
635  CALL errorhint(56, 'Names or order of columns in SUEWS_Soil.txt does not match model code.', notused, notused, notusedi)
636  ENDIF
637 
638  ELSEIF (filename == 'SUEWS_Conductance.txt') THEN
639  IF (any(headercond_file /= headercond_reqd)) THEN
640  WRITE (*, *) headercond_file == headercond_reqd
641  CALL errorhint(56, 'Names or order of columns in SUEWS_Cond.txt does not match model code.', notused, notused, notusedi)
642  ENDIF
643 
644  ELSEIF (filename == 'SUEWS_OHMCoefficients.txt') THEN
647  CALL errorhint(56, 'Names or order of columns in SUEWS_OHMCoefficients.txt does not match model code.', &
649  ENDIF
650 
651  ELSEIF (filename == 'SUEWS_ESTMCoefficients.txt') THEN
654  CALL errorhint(56, 'Names or order of columns in SUEWS_ESTMCoefficients.txt does not match model code.', &
656  ENDIF
657 
658  ELSEIF (filename == 'SUEWS_AnthropogenicEmission.txt') THEN
661  CALL errorhint(56, 'Names or order of columns in SUEWS_AnthropogenicEmission.txt does not match model code.', &
663  ENDIF
664 
665  ELSEIF (filename == 'SUEWS_Irrigation.txt') THEN
668  WRITE (*, *) headerirrigation_file
669  WRITE (*, *) headerirrigation_reqd
670  CALL errorhint(56, 'Names or order of columns in SUEWS_Irrigation.txt does not match model code.', &
672  ENDIF
673 
674  ELSEIF (filename == 'SUEWS_Profiles.txt') THEN
675  IF (any(headerprofiles_file /= headerprofiles_reqd)) THEN
677  CALL errorhint(56, 'Names or order of columns in SUEWS_Profiles.txt does not match model code.', &
679  ENDIF
680 
681  ELSEIF (filename == 'SUEWS_WithinGridWaterDist.txt') THEN
684  CALL errorhint(56, 'Names or order of columns in SUEWS_WithinGridWaterDist.txt does not match model code.', &
686  ENDIF
687 
688  ELSEIF (filename == 'SUEWS_BiogenCO2.txt') THEN
689  IF (any(headerbiogen_file /= headerbiogen_reqd)) THEN
690  WRITE (*, *) headerbiogen_file == headerbiogen_reqd
691  CALL errorhint(56, 'Names or order of columns in SUEWS_BiogenCO2.txt does not match model code.', &
693  ENDIF
694 
695  ELSE
696  WRITE (*, *) 'Problem in subroutine InputHeaderCheck. File header not specified in model code for ', filename
698 
699  ENDIF
700 
701 ENDSUBROUTINE inputheadercheck
702 
703 !-------------------------------------------------------------------------
704 !
705 ! TS 05 Jul 2018: No longer needed as interpolation is done through specific subroutines at each required instant
706 !Interpolates hourly profiles provided in SUEWS_Profiles.txt
707 ! to resolution of the model timestep
708 ! HCW 06 Feb 2015
709 !===================================================================================
710 ! SUBROUTINE SUEWS_InterpHourlyProfiles(Gridiv,TstepP_ID,SurfChar_HrProf)
711 !
712 ! USE allocateArray
713 ! USE ColNamesInputFiles
714 ! USE sues_data
715 !
716 ! IMPLICIT NONE
717 !
718 ! INTEGER:: i,j, ii !Used to count over hours and sub-hourly timesteps
719 ! INTEGER:: Gridiv, TstepP_ID
720 ! INTEGER,DIMENSION(24):: SurfChar_HrProf
721 ! REAL(KIND(1d0)):: deltaProf !Change in hourly profiles per model timestep
722 !
723 ! ! Copy value for first hour
724 ! TstepProfiles(Gridiv,TstepP_ID,1) = SurfaceChar(Gridiv,SurfChar_HrProf(1))
725 ! DO i=1,24
726 ! j = (i+1)
727 ! IF(i == 24) j = 1 !If last hour of day, loop round to first hour of day for interpolation
728 ! deltaProf = ((SurfaceChar(Gridiv,SurfChar_HrProf(j)) - SurfaceChar(Gridiv,SurfChar_HrProf(i))))/nsh_real
729 ! DO ii=1,nsh
730 ! IF((nsh*(i-1)+ii+1) < (23*nsh+nsh+1)) THEN
731 ! TstepProfiles(Gridiv,TstepP_ID,(nsh*(i-1)+ii+1)) = SurfaceChar(Gridiv,SurfChar_HrProf(i)) + deltaProf*ii
732 ! ENDIF
733 ! ENDDO
734 ! ENDDO
735 !
736 ! endsubroutine SUEWS_InterpHourlyProfiles
737 !===================================================================================
738 
739 !===================================================================================
740 ! get interpolated profile values at specified time
741 ! NO normalisation performed
742 FUNCTION get_prof_spectime_inst(Hour, Min, Sec, Prof_24h) RESULT(Prof_CurrTime)
744  IMPLICIT NONE
745 
746  INTEGER :: i, j !Used to count over hours and sub-hourly timesteps
747  INTEGER, INTENT(IN) :: Hour, Min, Sec
748  INTEGER :: total_sec, SecPerHour
749  REAL(KIND(1d0)), DIMENSION(0:23), INTENT(IN) :: Prof_24h
750  REAL(KIND(1d0)):: deltaProf !Change in hourly profiles per model timestep
751  REAL(KIND(1d0)) :: Prof_CurrTime
752 
753  total_sec = min*60 + sec
754  secperhour = 3600
755 
756  i = hour
757  j = i + 1
758  IF (j == 24) j = 0
759 
760  deltaprof = (prof_24h(j) - prof_24h(i))/secperhour
761  prof_currtime = prof_24h(hour) + deltaprof*total_sec
762 
763 END FUNCTION get_prof_spectime_inst
764 
765 !===================================================================================
766 ! get interpolated profile values at specified time
767 ! normalise so the AVERAGE of the multipliers is equal to 1
768 FUNCTION get_prof_spectime_mean(Hour, Min, Sec, Prof_24h) RESULT(Prof_CurrTime)
770  IMPLICIT NONE
771 
772  INTEGER :: i, j !Used to count over hours and sub-hourly timesteps
773  INTEGER, INTENT(IN) :: Hour, Min, Sec
774  INTEGER :: total_sec, SecPerHour
775  REAL(KIND(1d0)), DIMENSION(0:23), INTENT(IN) :: Prof_24h
776  REAL(KIND(1d0)), DIMENSION(0:23):: Prof_24h_mean
777  REAL(KIND(1d0)):: deltaProf !Change in hourly profiles per model timestep
778  REAL(KIND(1d0)) :: Prof_CurrTime
779 
780  total_sec = min*60 + sec
781  secperhour = 3600
782 
783  prof_24h_mean = merge(prof_24h/(sum(prof_24h)/size(prof_24h, dim=1)), 0.d0, sum(prof_24h) /= 0) ! prevent zero-division
784  ! print*, Prof_24h_mean
785 
786  i = hour
787  j = i + 1
788  IF (j == 24) j = 0
789 
790  deltaprof = (prof_24h_mean(j) - prof_24h_mean(i))/secperhour
791 
792  ! print*, deltaProf,total_sec
793  prof_currtime = prof_24h_mean(i) + deltaprof*total_sec
794 
795 END FUNCTION get_prof_spectime_mean
796 
797 !===================================================================================
798 ! get interpolated profile values at specified time
799 ! normalise so the SUM of the multipliers is equal to 1
800 FUNCTION get_prof_spectime_sum(Hour, Min, Sec, Prof_24h, dt) RESULT(Prof_CurrTime)
802  IMPLICIT NONE
803 
804  INTEGER :: i, j !Used to count over hours and sub-hourly timesteps
805  INTEGER, INTENT(IN) :: Hour, Min, Sec, dt
806  INTEGER :: total_sec, SecPerHour
807  REAL(KIND(1d0)), DIMENSION(0:23), INTENT(IN) :: Prof_24h
808  REAL(KIND(1d0)), DIMENSION(0:23):: Prof_24h_sum
809  REAL(KIND(1d0)):: deltaProf !Change in hourly profiles per model timestep
810  REAL(KIND(1d0)) :: Prof_CurrTime
811 
812  total_sec = min*60 + sec
813  secperhour = 3600
814 
815  prof_24h_sum = merge(prof_24h/(sum(prof_24h)), 0.d0, sum(prof_24h) /= 0) ! prevent zero-division
816 
817  i = hour
818  j = i + 1
819  IF (j == 24) j = 0
820 
821  deltaprof = (prof_24h_sum(j) - prof_24h_sum(i))/secperhour
822  prof_currtime = prof_24h_sum(hour) + deltaprof*total_sec
823  prof_currtime = prof_currtime*dt/secperhour
824 
825 END FUNCTION get_prof_spectime_sum
826 !===================================================================================
827 
828 ! Subroutines for matching codes in the input files
829 ! could re-write as a generic function later...
830 
831 SUBROUTINE codematchohm(Gridiv, is, SWWD)
832  ! Matches OHM coefficients a1, a2, a3 via OHM codes
833  ! for summer/winter wet/dry conditions
834  ! HCW 03 Nov 2014
835  ! ---------------------------------------------------------
836 
837  USE allocatearray
838  USE initial
840  USE defaultnotused
841 
842  IMPLICIT NONE
843 
844  INTEGER:: gridiv
845  INTEGER:: is
846  CHARACTER(len=4):: SWWD
847 
848  iv5 = 0 ! Reset iv5 to zero
849 
850  IF (swwd == 'SWet') THEN
851 
852  DO iv5 = 1, nlinesohmcoefficients
853  IF (ohmcoefficients_coeff(iv5, co_code) == surfacechar(gridiv, c_ohmcode_swet(is))) THEN
854  EXIT
855  ELSEIF (iv5 == nlinesohmcoefficients) THEN
856  WRITE (*, *) 'Program stopped! OHM code (summer wet)', surfacechar(gridiv, c_ohmcode_swet(is)), &
857  'not found in OHM_Coefficients.txt for surface', is, '.'
858  CALL errorhint(57, 'Cannot find OHM code (summer wet)', surfacechar(gridiv, c_ohmcode_swet(is)), notused, notusedi)
859  ENDIF
860  ENDDO
861 
862  ELSEIF (swwd == 'SDry') THEN
863 
864  DO iv5 = 1, nlinesohmcoefficients
865  IF (ohmcoefficients_coeff(iv5, co_code) == surfacechar(gridiv, c_ohmcode_sdry(is))) THEN
866  EXIT
867  ELSEIF (iv5 == nlinesohmcoefficients) THEN
868  WRITE (*, *) 'Program stopped! OHM code (summer dry)', surfacechar(gridiv, c_ohmcode_sdry(is)), &
869  'not found in OHM_Coefficients.txt for surface', is, '.'
870  CALL errorhint(57, 'Cannot find OHM code (summer dry)', surfacechar(gridiv, c_ohmcode_sdry(is)), notused, notusedi)
871  ENDIF
872  ENDDO
873 
874  ELSEIF (swwd == 'WWet') THEN
875 
876  DO iv5 = 1, nlinesohmcoefficients
877  IF (ohmcoefficients_coeff(iv5, co_code) == surfacechar(gridiv, c_ohmcode_wwet(is))) THEN
878  EXIT
879  ELSEIF (iv5 == nlinesohmcoefficients) THEN
880  WRITE (*, *) 'Program stopped! OHM code (winter wet)', surfacechar(gridiv, c_ohmcode_wwet(is)), &
881  'not found in OHM_Coefficients.txt for surface', is, '.'
882  CALL errorhint(57, 'Cannot find OHM code (winter wet)', surfacechar(gridiv, c_ohmcode_wwet(is)), notused, notusedi)
883  ENDIF
884  ENDDO
885 
886  ELSEIF (swwd == 'WDry') THEN
887 
888  DO iv5 = 1, nlinesohmcoefficients
889  IF (ohmcoefficients_coeff(iv5, co_code) == surfacechar(gridiv, c_ohmcode_wdry(is))) THEN
890  EXIT
891  ELSEIF (iv5 == nlinesohmcoefficients) THEN
892  WRITE (*, *) 'Program stopped! OHM code (winter dry)', surfacechar(gridiv, c_ohmcode_wdry(is)), &
893  'not found in OHM_Coefficients.txt for surface', is, '.'
894  CALL errorhint(57, 'Cannot find OHM code (winter dry)', surfacechar(gridiv, c_ohmcode_wdry(is)), notused, notusedi)
895  ENDIF
896  ENDDO
897 
898  ELSE
899  WRITE (*, *) 'Problem with CodeMatchOHM (in SUEWS_CodeMatch.f95). ', swwd, ' not recognised. Needs to be one of: ', &
900  'SWet = Summer Wet, SDry = Summer Dry, WWet = WinterWet, WDry = Winter Dry. N.B. Case sensitive. '
901  stop
902  ENDIF
903 
904  RETURN
905 ENDSUBROUTINE codematchohm
906 ! ---------------------------------------------------------
907 
908 SUBROUTINE codematchestm(Gridiv, is)
909  ! Matches ESTM coefficients via ESTM code
910  ! Modified HCW 16 Jun 2016 - for SUEWS surface types
911  ! - removed summer/winter wet/dry option
912  ! S.O. 04 Feb 2016
913  ! ---------------------------------------------------------
914 
915  USE allocatearray
916  USE initial
918  USE defaultnotused
919 
920  IMPLICIT NONE
921 
922  INTEGER:: gridiv
923  INTEGER:: is
924 
925  iv5 = 0 ! Reset iv5 to zero
926 
928  IF (estmcoefficients_coeff(iv5, ce_code) == surfacechar(gridiv, c_estmcode(is))) THEN
929  EXIT
930  ELSEIF (iv5 == nlinesestmcoefficients) THEN
931  WRITE (*, *) 'Program stopped! ESTM code', surfacechar(gridiv, c_estmcode(is)), &
932  'not found in ESTM_Coefficients.txt for surface', is, '.'
933  CALL errorhint(57, 'Cannot find ESTM code', surfacechar(gridiv, c_estmcode(is)), notused, notusedi)
934  ENDIF
935  ENDDO
936 
937  RETURN
938 ENDSUBROUTINE codematchestm
939 ! ---------------------------------------------------------
940 
941 SUBROUTINE codematchestm_class(Gridiv, is, ii)
942  ! Matches ESTM coefficients via ESTM codes in SiteSelect for Paved and Bldgs ESTM classes
943  ! HCW 16 Jun 2016
944  ! ---------------------------------------------------------
945 
946  USE allocatearray
947  USE initial
949  USE defaultnotused
950 
951  IMPLICIT NONE
952 
953  INTEGER:: gridiv
954  INTEGER:: is, ii
955 
956  iv5 = 0 ! Reset iv5 to zero
957 
958  IF (is == bldgsurf) THEN
961  EXIT
962  ELSEIF (iv5 == nlinesestmcoefficients) THEN
963  WRITE (*, *) 'Program stopped! ESTM code', surfacechar(gridiv, c_code_estmclass_bldgs(ii)), &
964  'not found in ESTM_Coefficients.txt for surface', is, '.'
965  CALL errorhint(57, 'Cannot find ESTM code', surfacechar(gridiv, c_code_estmclass_bldgs(ii)), notused, notusedi)
966  ENDIF
967  ENDDO
968  ELSEIF (is == pavsurf) THEN
971  EXIT
972  ELSEIF (iv5 == nlinesestmcoefficients) THEN
973  WRITE (*, *) 'Program stopped! ESTM code', surfacechar(gridiv, c_code_estmclass_paved(ii)), &
974  'not found in ESTM_Coefficients.txt for surface', is, '.'
975  CALL errorhint(57, 'Cannot find ESTM code', surfacechar(gridiv, c_code_estmclass_paved(ii)), notused, notusedi)
976  ENDIF
977  ENDDO
978  ELSE
979  WRITE (*, *) 'Problem with CodeMatchESTM_Class (in SUEWS_CodeMatch.f95). ', is, ' not correct. Needs to be either ', &
980  '1 = Paved surfaced, 2 = Bldgs surfaces.'
981  stop
982  ENDIF
983  RETURN
984 ENDSUBROUTINE codematchestm_class
985 ! ---------------------------------------------------------
986 
987 SUBROUTINE codematchprof(Gridiv, SurfaceCharCodeCol)
988  ! Matches Soil characteristics via codes in *SurfaceChar*
989  ! for energy use/water use/snow clearing
990  ! HCW 20 Nov 2014
991  ! ---------------------------------------------------------
992 
993  USE allocatearray
994  USE initial
996  USE defaultnotused
997 
998  IMPLICIT NONE
999 
1000  INTEGER:: Gridiv
1001  INTEGER:: SurfaceCharCodeCol
1002 
1003  iv5 = 0 ! Reset iv5 to zero
1004 
1005  DO iv5 = 1, nlinesprofiles
1006  IF (profiles_coeff(iv5, cpr_code) == surfacechar(gridiv, surfacecharcodecol)) THEN
1007  EXIT
1008  ELSEIF (iv5 == nlinesprofiles) THEN
1009  WRITE (*, *) 'Program stopped! Profile code ', surfacechar(gridiv, surfacecharcodecol), 'not found in SUEWS_Profiles.txt.'
1010  CALL errorhint(57, 'Cannot find code in SUEWS_Profiles.txt', surfacechar(gridiv, surfacecharcodecol), notused, notusedi)
1011  ENDIF
1012  ENDDO
1013 
1014  RETURN
1015 ENDSUBROUTINE codematchprof
1016 ! ---------------------------------------------------------
1017 
1018 SUBROUTINE codematchdist(rr, CodeCol, codeColSameSurf)
1019  ! Matches within-grid water distribution via codes
1020  ! Checks water cannot flow from one surface to the same surface
1021  ! Checks water distribution fractions sum to 1.
1022  ! HCW 10 Nov 2014
1023  ! ---------------------------------------------------------
1024 
1025  USE allocatearray
1026  USE initial
1027  USE colnamesinputfiles
1028  USE defaultnotused
1029 
1030  IMPLICIT NONE
1031 
1032  INTEGER:: rr
1033  INTEGER:: codeCol, codeColSameSurf
1034 
1035  iv5 = 0 ! Reset iv5 to zero
1036 
1037  DO iv5 = 1, nlineswgwaterdist
1038  IF (wgwaterdist_coeff(iv5, cwg_code) == siteselect(rr, codecol)) THEN
1039  EXIT
1040  ELSEIF (iv5 == nlineswgwaterdist) THEN
1041  WRITE (*, *) 'Program stopped! Within-grid water distribution code ', siteselect(rr, codecol), &
1042  'not found in SUEWS_WaterDistWithinGrid.txt.'
1043  CALL errorhint(57, 'Cannot find code in SUEWS_WaterDistWithinGrid.txt', siteselect(rr, codecol), notused, notusedi)
1044  ENDIF
1045  ENDDO
1046 
1047  ! Check water flow to same surface is zero (previously in RunControlByGridByYear in SUEWS_Initial.f95)
1048  IF (wgwaterdist_coeff(iv5, codecolsamesurf) /= 0) THEN
1049  CALL errorhint(8, 'Diagonal elements should be zero as water cannot move from one surface to the same surface.', &
1050  wgwaterdist_coeff(iv5, codecolsamesurf), notused, notusedi)
1051  ENDIF
1052 
1053  !! QUESTION: MODIFY THIS?
1054  ! Check water either moves to runoff or soilstore, but not to both
1055  ! Model returns an error if both ToRunoff and ToSoilStore are non-zero.
1056  !! - Probably should remove this...
1057  !! - Also look at SUEWS_translate, as the non-zero value goes into WaterDist
1059  CALL errorhint(9, 'One of these (ToRunoff,ToSoilStore) should be zero.', &
1061  ENDIF
1062 
1063  !! Also do for water surface once implemented
1064  IF (codecol /= c_wgwatercode) THEN ! Except for Water surface
1065  ! Check total water distribution from each surface adds up to 1
1066  IF (sum(wgwaterdist_coeff(iv5, cwg_topaved:cwg_tosoilstore)) > 1.0000001 &
1067  .OR. sum(wgwaterdist_coeff(iv5, cwg_topaved:cwg_tosoilstore)) < 0.9999999) THEN
1068  CALL errorhint(8, 'Total water distribution from each surface should add up to 1.', &
1070  ENDIF
1071  ENDIF
1072 
1073  RETURN
1074 ENDSUBROUTINE codematchdist
1075 ! ---------------------------------------------------------
1076 
1077 SUBROUTINE codematchnonveg(rr, CodeCol)
1078  ! Matches Impervious characteristics via codes in SiteSelect
1079  ! HCW 20 Nov 2014
1080  ! ---------------------------------------------------------
1081 
1082  USE allocatearray
1083  USE initial
1084  USE colnamesinputfiles
1085  USE defaultnotused
1086 
1087  IMPLICIT NONE
1088 
1089  INTEGER:: rr
1090  INTEGER:: codeCol
1091 
1092  iv5 = 0 ! Reset iv5 to zero
1093 
1094  DO iv5 = 1, nlinesnonveg
1095  IF (nonveg_coeff(iv5, ci_code) == siteselect(rr, codecol)) THEN
1096  EXIT
1097  ELSEIF (iv5 == nlinesnonveg) THEN
1098  WRITE (*, *) 'Program stopped! NonVeg code ', siteselect(rr, codecol), 'not found in SUEWS_NonVeg.txt.'
1099  CALL errorhint(57, 'Cannot find code in SUEWS_NonVeg.txt', siteselect(rr, codecol), notused, notusedi)
1100  ENDIF
1101  ENDDO
1102 
1103  RETURN
1104 ENDSUBROUTINE codematchnonveg
1105 ! ---------------------------------------------------------
1106 
1107 SUBROUTINE codematchveg(rr, CodeCol)
1108  ! Matches Pervious characteristics via codes in SiteSelect
1109  ! HCW 20 Nov 2014
1110  ! ---------------------------------------------------------
1111 
1112  USE allocatearray
1113  USE initial
1114  USE colnamesinputfiles
1115  USE defaultnotused
1116 
1117  IMPLICIT NONE
1118 
1119  INTEGER:: rr
1120  INTEGER:: codeCol
1121 
1122  iv5 = 0 ! Reset iv5 to zero
1123 
1124  DO iv5 = 1, nlinesveg
1125  IF (veg_coeff(iv5, cp_code) == siteselect(rr, codecol)) THEN
1126  EXIT
1127  ELSEIF (iv5 == nlinesveg) THEN
1128  WRITE (*, *) 'Program stopped! Veg code ', siteselect(rr, codecol), 'not found in SUEWS_Vegs.txt.'
1129  CALL errorhint(57, 'Cannot find code in SUEWS_Veg.txt', siteselect(rr, codecol), notused, notusedi)
1130  ENDIF
1131  ENDDO
1132 
1133  RETURN
1134 ENDSUBROUTINE codematchveg
1135 ! ---------------------------------------------------------
1136 
1137 SUBROUTINE codematchwater(rr, CodeCol)
1138  ! Matches Water characteristics via codes in SiteSelect
1139  ! HCW 20 Nov 2014
1140  ! ---------------------------------------------------------
1141 
1142  USE allocatearray
1143  USE initial
1144  USE colnamesinputfiles
1145  USE defaultnotused
1146 
1147  IMPLICIT NONE
1148 
1149  INTEGER:: rr
1150  INTEGER:: codeCol
1151 
1152  iv5 = 0 ! Reset iv5 to zero
1153 
1154  DO iv5 = 1, nlineswater
1155  IF (water_coeff(iv5, cw_code) == siteselect(rr, codecol)) THEN
1156  EXIT
1157  ELSEIF (iv5 == nlineswater) THEN
1158  WRITE (*, *) 'Program stopped! Water code ', siteselect(rr, codecol), 'not found in SUEWS_Water.txt.'
1159  CALL errorhint(57, 'Cannot find code in SUEWS_Water.txt', siteselect(rr, codecol), notused, notusedi)
1160  ENDIF
1161  ENDDO
1162 
1163  RETURN
1164 ENDSUBROUTINE codematchwater
1165 ! ---------------------------------------------------------
1166 
1167 SUBROUTINE codematchsnow(rr, CodeCol)
1168  ! Matches Snow characteristics via codes in SiteSelect
1169  ! HCW 20 Nov 2014
1170  ! ---------------------------------------------------------
1171 
1172  USE allocatearray
1173  USE initial
1174  USE colnamesinputfiles
1175  USE defaultnotused
1176 
1177  IMPLICIT NONE
1178 
1179  INTEGER:: rr
1180  INTEGER:: codeCol
1181 
1182  iv5 = 0 ! Reset iv5 to zero
1183 
1184  DO iv5 = 1, nlinessnow
1185  IF (snow_coeff(iv5, cs_code) == siteselect(rr, codecol)) THEN
1186  EXIT
1187  ELSEIF (iv5 == nlinessnow) THEN
1188  WRITE (*, *) 'Program stopped! Snow code ', siteselect(rr, codecol), 'not found in SUEWS_Snow.txt.'
1189  CALL errorhint(57, 'Cannot find code in SUEWS_Snow.txt', siteselect(rr, codecol), notused, notusedi)
1190  ENDIF
1191  ENDDO
1192 
1193  RETURN
1194 ENDSUBROUTINE codematchsnow
1195 ! ---------------------------------------------------------
1196 
1197 SUBROUTINE codematchconductance(rr, CodeCol)
1198  ! Matches Conductance characteristics via codes in SiteSelect
1199  ! HCW 20 Nov 2014
1200  ! ---------------------------------------------------------
1201 
1202  USE allocatearray
1203  USE initial
1204  USE colnamesinputfiles
1205  USE defaultnotused
1206 
1207  IMPLICIT NONE
1208 
1209  INTEGER:: rr
1210  INTEGER:: codeCol
1211 
1212  iv5 = 0 ! Reset iv5 to zero
1213 
1214  DO iv5 = 1, nlinesconductance
1215  IF (conductance_coeff(iv5, cc_code) == siteselect(rr, codecol)) THEN
1216  EXIT
1217  ELSEIF (iv5 == nlinesconductance) THEN
1218  WRITE (*, *) 'Program stopped! Conductance code ', siteselect(rr, codecol), 'not found in SUEWS_Conductance.txt.'
1219  CALL errorhint(57, 'Cannot find code in SUEWS_Conductance.txt', siteselect(rr, codecol), notused, notusedi)
1220  ENDIF
1221  ENDDO
1222 
1223  RETURN
1224 ENDSUBROUTINE codematchconductance
1225 ! ---------------------------------------------------------
1226 
1227 SUBROUTINE codematchanthropogenic(rr, CodeCol)
1228  ! Matches AnthropogenicHeat characteristics via codes in SiteSelect
1229  ! HCW 20 Nov 2014
1230  ! MH 21 Jun 2017
1231  ! ---------------------------------------------------------
1232 
1233  USE allocatearray
1234  USE initial
1235  USE colnamesinputfiles
1236  USE defaultnotused
1237 
1238  IMPLICIT NONE
1239 
1240  INTEGER:: rr
1241  INTEGER:: codeCol
1242 
1243  iv5 = 0 ! Reset iv5 to zero
1244 
1245  DO iv5 = 1, nlinesanthropogenic
1246  IF (anthropogenic_coeff(iv5, ca_code) == siteselect(rr, codecol)) THEN
1247  EXIT
1248  ELSEIF (iv5 == nlinesanthropogenic) THEN
1249  WRITE (*, *) 'Program stopped! Anthropogenic code ', siteselect(rr, codecol), &
1250  'not found in SUEWS_AnthropogenicEmission.txt.'
1251  CALL errorhint(57, 'Cannot find code in SUEWS_AnthropogenicEmission.txt', &
1252  siteselect(rr, codecol), notused, notusedi)
1253  ENDIF
1254  ENDDO
1255 
1256  RETURN
1257 ENDSUBROUTINE codematchanthropogenic
1258 ! ---------------------------------------------------------
1259 
1260 SUBROUTINE codematchirrigation(rr, CodeCol)
1261  ! Matches Irrigation characteristics via codes in SiteSelect
1262  ! HCW 20 Nov 2014
1263  ! ---------------------------------------------------------
1264 
1265  USE allocatearray
1266  USE initial
1267  USE colnamesinputfiles
1268  USE defaultnotused
1269 
1270  IMPLICIT NONE
1271 
1272  INTEGER:: rr
1273  INTEGER:: codeCol
1274 
1275  iv5 = 0 ! Reset iv5 to zero
1276 
1277  DO iv5 = 1, nlinesirrigation
1278  IF (irrigation_coeff(iv5, cir_code) == siteselect(rr, codecol)) THEN
1279  EXIT
1280  ELSEIF (iv5 == nlinesirrigation) THEN
1281  WRITE (*, *) 'Program stopped! Irrigation code ', siteselect(rr, codecol), 'not found in SUEWS_Irrigation.txt.'
1282  CALL errorhint(57, 'Cannot find code in SUEWS_Irrigation.txt', siteselect(rr, codecol), notused, notusedi)
1283  ENDIF
1284  ENDDO
1285 
1286  RETURN
1287 ENDSUBROUTINE codematchirrigation
1288 ! ---------------------------------------------------------
1289 
1290 SUBROUTINE codematchsoil(Gridiv, SurfaceCharCodeCol)
1291  ! Matches Soil characteristics via codes in *SurfaceChar*
1292  ! HCW 20 Nov 2014
1293  ! ---------------------------------------------------------
1294 
1295  USE allocatearray
1296  USE initial
1297  USE colnamesinputfiles
1298  USE defaultnotused
1299 
1300  IMPLICIT NONE
1301 
1302  INTEGER:: Gridiv
1303  INTEGER:: SurfaceCharCodeCol
1304 
1305  iv5 = 0 ! Reset iv5 to zero
1306 
1307  DO iv5 = 1, nlinessoil
1308  IF (soil_coeff(iv5, cso_code) == surfacechar(gridiv, surfacecharcodecol)) THEN
1309  EXIT
1310  ELSEIF (iv5 == nlinessoil) THEN
1311  WRITE (*, *) 'Program stopped! Soil code ', surfacechar(gridiv, surfacecharcodecol), 'not found in SUEWS_Soil.txt.'
1312  CALL errorhint(57, 'Cannot find code in SUEWS_Soil.txt', surfacechar(gridiv, surfacecharcodecol), notused, notusedi)
1313  ENDIF
1314  ENDDO
1315 
1316  RETURN
1317 ENDSUBROUTINE codematchsoil
1318 ! ---------------------------------------------------------
1319 
1320 SUBROUTINE codematchbiogen(Gridiv, SurfaceCharCodeCol)
1321  ! Matches Biogen characteristics via codes in *SuraceChar*
1322  ! MH 16 Jun 2017
1323  ! ---------------------------------------------------------
1324 
1325  USE allocatearray
1326  USE initial
1327  USE colnamesinputfiles
1328  USE defaultnotused
1329 
1330  IMPLICIT NONE
1331 
1332  INTEGER:: Gridiv
1333  INTEGER:: SurfaceCharCodeCol
1334 
1335  iv5 = 0 ! Reset iv5 to zero
1336 
1337  DO iv5 = 1, nlinesbiogen
1338  IF (biogen_coeff(iv5, cb_code) == surfacechar(gridiv, surfacecharcodecol)) THEN
1339  EXIT
1340  ELSEIF (iv5 == nlinesbiogen) THEN
1341  WRITE (*, *) 'Program stopped! Biogen code ', surfacechar(gridiv, surfacecharcodecol), 'not found in SUEWS_BiogenCO2.txt.'
1342  CALL errorhint(57, 'Cannot find code in SUEWS_BiogenCO2.txt', surfacechar(gridiv, surfacecharcodecol), notused, notusedi)
1343  ENDIF
1344  ENDDO
1345 
1346  RETURN
1347 ENDSUBROUTINE codematchbiogen
1348 ! ---------------------------------------------------------
1349 
1351  !========================================================================================
1352  ! Disaggregation of meteorological forcing data
1353  ! Code to disaggregate met forcing data from resolution provided to the model time-step
1354  ! Created: HCW 10 Feb 2017
1355  !
1356  ! ---- Key for MetDisaggMethod settings ----
1357  !10 -> linear disaggregation for timestamps representing the end of the averaging period
1358  !20 -> linear disaggregation for instantaneous variables (e.g. air temperature, humidity, pressure, wind speed in WFDEI dataset)
1359  !100 -> evenly distribute rainfall among all subintervals in a rainy interval
1360  !101 -> evenly distribute rainfall among RainAmongN subintervals in a rainy interval
1361  ! - requires RainAmongN to be set in RunControl.nml
1362  ! If KdownZen = 1 -> include additional zenith check in kdown disaggregation
1363  !
1364  ! N.B. wdir downscaling is currently not implemented
1365  !
1366  !========================================================================================
1367 
1368  USE allocatearray
1369  USE colnamesinputfiles
1370  USE data_in
1371  USE initial
1373 
1374  IMPLICIT NONE
1375 
1376 CONTAINS
1377 
1378  !======================================================================================
1379  SUBROUTINE disaggregatemet(iBlock, igrid)
1380  ! Subroutine to disaggregate met forcing data to model time-step
1381  ! HCW 10 Feb 2017
1382  !======================================================================================
1383 
1384  USE sues_data
1385  USE defaultnotused
1386 
1387  IMPLICIT NONE
1388 
1389  INTEGER:: lunit = 100
1390  INTEGER:: tdiff !Time difference (in minutes) between first and second rows of original met forcing file
1391  INTEGER:: i, ii !counter
1392  INTEGER:: iBlock, igrid
1393  INTEGER, DIMENSION(Nper):: seq1Nper
1394  INTEGER, DIMENSION(nsd):: seq1nsd
1395  INTEGER, DIMENSION(nColumnsMetForcingData):: MetDisaggMethod ! Stores method to use for disaggregating met data
1396  REAL(KIND(1d0)), DIMENSION(nColumnsMetForcingData):: MetArrayOrig
1397  REAL(KIND(1d0)), DIMENSION(ReadLinesOrigMetData*Nper, ncolumnsMetForcingData):: Met_tt
1398  REAL(KIND(1d0)), DIMENSION(ReadLinesOrigMetData*Nper):: Met_tt_kdownAdj
1399  CHARACTER(LEN=9), DIMENSION(ncolumnsMetForcingData):: HeaderMet
1400  CHARACTER(LEN=10*ncolumnsMetForcingData):: HeaderMetOut
1401  ! REAL(KIND(1d0)),DIMENSION(ReadLinesOrigMetData):: dectimeOrig
1402  ! REAL(KIND(1d0)),DIMENSION(ReadLinesOrigMetData*Nper):: dectimeDscd, dectimeFast
1403  REAL(KIND(1d0)), DIMENSION(ReadLinesOrigMetData*Nper):: dectimeFast
1404  REAL(KIND(1d0)), DIMENSION(ReadLinesOrigMetData*Nper):: idectime ! sun position at middle of time-step before
1405 
1406  ! INTEGER, DIMENSION(Nper):: temp_iy, temp_id, temp_ih, temp_im, temp_ihm ! temorary varaibles for disaggragation
1407  ! REAL(KIND(1d0)), DIMENSION(Nper):: temp_dectime ! temorary varaibles for disaggragation
1408  ! INTEGER :: Days_of_Year
1409  ! INTEGER, DIMENSION(Nper)::ndays_iy ! number of days in iy
1410 
1411  ! Allocate and initialise arrays to receive original forcing data --------------------
1415  metfordisagg(:, :) = -999
1416  metfordisaggprev(:) = -999
1417  metfordisaggnext(:) = -999
1418  ! Initialise array to receive disaggregated data
1419  met_tt = -999
1420  ! Intialise disaggregation method
1421  metdisaggmethod = -999
1422 
1423  ! Generate useful sequences
1424  seq1nper = (/(i, i=1, nper, 1)/)
1425  seq1nsd = (/(i, i=1, nsd, 1)/)
1426 
1427  ! Get methods to use for disaggregation from RunControl
1428  IF (diagnosedisagg == 1) WRITE (*, *) 'DisaggMethod: ', disaggmethod, 'RainDisaggMethod:', raindisaggmethod
1429  IF (disaggmethod == 1) THEN
1430  metdisaggmethod(:) = 10 !linear disaggregation of averages
1431  ELSEIF (disaggmethod == 2) THEN
1432  metdisaggmethod(:) = 20 !linear disaggregation of instantaneous values
1433  ELSEIF (disaggmethod == 3) THEN !WFDEI set up, where T, Q, pres, U are instantaneous
1434  metdisaggmethod(:) = 10 !linear disaggregation of averages
1435  metdisaggmethod(10:13) = 20 !linear disagg instantaneous values for U, RH, Tair, pres
1436  ELSE
1437  CALL errorhint(2, 'Problem in SUEWS_MetDisagg: DisaggMethod value should be 1, 2, or 3', &
1439  ENDIF
1440  ! Set rainfall
1441  metdisaggmethod(14) = raindisaggmethod
1442 
1443  ! Read data ---------------------------------------------------------------------
1444  IF (diagnosedisagg == 1) WRITE (*, *) 'Reading file: ', trim(fileorigmet)
1445  OPEN (lunit, file=trim(fileorigmet), status='old')
1446  ! CALL skipHeader(lunit,SkipHeaderMet) !Skip header -> read header instead
1447  READ (lunit, *) headermet
1448  !write(*,*) HeaderMet
1449  ! Skip over lines that have already been read and downscaled
1450  IF (skippedlinesorig > 0) THEN
1451  DO i = 1, skippedlinesorig - 1 ! minus 1 here because last line of last block needs to be read again
1452  READ (lunit, *)
1453  ENDDO
1454  ! Read in last line of previous block
1455  CALL metread(lunit, metarrayorig, inputmetformat, ldown_option, netradiationmethod, &
1457  metfordisaggprev(1:ncolumnsmetforcingdata) = metarrayorig
1458  ENDIF
1459  ! print*, 'MetForDisagg',MetForDisagg(1:3,1:4)
1460  ! print*, 'ReadLinesOrigMetDataMax',ReadLinesOrigMetDataMax
1461  ! Read in current block
1462  DO i = 1, readlinesorigmetdatamax
1463  CALL metread(lunit, metarrayorig, inputmetformat, ldown_option, netradiationmethod, &
1465  metfordisagg(i, 1:ncolumnsmetforcingdata) = metarrayorig
1466  ENDDO
1467  ! print*, 'MetForDisagg',MetForDisagg(1:3,1:4)
1468  ! Read in first line of next block (except for last block)
1469  IF (iblock /= readblocksorigmetdata) THEN
1470  CALL metread(lunit, metarrayorig, inputmetformat, ldown_option, netradiationmethod, &
1472  metfordisaggnext(1:ncolumnsmetforcingdata) = metarrayorig
1473  ENDIF
1474  CLOSE (lunit)
1475 
1476  ! Check resolution of original met forcing data -------------------------------------
1477  ! Find time difference (in minutes) between first and second row
1478  tdiff = int(metfordisagg(2, 4) - metfordisagg(1, 4)) !Try using minutes
1479  IF (tdiff == 0) tdiff = int(metfordisagg(2, 3) - metfordisagg(1, 3))*60 !If no difference in minutes, try using hours
1480  IF (tdiff < 0) THEN !If time difference is negative (e.g. change of day), instead use second and third row
1481  tdiff = int(metfordisagg(3, 4) - metfordisagg(2, 4))
1482  IF (tdiff == 0) tdiff = int(metfordisagg(3, 3) - metfordisagg(2, 3))*60 !If no difference in minutes, try using hours
1483  ENDIF
1484  ! Check actual resolution matches specified input resolution
1485  IF (tdiff /= resolutionfilesin/60) THEN
1486  CALL errorhint(2, 'Problem in SUEWS_MetDisagg: timestamps in met forcing file inconsistent with ResolutionFilesIn', &
1487  REAL(ResolutionFilesIn, KIND(1d0)), NotUsed, tdiff*60)
1488  ENDIF
1489 
1490  ! Check file only contains a single year --------------------------------------------
1491  ! Very last data point is allowed to be (should be) timestamped with following year
1492  IF (any(metfordisagg(1:(readlinesorigmetdatamax - 1), 1) /= metfordisagg(1, 1))) THEN
1493  CALL errorhint(3, 'Problem in SUEWS_MetDisagg: multiple years found in original met forcing file.', &
1494  metfordisagg(1, 1), notused, notusedi)
1495  ENDIF
1496 
1497  ! Disaggregate time columns ---------------------------------------------------------
1498  IF (diagnose == 1) WRITE (*, *) 'Disaggregating met forcing data (', trim(fileorigmet), ') to model time-step...'
1499 
1500  CALL disaggregatedatetime(metfordisagg(:, 1:4), tstep, nper, readlinesorigmetdatamax, met_tt(:, 1:4))
1501 
1502  ! Disaggregate other columns --------------------------------------------------------
1503  DO ii = 5, ncolumnsmetforcingdata
1504  IF (ii == 14) THEN !Do something different for rainfall and snowfall (if present)
1505  IF (metdisaggmethod(14) == 100) THEN
1507  IF (all(metfordisagg(:, 16) == -999)) THEN
1508  met_tt(:, 16) = -999
1509  ELSE
1511  ENDIF
1512  ELSEIF (metdisaggmethod(14) == 101) THEN
1513  IF (rainamongn == -999) THEN
1514  CALL errorhint(2, 'Problem in SUEWS_MetDisagg: RainDisaggMethod requires RainAmongN', &
1515  REAL(RainAmongN, KIND(1d0)), NotUsed, RainDisaggMethod)
1516  ELSEIF (rainamongn > nper) THEN
1517  CALL errorhint(2, 'Problem in SUEWS_MetDisagg: RainAmongN > Nper', REAL(Nper, KIND(1d0)), NotUsed, RainAmongN)
1518  ELSE
1519  met_tt(:, 14) = disaggp_amongn(metfordisagg(:, 14), &
1521  IF (all(metfordisagg(:, 16) == -999)) THEN
1522  met_tt(:, 16) = -999
1523  ELSE
1524  met_tt(:, 16) = disaggp_amongn(metfordisagg(:, 16), &
1526  ENDIF
1527  ENDIF
1528  ELSEIF (metdisaggmethod(14) == 102) THEN
1529  IF (all(multrainamongn == -999)) THEN
1530  CALL errorhint(2, 'Problem in SUEWS_MetDisagg: RainDisaggMethod requires MultRainAmongN', &
1531  REAL(MultRainAmongN(1), KIND(1d0)), NotUsed, RainDisaggMethod)
1532  ELSEIF (all(multrainamongnupperi == -999)) THEN
1533  CALL errorhint(2, 'Problem in SUEWS_MetDisagg: RainDisaggMethod requires MultRainAmongNUpperI', &
1534  multrainamongnupperi(1), notused, raindisaggmethod)
1535  ELSEIF (any(multrainamongn > nper)) THEN
1536  CALL errorhint(2, 'Problem in SUEWS_MetDisagg: MultRainAmongN > Nper', REAL(Nper, KIND(1d0)), NotUsed, &
1537  MAXVAL(multrainamongn))
1538  ELSE
1541  IF (all(metfordisagg(:, 16) == -999)) THEN
1542  met_tt(:, 16) = -999
1543  ELSE
1546  ENDIF
1547  ENDIF
1548  ELSE
1549  WRITE (*, *) 'Disaggregation code for rain not recognised'
1550  ENDIF
1551  ELSEIF (ii == 24) THEN !wind direction disaggregation not coded yet...
1552  IF (any(metfordisagg(:, ii) /= -999)) THEN
1553  WRITE (*, *) 'Disaggregation of wind direction not currently implemented!'
1554  ENDIF
1555  ELSE
1556  IF (all(metfordisagg(:, ii) == -999)) THEN
1557  !IF(DiagnoseDisagg==1) write(*,*) 'No data for col.', ii
1558  met_tt(:, ii) = -999
1559  ELSE
1560  met_tt(:, ii) = disagg_lin(metfordisagg(:, ii), metfordisaggprev(ii), metfordisaggnext(ii), metdisaggmethod(ii), &
1562  ENDIF
1563  ENDIF
1564  ENDDO
1565 
1566  ! Adjust kdown disaggregation using zenith angle
1567  IF (kdownzen == 1) THEN
1568  IF (diagnosedisagg == 1) WRITE (*, *) 'Adjusting disaggregated kdown using zenith angle'
1569  met_tt_kdownadj(:) = met_tt(:, 15)
1570  ! Translate location data from SurfaceChar to find solar angles
1571  lat = surfacechar(igrid, c_lat)
1572  lng = surfacechar(igrid, c_lng)
1573  timezone = surfacechar(igrid, c_tz)
1574  alt = surfacechar(igrid, c_alt)
1575  ! Calculate dectime at downscaled time-step
1576  dectimefast(:) = met_tt(:, 2) + met_tt(:, 3)/24.0 + met_tt(:, 4)/(60.0*24.0)
1577  idectime = dectimefast - halftimestep! sun position at middle of timestep before
1578  DO i = 1, (readlinesorigmetdatamax*nper)
1579  CALL narp_cal_sunposition(met_tt(i, 2), idectime(i), timezone, lat, lng, alt, azimuth, zenith_deg)
1580  ! If sun below horizon, set disaggregated kdown to zero
1581  IF (zenith_deg > 90) THEN
1582  !write(*,*) Met_tt(i,1:4)
1583  met_tt_kdownadj(i) = 0.0
1584  ENDIF
1585  ENDDO
1586  ! Redistribute kdown over each day
1587  DO i = 1, (readlinesorigmetdatamax*nper/nsd) ! Loop over each day
1588  met_tt_kdownadj((i - 1)*nsd + seq1nsd) = &
1589  met_tt_kdownadj((i - 1)*nsd + seq1nsd) &
1590  *sum(met_tt((i - 1)*nsd + seq1nsd, 15)) &
1591  /sum(met_tt_kdownadj((i - 1)*nsd + seq1nsd))
1592  ENDDO
1593  ! Copy adjusted kdown back to Met_tt array
1594  met_tt(:, 15) = met_tt_kdownadj(:)
1595  ENDIF
1596 
1597  ! Copy disaggregated data to MetForcingDataArray
1598  metforcingdata(:, 1:24, gridcounter) = met_tt(:, 1:24)
1599 
1600  ! If snow is -999, set to zero (also in LUMPS_metRead.f95)
1601  IF (all(metforcingdata(:, 16, gridcounter) == -999)) metforcingdata(:, 16, gridcounter) = 0
1602 
1603  ! Undo pressure conversion again for writing out
1604  met_tt(:, 13) = met_tt(:, 13)/10.0
1605 
1606  ! Write out disaggregated file ------------------------------------------------------
1607  IF (keeptstepfilesin == 1) THEN
1608  IF (iblock == 1) THEN
1609  ! Prepare header
1610  DO i = 1, ncolumnsmetforcingdata
1611  IF (i == 1) THEN
1612  headermetout = adjustl(headermet(i))
1613  ELSE
1614  headermetout = trim(headermetout)//' '//adjustl(headermet(i))
1615  ENDIF
1616  ENDDO
1617  ! Write out header
1618  OPEN (78, file=trim(filedscdmet), err=112)
1619  WRITE (78, '(a)') headermetout
1620  ELSE
1621  OPEN (78, file=trim(filedscdmet), position='append')!,err=112)
1622  ENDIF
1623  ! Write out data
1624  DO i = 1, (readlinesorigmetdatamax*nper)
1625  WRITE (78, 303) (int(met_tt(i, ii)), ii=1, 4), met_tt(i, 5:ncolumnsmetforcingdata)
1626  ENDDO
1627  !IF(iBlock == ReadBlocksOrigMetData) THEN
1628  ! WRITE(78,'(i2)') -9
1629  ! WRITE(78,'(i2)') -9
1630  !ENDIF
1631  CLOSE (78) !Close output file
1632  ENDIF
1633 
1634 303 FORMAT((i4, 1x), 3(i3, 1x), 9(f12.6, 1x), (f9.4, 1x), 10(f9.4, 1x)) !Allows 4 dp for rainfall
1635 
1636  ! Deallocate arrays -----------------------------------------------------------------
1637  DEALLOCATE (metfordisagg)
1638  DEALLOCATE (metfordisaggprev)
1639  DEALLOCATE (metfordisaggnext)
1640 
1641  RETURN
1642 
1643 112 CALL errorhint(52, trim(filedscdmet), notused, notused, notusedi)
1644 
1645  ENDSUBROUTINE disaggregatemet
1646  !======================================================================================
1647 
1648  !======================================================================================
1649  SUBROUTINE disaggregateestm(iBlock)
1650  ! Subroutine to disaggregate met forcing data to model time-step
1651  ! HCW 10 Feb 2017
1652  !======================================================================================
1653 
1654  USE sues_data
1655  USE defaultnotused
1656 
1657  IMPLICIT NONE
1658 
1659  INTEGER:: lunit = 101
1660  INTEGER:: tdiff !Time difference (in minutes) between first and second rows of original met forcing file
1661  INTEGER:: i, ii !counter
1662  INTEGER:: iBlock
1663  INTEGER, DIMENSION(NperESTM):: seq1NperESTM
1664  INTEGER, DIMENSION(nsd):: seq1nsd
1665  INTEGER, DIMENSION(ncolsESTMdata):: ESTMDisaggMethod ! Stores method to use for disaggregating met data
1666  REAL(KIND(1d0)), DIMENSION(ncolsESTMdata):: ESTMArrayOrig
1667  REAL(KIND(1d0)), DIMENSION(ReadLinesOrigESTMData*NperESTM, ncolsESTMdata):: ESTM_tt
1668  CHARACTER(LEN=9), DIMENSION(ncolsESTMdata):: HeaderESTM
1669  CHARACTER(LEN=10*ncolsESTMdata):: HeaderESTMOut
1670  ! REAL(KIND(1d0)),DIMENSION(ReadLinesOrigESTMData):: dectimeOrig
1671  ! REAL(KIND(1d0)),DIMENSION(ReadLinesOrigESTMData*NperESTM):: dectimeDscd
1672  INTEGER::iostat_var
1673 
1674  ! INTEGER, DIMENSION(NperESTM):: temp_iy, temp_id, temp_ih, temp_im, temp_ihm
1675 
1676  ! Allocate and initialise arrays to receive original forcing data --------------------
1678  ALLOCATE (estmfordisaggprev(ncolsestmdata))
1679  ALLOCATE (estmfordisaggnext(ncolsestmdata))
1680  estmfordisagg(:, :) = -999
1681  estmfordisaggprev(:) = -999
1682  estmfordisaggnext(:) = -999
1683  ! Initialise array to receive disaggregated data
1684  estm_tt = -999
1685  ! Intialise disaggregation method
1686  estmdisaggmethod = -999
1687 
1688  ! Generate useful sequences
1689  seq1nperestm = (/(i, i=1, nperestm, 1)/)
1690  seq1nsd = (/(i, i=1, nsd, 1)/)
1691 
1692  ! Get methods to use for disaggregation from RunControl
1693  ! (N.B.DisaggMethodESTM is set as 1 or 2 in RunControl; ESTMDisaggMethod is array of ncolsESTMdata used here)
1694  IF (disaggmethodestm == 1) THEN
1695  estmdisaggmethod(:) = 10 !linear disaggregation of averages
1696  ELSEIF (disaggmethodestm == 2) THEN
1697  estmdisaggmethod(:) = 20 !linear disaggregation of instantaneous values
1698  ELSE
1699  CALL errorhint(2, 'Problem in SUEWS_ESTMDisagg: DisaggMethodESTM value should be 1 or 2', &
1701  ENDIF
1702 
1703  ! Read data ---------------------------------------------------------------------
1704  IF (diagnosedisaggestm == 1) WRITE (*, *) 'Reading file: ', trim(fileorigestm)
1705  OPEN (lunit, file=trim(fileorigestm), status='old')
1706  ! CALL skipHeader(lunit,SkipHeaderMet) !Skip header -> read header instead
1707  READ (lunit, *) headerestm
1708  !write(*,*) HeaderMet
1709  ! Skip over lines that have already been read and downscaled
1710  IF (skippedlinesorigestm > 0) THEN
1711  DO i = 1, skippedlinesorigestm - 1 ! minus 1 here because last line of last block needs to be read again
1712  READ (lunit, *)
1713  ENDDO
1714  ! Read in last line of previous block
1715  READ (lunit, *, iostat=iostat_var) estmarrayorig
1716  estmfordisaggprev(1:ncolsestmdata) = estmarrayorig
1717  ENDIF
1718  ! Read in current block
1719  DO i = 1, readlinesorigestmdatamax
1720  READ (lunit, *, iostat=iostat_var) estmarrayorig
1721  estmfordisagg(i, 1:ncolsestmdata) = estmarrayorig
1722  ENDDO
1723  ! Read in first line of next block (except for last block)
1724  IF (iblock /= readblocksorigmetdata) THEN
1725  READ (lunit, *, iostat=iostat_var) estmarrayorig
1726  estmfordisaggnext(1:ncolsestmdata) = estmarrayorig
1727  ENDIF
1728  CLOSE (lunit)
1729 
1730  ! Check resolution of original met forcing data -------------------------------------
1731  ! Find time difference (in minutes) between first and second row
1732  tdiff = int(estmfordisagg(2, 4) - estmfordisagg(1, 4)) !Try using minutes
1733  IF (tdiff == 0) tdiff = int(estmfordisagg(2, 3) - estmfordisagg(1, 3))*60 !If no difference in minutes, try using hours
1734  IF (tdiff < 0) THEN !If time difference is negative (e.g. change of day), instead use second and third row
1735  tdiff = int(estmfordisagg(3, 4) - estmfordisagg(2, 4))
1736  IF (tdiff == 0) tdiff = int(estmfordisagg(3, 3) - estmfordisagg(2, 3))*60 !If no difference in minutes, try using hours
1737  ENDIF
1738  ! Check actual resolution matches specified input resolution
1739  IF (tdiff /= resolutionfilesinestm/60) THEN
1740  CALL errorhint(2, 'Problem in SUEWS_ESTMDisagg: timestamps in ESTM forcing file inconsistent with ResolutionFilesInESTM', &
1741  REAL(ResolutionFilesInESTM, KIND(1d0)), NotUsed, tdiff*60)
1742  ENDIF
1743 
1744  ! Disaggregate time columns ---------------------------------------------------------
1745  IF (diagnose == 1) THEN
1746  WRITE (*, *) 'Disaggregating ESTM forcing data (', trim(fileorigestm), ') to model time-step...'
1747  END IF
1748  CALL disaggregatedatetime(estmfordisagg(:, 1:4), tstep, nperestm, readlinesorigmetdatamax, estm_tt(:, 1:4))
1749  ! ! Convert to dectime
1750  ! dectimeOrig = ESTMForDisagg(:,2) + ESTMForDisagg(:,3)/24.0 + ESTMForDisagg(:,4)/(60.0*24.0)
1751  !
1752  ! DO i=1,ReadLinesOrigESTMDataMax
1753  ! ! Downscale dectime using dectimeOrig(i) [becomes timestamp of last subinterval]
1754  ! dectimeDscd(NperESTM*(i-1)+Seq1NperESTM) = dectimeOrig(i) - (tstep/60.0)/(60.0*24.0)*(/(ii, ii=(NperESTM-1),0, -1)/)
1755  ! ! Convert to required formats
1756  ! temp_iy = INT(ESTMForDisagg(i,1)) !Copy year
1757  ! temp_id = FLOOR(dectimeDscd(NperESTM*(i-1)+Seq1NperESTM)) !DOY
1758  ! ! To avoid precision errors, round here
1759  ! ! - this should not be a problem as a difference of 1 = 1 min, so a difference of 0.001 << 1 min
1760  ! 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)
1761  ! temp_ih = (temp_ihm-MOD(temp_ihm,60))/60 !Hours
1762  ! temp_im = MOD(temp_ihm,60) !Minutes
1763  !
1764  ! IF(dectimeOrig(i) == 1.0000 .AND. i > 1) THEN !If year changes and it is not the beginning of the dataset
1765  ! WRITE(*,*) 'Year change encountered: ', dectimeOrig(i), dectimeOrig(i-1)
1766  ! ! Re-downscale dectime using dectimeOrig(i-1)
1767  ! dectimeDscd(NperESTM*(i-1)+Seq1NperESTM) = dectimeOrig(i-1) + (tstep/60.0)/(60.0*24.0)*Seq1NperESTM
1768  ! ! Convert to required formats
1769  ! temp_iy = INT(ESTMForDisagg(i,1)) !Copy year
1770  ! temp_id = FLOOR(dectimeDscd(NperESTM*(i-1)+Seq1NperESTM)) !DOY
1771  ! 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)
1772  ! temp_ih = (temp_ihm-MOD(temp_ihm,60))/60 !Hours
1773  ! temp_im = MOD(temp_ihm,60) !Minutes
1774  ! ! Adjust year and DOY to account for year change
1775  ! temp_iy(1:(NperESTM-1)) = temp_iy(1:(NperESTM-1)) - 1 !Subtract 1 from year for all except final timestamp
1776  ! temp_id(NperESTM) = 1 !Set day for final timestamp to 1
1777  ! ENDIF
1778  !
1779  ! !IF(i==1 .or. i== ReadlinesOrigESTMDataMax) THEN
1780  ! ! write(*,*) temp_iy
1781  ! ! write(*,*) temp_id
1782  ! ! !write(*,*) temp_ihm
1783  ! ! write(*,*) temp_ih
1784  ! ! write(*,*) temp_im
1785  ! !ENDIF
1786  !
1787  ! ! Copy to ESTM_tt array
1788  ! ESTM_tt(NperESTM*(i-1)+Seq1NperESTM,1) = temp_iy
1789  ! ESTM_tt(NperESTM*(i-1)+Seq1NperESTM,2) = temp_id
1790  ! ESTM_tt(NperESTM*(i-1)+Seq1NperESTM,3) = temp_ih
1791  ! ESTM_tt(NperESTM*(i-1)+Seq1NperESTM,4) = temp_im
1792  !
1793  ! ENDDO
1794 
1795  ! Disaggregate other columns --------------------------------------------------------
1796  ! All other columns are temperatures
1797  DO ii = 5, ncolsestmdata
1798  IF (all(estmfordisagg(:, ii) == -999)) THEN
1799  !IF(DiagnoseDisaggESTM==1) write(*,*) 'No data for col.', ii
1800  estm_tt(:, ii) = -999
1801  ELSE
1802  estm_tt(:, ii) = disagg_lin(estmfordisagg(:, ii), estmfordisaggprev(ii), estmfordisaggnext(ii), estmdisaggmethod(ii), &
1804  ENDIF
1805  ENDDO
1806 
1807  ! Copy disaggregated data to MetForcingDataArray
1809 
1810  ! Write out disaggregated file ------------------------------------------------------
1811  IF (keeptstepfilesin == 1) THEN
1812  IF (iblock == 1) THEN
1813  ! Prepare header
1814  DO i = 1, ncolsestmdata
1815  IF (i == 1) THEN
1816  headerestmout = adjustl(headerestm(i))
1817  ELSE
1818  headerestmout = trim(headerestmout)//' '//adjustl(headerestm(i))
1819  ENDIF
1820  ENDDO
1821  ! Write out header
1822  OPEN (78, file=trim(filedscdestm), err=113)
1823  WRITE (78, '(a)') headerestmout
1824  ELSE
1825  OPEN (78, file=trim(filedscdestm), position='append')!,err=113)
1826  ENDIF
1827  ! Write out data
1829  WRITE (78, 304) (int(estm_tt(i, ii)), ii=1, 4), estm_tt(i, 5:ncolsestmdata)
1830  ENDDO
1831  !IF(iBlock == ReadBlocksOrigMetData) THEN
1832  ! WRITE(78,'(i2)') -9
1833  ! WRITE(78,'(i2)') -9
1834  !ENDIF
1835  CLOSE (78) !Close output file
1836  ENDIF
1837 
1838 304 FORMAT((i4, 1x), 3(i3, 1x), 9(f9.4, 1x))
1839 
1840  ! Deallocate arrays -----------------------------------------------------------------
1841  DEALLOCATE (estmfordisagg)
1842  DEALLOCATE (estmfordisaggprev)
1843  DEALLOCATE (estmfordisaggnext)
1844 
1845  RETURN
1846 
1847 113 CALL errorhint(52, trim(filedscdestm), notused, notused, notusedi)
1848 
1849  ENDSUBROUTINE disaggregateestm
1850  !======================================================================================
1851 
1852  !======================================================================================
1853  SUBROUTINE disaggregatedatetime(DateTimeForDisagg, tstep, Nper, ReadLinesOrigMetDataMax, DateTimeDscd)
1854  USE datetime_module, ONLY: daysinyear
1855  IMPLICIT NONE
1856  INTEGER, INTENT(in) :: tstep, Nper, ReadLinesOrigMetDataMax
1857  REAL(KIND(1d0)), DIMENSION(ReadLinesOrigMetData, 4), INTENT(in):: DateTimeForDisagg
1858  REAL(KIND(1d0)), DIMENSION(ReadLinesOrigMetData*Nper, 4), INTENT(out):: DateTimeDscd
1859 
1860  REAL(KIND(1d0)), DIMENSION(ReadLinesOrigMetData):: dectimeOrig
1861  REAL(KIND(1d0)), DIMENSION(Nper) :: temp_dectime !temorary varaibles for disaggragation
1862  INTEGER, DIMENSION(Nper):: temp_iy, temp_id, temp_ih, temp_im, temp_ihm ! temorary varaibles for disaggragation
1863  INTEGER, DIMENSION(Nper)::ndays_iy ! number of days in iy
1864 
1865  INTEGER :: i, ii
1866  INTEGER, DIMENSION(Nper):: seq1Nper
1867 
1868  ! Generate useful sequences
1869  seq1nper = (/(i, i=1, nper, 1)/)
1870  ! Convert to dectime
1871  ! dectimeOrig = MetForDisagg(:,2) + MetForDisagg(:,3)/24.0 + MetForDisagg(:,4)/(60.0*24.0)
1872  ! 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
1873  dectimeorig = (datetimefordisagg(:, 2) - 1) + datetimefordisagg(:, 3)/24.0 + datetimefordisagg(:, 4)/(60.0*24.0)
1874 
1875  DO i = 1, readlinesorigmetdatamax
1876 
1877  ! Downscale dectime using dectimeOrig(i) [becomes timestamp of last subinterval]
1878  temp_dectime = dectimeorig(i) - (tstep/60.0)/(60.0*24.0)*(/(ii, ii=(nper - 1), 0, -1)/)
1879  temp_dectime = nint(temp_dectime*60*60*24)/(60*60*24*1.)
1880 
1881  ! Convert to required formats
1882  ! get year: year-1 if dectime <0, copy `year` from metforcing otherwise
1883  temp_iy = merge(int(datetimefordisagg(i, 1)) - 1, int(datetimefordisagg(i, 1)), temp_dectime < 0)
1884 
1885  ! get day of year:
1886  ndays_iy = daysinyear(temp_iy) ! get days of year for DOY correction when temp_dectime<0 during year-crossing
1887  temp_dectime = merge(temp_dectime + ndays_iy, temp_dectime, temp_dectime < 0) ! correct minus temp_dectime to positive values
1888  temp_id = floor(temp_dectime) + 1 !DOY
1889 
1890  temp_ihm = nint((temp_dectime + 1 - temp_id/1.0)*60.0*24.0) !Minutes of the day (1440 max)
1891  temp_ih = (temp_ihm - mod(temp_ihm, 60))/60 !Hours
1892  temp_ih = merge(temp_ih, 0, mask=(temp_ih < 24))
1893  temp_im = mod(temp_ihm, 60) !Minutes
1894 
1895  ! Copy to Met_tt array
1896  datetimedscd(nper*(i - 1) + seq1nper, 1) = temp_iy
1897  datetimedscd(nper*(i - 1) + seq1nper, 2) = temp_id
1898  datetimedscd(nper*(i - 1) + seq1nper, 3) = temp_ih
1899  datetimedscd(nper*(i - 1) + seq1nper, 4) = temp_im
1900 
1901  ENDDO
1902 
1903  END SUBROUTINE disaggregatedatetime
1904  !======================================================================================
1905 
1906  ! Define functions here:
1907  FUNCTION disagg_lin(Slow, SlowPrev, SlowNext, DisaggType, Nper_loc, ReadLinesOrig_loc, ReadLinesOrigMax_loc, iBlock) RESULT(Fast)
1909  USE defaultnotused
1910  USE sues_data
1911 
1912  IMPLICIT NONE
1913 
1914  INTEGER:: DisaggType !Type of disaggregation: 10 for averaged variables; 20 for instantaneous variables
1915  INTEGER:: Nper_loc !Number of subintervals per interval (local Nper)
1916  INTEGER:: ReadLinesOrig_loc, ReadLinesOrigMax_loc !Number of lines to read in original file (local)
1917  INTEGER:: iBlock
1918  REAL(KIND(1d0)), DIMENSION(ReadLinesOrig_loc*Nper_loc):: Fast !Array to receive disaggregated data
1919  REAL(KIND(1d0)), DIMENSION(ReadLinesOrig_loc):: Slow !Array to disaggregate
1920  REAL(KIND(1d0)):: SlowPrev, SlowNext
1921  INTEGER, DIMENSION(Nper_loc):: FastRows !Group of rows that are filled with each iteration
1922  INTEGER, DIMENSION(FLOOR(Nper_loc/2.0)):: FirstRows10 !Rows at the beginning that are not filled during iteration (for averages)
1923  INTEGER, DIMENSION(Nper_loc - FLOOR(Nper_loc/2.0)):: LastRows10 !Rows at the end that are not filled during iteration
1924  INTEGER, DIMENSION(Nper_loc):: FirstRows20 !Rows at the beginning that are not filled during iteration (for instantaneous)
1925  INTEGER, DIMENSION(Nper_loc):: seq1Nper_loc !1 to Nper_loc
1926  INTEGER:: XNper_loc !XNper_loc = 2 for even Nper_loc; XNper_loc=1 for odd Nper_loc
1927  INTEGER:: i, ii !counters
1928 
1929  ! Calculate XNper_loc (differentiates between disaggregations with odd and even Nper_loc)
1930  IF (mod(nper_loc, 2) == 0) xnper_loc = 2
1931  IF (mod(nper_loc, 2) == 1) xnper_loc = 1
1932 
1933  seq1nper_loc = (/(i, i=1, nper_loc, 1)/)
1934 
1935  ! Setup counters for iteration
1936  IF (disaggtype == 10) THEN
1937  fastrows = floor(nper_loc/2.0) + seq1nper_loc ! Rows to create at model time-step
1938  firstrows10 = (/(i, i=1, (fastrows(1) - 1), 1)/) !For start of dataset
1939  lastrows10 = &
1940  (/(i, i=nper_loc*(readlinesorigmax_loc - 1 - 1) + fastrows(nper_loc) + 1, &
1941  (readlinesorigmax_loc*nper_loc), 1)/) ! For end of dataset
1942  ELSEIF (disaggtype == 20) THEN
1943  fastrows = nper_loc + seq1nper_loc !Rows to create at model time-step
1944  firstrows20 = (/(i, i=1, (fastrows(1) - 1), 1)/) !For start of dataset
1945  ENDIF
1946 
1947  ! Initialise fast array to -999
1948  fast = -999
1949  ! Linearly disaggregate
1950  IF (disaggtype == 10) THEN !Averaged variables
1951  IF (diagnosedisagg == 1) WRITE (*, *) 'Linearly disaggregating averaged variable'
1952  DO i = 1, (readlinesorigmax_loc - 1)
1953  fast(nper_loc*(i - 1) + fastrows) = slow(i) - &
1954  (slow(i + 1) - slow(i))/(xnper_loc*nper_loc) + &
1955  (slow(i + 1) - slow(i))/nper_loc*(/(ii, ii=1, nper_loc, 1)/)
1956  ENDDO
1957 
1958  ! For first few rows, use previous met block
1959  IF (iblock == 1) THEN
1960  fast(firstrows10) = fast(fastrows(1)) !Use repeat values at the start of the year
1961  ELSE
1962  fast(firstrows10) = slowprev - &
1963  (slow(1) - slowprev)/(xnper_loc*nper_loc) + &
1964  (slow(1) - slowprev)/nper_loc* &
1965  (/(ii, ii=(nper_loc - SIZE(firstrows10) + 1), nper_loc, 1)/)
1966  ENDIF
1967  ! For last few rows, use next met block
1968  IF (iblock == readblocksorigmetdata) THEN
1969  fast(lastrows10) = fast(nper_loc*(readlinesorigmax_loc - 1 - 1) + fastrows(nper_loc)) !Use repeat values at the end of the year
1970  ELSE
1971  fast(lastrows10) = slow(readlinesorigmax_loc) - &
1972  (slownext - slow(readlinesorigmax_loc))/(xnper_loc*nper_loc) + &
1973  (slownext - slow(readlinesorigmax_loc))/nper_loc* &
1974  (/(ii, ii=1, SIZE(lastrows10), 1)/)
1975  ENDIF
1976  ELSEIF (disaggtype == 20) THEN !Instantaneous variables
1977  IF (diagnosedisagg == 1) WRITE (*, *) 'Linearly disaggregating instantaneous variable'
1978  DO i = 1, (readlinesorigmax_loc - 1)
1979  fast(nper_loc*(i - 1) + fastrows) = (slow(i) + &
1980  (slow(i + 1) - slow(i))/nper_loc*2*(seq1nper_loc - 1) + &
1981  slow(i))/2
1982  ENDDO
1983  ! For first few rows, use previous met block
1984  IF (iblock == 1) THEN
1985  fast(firstrows20) = fast(fastrows(1)) !Use repeat values at the start of the year
1986  ELSE
1987  fast(firstrows20) = (slowprev + &
1988  (slow(1) - slowprev)/nper_loc*2* &
1989  ((/(ii, ii=(nper_loc - SIZE(firstrows20) + 1), nper_loc, 1)/) - 1) + &
1990  slowprev)/2
1991  ENDIF
1992  !! Last few rows are already filled for the instantaneous value disaggregation
1993  !IF(iBlock==ReadBlocksOrigMetData) THEN
1994  ! Fast(LastRows20) = Fast(Nper_loc*(ReadLinesOrigMax_loc-1-1)+FastRows(Nper_loc)) !Use repeat values at the end of the year
1995  !ELSE
1996  ! Fast(LastRows20) = (Slow(ReadLinesOrigMax_loc) + &
1997  ! (SlowNext-Slow(ReadLinesOrigMax_loc))/Nper_loc*2 * &
1998  ! ((/(ii, ii=1,SIZE(LastRows20), 1)/)-1) + &
1999  ! Slow(ReadLinesOrigMax_loc))/2
2000  !ENDIF
2001  ENDIF
2002 
2003  IF (any(fast(1:readlinesorigmax_loc*nper_loc) == -999)) THEN
2004  WRITE (*, *) 'Problem: -999s (', count(fast(1:readlinesorigmax_loc*nper_loc) == -999), ') in disaggregated data.'
2005  CALL errorhint(13, 'Problem in SUEWS_MetDisagg: -999 values in disaggregated data.', notused, notused, notusedi)
2006  ENDIF
2007 
2008  ENDFUNCTION disagg_lin
2009  !======================================================================================
2010 
2011  !======================================================================================
2012  FUNCTION disaggp_amongn(Slow, amongN, Nper_loc, ReadLinesOrig_loc, ReadLinesOrigMax_loc) RESULT(Fast)
2013  ! Subroutine to disaggregate precipitation by evenly distributing among N subintervals
2014  ! (i.e. equal intensity in N subintervals)
2015  ! See Ward et al. (in review), meanN, 0.5N or 0.25N approach
2016  ! HCW 10 Feb 2017
2017  !======================================================================================
2018 
2019  USE defaultnotused
2020  USE sues_data
2021 
2022  IMPLICIT NONE
2023 
2024  INTEGER:: amongN !Number of subintervals over which rain will be distributed
2025  INTEGER:: Nper_loc !Number of subintervals per interval (local Nper)
2026  INTEGER:: ReadLinesOrig_loc, ReadLinesOrigMax_loc !Number of lines to read in original file (local)
2027  REAL(KIND(1d0)), DIMENSION(ReadLinesOrig_loc*Nper_loc):: Fast !Array to receive disaggregated data
2028  REAL(KIND(1d0)), DIMENSION(ReadLinesOrig_loc):: Slow !Array to disaggregate
2029  INTEGER, DIMENSION(:), ALLOCATABLE:: Subintervals !Array of subintervals that contain rain
2030  INTEGER, DIMENSION(Nper_loc):: seq1Nper_loc !1 to Nper_loc
2031  INTEGER:: i
2032 
2033  ! For each averaging period, get subintervals which will receive rain
2034  ALLOCATE (subintervals(amongn))
2035  subintervals(:) = -999
2036 
2037  seq1nper_loc = (/(i, i=1, nper_loc, 1)/)
2038 
2039  IF (diagnosedisagg == 1) WRITE (*, *) 'Distributing over ', amongn, ' subintervals for variable'
2040  ! If all subintervals are to contain rain, don't need to generate random numbers
2041  IF (amongn == nper_loc) THEN
2042  subintervals(:) = seq1nper_loc
2043  ENDIF
2044  IF (amongn > nper_loc) &
2045  CALL errorhint(2, 'Problem in SUEWS_MetDisagg: no. of rainy periods cannot exceed number of subintervals', &
2046  REAL(Nper_loc, KIND(1d0)), NotUsed, amongN)
2047 
2048  ! Initialise fast array to -999
2049  fast = -999
2050  DO i = 1, readlinesorigmax_loc
2051  fast(nper_loc*(i - 1) + seq1nper_loc) = 0 !Fill all subintervals with zeros initially
2052  IF (slow(i) > 0) THEN !If there is some rainfall during this interval...
2053  IF (amongn < nper_loc) THEN
2054  subintervals(:) = -999
2055  subintervals = randomsamples(amongn, nper_loc)
2056  ENDIF
2057  fast(nper_loc*(i - 1) + subintervals) = slow(i)/amongn
2058  ENDIF
2059  ENDDO
2060 
2061  IF (any(fast(1:readlinesorigmax_loc*nper_loc) == -999)) THEN
2062  WRITE (*, *) 'Problem: -999s (', count(fast(1:readlinesorigmax_loc*nper_loc) == -999), ') in disaggregated data'
2063  CALL errorhint(13, 'Problem in SUEWS_MetDisagg: -999 values in disaggregated data.', notused, notused, notusedi)
2064  ENDIF
2065 
2066  ENDFUNCTION disaggp_amongn
2067  !======================================================================================
2068 
2069  !======================================================================================
2070  FUNCTION disaggp_amongnmult(Slow, multupperI, multamongN, Nper_loc, ReadLinesOrig_loc, ReadLinesOrigMax_loc) RESULT(Fast)
2071  ! Subroutine to disaggregate precipitation by evenly distributing among N subintervals
2072  ! (i.e. equal intensity in N subintervals) for different intensity bins
2073  ! Based on analsysis by Wen Gu
2074  ! HCW 21 Apr 2017
2075  !======================================================================================
2076 
2077  USE defaultnotused
2078  USE sues_data
2079 
2080  IMPLICIT NONE
2081 
2082  REAL(KIND(1d0)), DIMENSION(5):: multupperI !Upper bound of intensity bin
2083  INTEGER, DIMENSION(5):: multamongN !Number of subintervals over which rain will be distributed (array)
2084  INTEGER:: thisamongN !Number of subintervals over which rain will be distributed
2085  INTEGER:: Nper_loc !Number of subintervals per interval (local Nper)
2086  INTEGER:: ReadLinesOrig_loc, ReadLinesOrigMax_loc !Number of lines to read in original file (local)
2087  REAL(KIND(1d0)), DIMENSION(ReadLinesOrig_loc*Nper_loc):: Fast !Array to receive disaggregated data
2088  REAL(KIND(1d0)), DIMENSION(ReadLinesOrig_loc):: Slow !Array to disaggregate
2089  INTEGER, DIMENSION(:), ALLOCATABLE:: Subintervals !Array of subintervals that contain rain
2090  INTEGER, DIMENSION(Nper_loc):: seq1Nper_loc !1 to Nper_loc
2091  INTEGER:: i
2092 
2093  seq1nper_loc = (/(i, i=1, nper_loc, 1)/)
2094 
2095  IF (diagnosedisagg == 1) WRITE (*, *) 'Distributing over variable subintervals depending on intensity for variable'
2096 
2097  ! Initialise fast array to -999
2098  fast = -999
2099  DO i = 1, readlinesorigmax_loc
2100  fast(nper_loc*(i - 1) + seq1nper_loc) = 0 !Fill all subintervals with zeros initially
2101  IF (slow(i) > 0) THEN !If there is some rainfall during this interval...
2102  !Use intensity in this interval to decide number of subintervals to fill with rain
2103  IF (slow(i) <= multupperi(1)) THEN
2104  thisamongn = multamongn(1)
2105  ELSEIF (slow(i) > multupperi(1) .AND. slow(i) <= multupperi(2)) THEN
2106  thisamongn = multamongn(2)
2107  ELSEIF (slow(i) > multupperi(2) .AND. slow(i) <= multupperi(3)) THEN
2108  thisamongn = multamongn(3)
2109  ELSEIF (slow(i) > multupperi(3) .AND. slow(i) <= multupperi(4)) THEN
2110  thisamongn = multamongn(4)
2111  ELSEIF (slow(i) > multupperi(4) .AND. slow(i) <= multupperi(5)) THEN
2112  thisamongn = multamongn(5)
2113  ELSEIF (slow(i) > multupperi(5)) THEN
2114  thisamongn = multamongn(5)
2115  CALL errorhint(4, 'Precip in met forcing file exceeds maxiumum MultRainAmongNUpperI', &
2116  slow(i), multrainamongnupperi(5), notusedi)
2117  ENDIF
2118 
2119  ! For each averaging period, get subintervals which will receive rain
2120  ALLOCATE (subintervals(thisamongn))
2121  subintervals(:) = -999
2122 
2123  IF (thisamongn > nper_loc) CALL errorhint(2, 'Problem in SUEWS_MetDisagg: no. of rainy periods cannot exceed '// &
2124  'number of subintervals', REAL(Nper_loc, KIND(1d0)), NotUsed, thisamongN)
2125 
2126  IF (thisamongn == nper_loc) THEN ! If all subintervals are to contain rain, don't need to generate random numbers
2127  subintervals(:) = seq1nper_loc
2128  ELSEIF (thisamongn < nper_loc) THEN
2129  subintervals = randomsamples(thisamongn, nper_loc)
2130  ENDIF
2131  fast(nper_loc*(i - 1) + subintervals) = slow(i)/thisamongn
2132  !write(*,*) Slow(i), thisamongN
2133  DEALLOCATE (subintervals)
2134  ENDIF
2135  ENDDO
2136 
2137  IF (any(fast(1:readlinesorigmax_loc*nper_loc) == -999)) THEN
2138  WRITE (*, *) 'Problem: -999s (', count(fast(1:readlinesorigmax_loc*nper_loc) == -999), ') in disaggregated data'
2139  CALL errorhint(13, 'Problem in SUEWS_MetDisagg: -999 values in disaggregated data.', notused, notused, notusedi)
2140  ENDIF
2141 
2142  ENDFUNCTION disaggp_amongnmult
2143  !======================================================================================
2144 
2145  !======================================================================================
2146  FUNCTION randomsamples(N, OutOf) RESULT(Samples)
2147  ! Generates N/OutOf random samples without repeats
2148  ! e.g. for N = 3 and OutOf = 12, a possibility for Samples = 7,3,11
2149  ! HCW 10 Feb 2017
2150  !======================================================================================
2151 
2152  IMPLICIT NONE
2153 
2154  INTEGER:: i !counter
2155  INTEGER:: N !number of samples to return
2156  INTEGER:: OutOf !number to sample from
2157  INTEGER:: X !next sample to be added
2158  REAL(KIND(1D0)):: r !random number
2159  INTEGER, DIMENSION(:), ALLOCATABLE:: Samples !Array to receive random samples
2160 
2161  ! Allocate and initialise Samples
2162  ALLOCATE (samples(n))
2163  samples(:) = -999
2164 
2165  ! Generate random sample (no repeats)
2166  i = 0 !Set counter to zero initially
2167  DO WHILE (any(samples == -999))
2168  CALL random_number(r)
2169  x = int(r*outof) + 1
2170  !write(*,*) X
2171  !write(*,*) COUNT(Samples == X)
2172  IF (count(samples == x) == 0) THEN
2173  ! Only keep if this subinterval has not already been selected
2174  i = i + 1
2175  samples(i) = x
2176  ENDIF
2177  !write(*,*) Samples
2178  ENDDO
2179 
2180  ENDFUNCTION randomsamples
2181  !======================================================================================
2182 
2183 ENDMODULE metdisagg
2184 !========================================================================================
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