SUEWS API Site
Documentation of SUEWS source code
suews_ctrl_driver.f95
Go to the documentation of this file.
1 !========================================================================================
2 ! SUEWS driver subroutines
3 ! TS 31 Aug 2017: initial version
4 ! TS 02 Oct 2017: added `SUEWS_cal_Main` as the generic wrapper
5 ! TS 03 Oct 2017: added `SUEWS_cal_AnthropogenicEmission`
7  USE meteo, ONLY: qsatf, rh2qa, qa2rh
10  USE anohm_module, ONLY: anohm
13  USE estm_module, ONLY: estm
16  USE waterdist_module, ONLY: drainage, soilstore, &
20  USE ctrl_output, ONLY: varlistall
22  use lumps_module, only: lumps_cal_qhqe
23  use evap_module, only: evap_suews
24  use rsl_module, only: rslprofile
26  use co2_module, only: co2_biogen
27  use evap_module, only: evap_suews
28  USE allocatearray, ONLY: &
29  nsurf, nvegsurf, &
35  use moist, only: avcp, avdens, lv_j_kg
36 
37  IMPLICIT NONE
38 
39 CONTAINS
40  ! ===================MAIN CALCULATION WRAPPER FOR ENERGY AND WATER FLUX===========
41  SUBROUTINE suews_cal_main( &
42  AerodynamicResistanceMethod, AH_MIN, AHProf_24hr, AH_SLOPE_Cooling, & ! input&inout in alphabetical order
43  AH_SLOPE_Heating, &
44  alb, AlbMax_DecTr, AlbMax_EveTr, AlbMax_Grass, &
45  AlbMin_DecTr, AlbMin_EveTr, AlbMin_Grass, &
46  alpha_bioCO2, alpha_enh_bioCO2, alt, avkdn, avRh, avU1, BaseT, BaseTe, &
47  BaseTHDD, beta_bioCO2, beta_enh_bioCO2, bldgH, CapMax_dec, CapMin_dec, &
48  chAnOHM, CO2PointSource, cpAnOHM, CRWmax, CRWmin, DayWat, DayWatPer, &
49  DecTreeH, Diagnose, DiagQN, DiagQS, DRAINRT, &
50  dt_since_start, dqndt, qn1_av, dqnsdt, qn1_s_av, &
51  EF_umolCO2perJ, emis, EmissionsMethod, EnEF_v_Jkm, endDLS, EveTreeH, FAIBldg, &
52  FAIDecTree, FAIEveTree, Faut, FcEF_v_kgkm, fcld_obs, FlowChange, &
53  FrFossilFuel_Heat, FrFossilFuel_NonHeat, G1, G2, G3, G4, G5, G6, GDD_id, &
54  GDDFull, Gridiv, gsModel, HDD_id, HumActivity_24hr, &
55  IceFrac, id, Ie_a, Ie_end, Ie_m, Ie_start, imin, &
56  InternalWaterUse_h, IrrFracConif, IrrFracDecid, IrrFracGrass, isec, it, EvapMethod, &
57  iy, kkAnOHM, Kmax, LAI_id, LAICalcYes, LAIMax, LAIMin, LAI_obs, &
58  LAIPower, LAIType, lat, lenDay_id, ldown_obs, lng, MaxConductance, MaxFCMetab, MaxQFMetab, &
59  SnowWater, MetForcingData_grid, MinFCMetab, MinQFMetab, min_res_bioCO2, &
60  NARP_EMIS_SNOW, NARP_TRANS_SITE, NetRadiationMethod, &
61  OHM_coef, OHMIncQF, OHM_threshSW, &
62  OHM_threshWD, PipeCapacity, PopDensDaytime, &
63  PopDensNighttime, PopProf_24hr, PorMax_dec, PorMin_dec, &
64  Precip, PrecipLimit, PrecipLimitAlb, Press_hPa, &
65  QF0_BEU, Qf_A, Qf_B, Qf_C, &
66  qn1_obs, qs_obs, qf_obs, &
67  RadMeltFact, RAINCOVER, RainMaxRes, resp_a, resp_b, &
68  RoughLenHeatMethod, RoughLenMomMethod, RunoffToWater, S1, S2, &
69  SatHydraulicConduct, SDDFull, SDD_id, sfr, SMDMethod, SnowAlb, SnowAlbMax, &
70  SnowAlbMin, SnowPackLimit, SnowDens, SnowDensMax, SnowDensMin, SnowfallCum, SnowFrac, &
71  SnowLimBldg, SnowLimPaved, snowFrac_obs, SnowPack, SnowProf_24hr, snowUse, SoilDepth, &
72  soilstore_id, SoilStoreCap, StabilityMethod, startDLS, state_id, StateLimit, &
73  StorageHeatMethod, StoreDrainPrm, SurfaceArea, Tair_av, tau_a, tau_f, tau_r, &
74  Tmax_id, Tmin_id, &
75  T_CRITIC_Cooling, T_CRITIC_Heating, Temp_C, TempMeltFact, TH, &
76  theta_bioCO2, timezone, TL, TrafficRate, TrafficUnits, &
77  TraffProf_24hr, Ts5mindata_ir, tstep, tstep_prev, veg_type, &
78  WaterDist, WaterUseMethod, WetThresh, wu_m3, &
79  WUDay_id, DecidCap_id, albDecTr_id, albEveTr_id, albGrass_id, porosity_id, &
80  WUProfA_24hr, WUProfM_24hr, xsmd, Z, z0m_in, zdm_in, &
81  datetimeLine, dataOutLineSUEWS, dataOutLineSnow, dataOutLineESTM, dataoutLineRSL, &!output
82  DailyStateLine)!output
83 
84  IMPLICIT NONE
85 
86  ! ########################################################################################
87  ! input variables
88  INTEGER, INTENT(IN)::AerodynamicResistanceMethod
89  INTEGER, INTENT(IN)::Diagnose
90  INTEGER, INTENT(IN)::DiagQN
91  INTEGER, INTENT(IN)::DiagQS
92  INTEGER, INTENT(IN)::startDLS
93  INTEGER, INTENT(IN)::endDLS
94  INTEGER, INTENT(IN)::EmissionsMethod
95  INTEGER, INTENT(IN)::Gridiv
96  INTEGER, INTENT(IN)::gsModel
97  INTEGER, INTENT(IN)::id
98  INTEGER, INTENT(IN)::Ie_end
99  INTEGER, INTENT(IN)::Ie_start
100  INTEGER, INTENT(IN)::isec
101  INTEGER, INTENT(IN)::imin
102  INTEGER, INTENT(IN)::it
103  INTEGER, INTENT(IN)::EvapMethod
104  INTEGER, INTENT(IN)::iy
105  INTEGER, INTENT(IN)::LAICalcYes
106  INTEGER, INTENT(IN)::NetRadiationMethod
107  INTEGER, INTENT(IN)::OHMIncQF
108  INTEGER, INTENT(IN)::RoughLenHeatMethod
109  INTEGER, INTENT(IN)::RoughLenMomMethod
110  INTEGER, INTENT(IN)::SMDMethod
111  INTEGER, INTENT(IN)::snowUse
112  INTEGER, INTENT(IN)::StabilityMethod
113  INTEGER, INTENT(IN)::StorageHeatMethod
114  INTEGER, INTENT(IN)::tstep
115  INTEGER, INTENT(IN)::tstep_prev ! tstep size of the previous step
116  INTEGER, INTENT(in)::dt_since_start ! time since simulation starts [s]
117  INTEGER, INTENT(IN)::veg_type
118  INTEGER, INTENT(IN)::WaterUseMethod
119 
120  REAL(KIND(1D0)), INTENT(IN)::AlbMax_DecTr
121  REAL(KIND(1D0)), INTENT(IN)::AlbMax_EveTr
122  REAL(KIND(1D0)), INTENT(IN)::AlbMax_Grass
123  REAL(KIND(1D0)), INTENT(IN)::AlbMin_DecTr
124  REAL(KIND(1D0)), INTENT(IN)::AlbMin_EveTr
125  REAL(KIND(1D0)), INTENT(IN)::AlbMin_Grass
126  REAL(KIND(1D0)), INTENT(IN)::alt
127  REAL(KIND(1D0)), INTENT(IN)::avkdn
128  REAL(KIND(1D0)), INTENT(IN)::avRh
129  REAL(KIND(1D0)), INTENT(IN)::avU1
130  REAL(KIND(1D0)), INTENT(IN)::BaseTHDD
131  REAL(KIND(1D0)), INTENT(IN)::bldgH
132  REAL(KIND(1D0)), INTENT(IN)::CapMax_dec
133  REAL(KIND(1D0)), INTENT(IN)::CapMin_dec
134  REAL(KIND(1D0)), INTENT(IN)::CO2PointSource
135  REAL(KIND(1D0)), INTENT(IN)::CRWmax
136  REAL(KIND(1D0)), INTENT(IN)::CRWmin
137  REAL(KIND(1D0)), INTENT(IN)::DecTreeH
138  REAL(KIND(1D0)), INTENT(IN)::DRAINRT
139  REAL(KIND(1D0)), INTENT(IN)::EF_umolCO2perJ
140  REAL(KIND(1D0)), INTENT(IN)::EnEF_v_Jkm
141  REAL(KIND(1D0)), INTENT(IN)::EveTreeH
142  REAL(KIND(1D0)), INTENT(IN)::FAIBldg
143  REAL(KIND(1D0)), INTENT(IN)::FAIDecTree
144  REAL(KIND(1D0)), INTENT(IN)::FAIEveTree
145  REAL(KIND(1D0)), INTENT(IN)::Faut
146  REAL(KIND(1D0)), INTENT(IN)::fcld_obs
147  REAL(KIND(1D0)), INTENT(IN)::FlowChange
148  REAL(KIND(1D0)), INTENT(IN)::FrFossilFuel_Heat
149  REAL(KIND(1D0)), INTENT(IN)::FrFossilFuel_NonHeat
150  REAL(KIND(1D0)), INTENT(IN)::G1
151  REAL(KIND(1D0)), INTENT(IN)::G2
152  REAL(KIND(1D0)), INTENT(IN)::G3
153  REAL(KIND(1D0)), INTENT(IN)::G4
154  REAL(KIND(1D0)), INTENT(IN)::G5
155  REAL(KIND(1D0)), INTENT(IN)::G6
156  REAL(KIND(1D0)), INTENT(IN)::InternalWaterUse_h
157  REAL(KIND(1D0)), INTENT(IN)::IrrFracConif
158  REAL(KIND(1D0)), INTENT(IN)::IrrFracDecid
159  REAL(KIND(1D0)), INTENT(IN)::IrrFracGrass
160  REAL(KIND(1D0)), INTENT(IN)::Kmax
161  REAL(KIND(1D0)), INTENT(IN)::LAI_obs
162  REAL(KIND(1D0)), INTENT(IN)::lat
163  REAL(KIND(1D0)), INTENT(IN)::ldown_obs
164  REAL(KIND(1D0)), INTENT(IN)::lng
165  REAL(KIND(1D0)), INTENT(IN)::MaxFCMetab
166  REAL(KIND(1D0)), INTENT(IN)::MaxQFMetab
167  REAL(KIND(1D0)), INTENT(IN)::MinFCMetab
168  REAL(KIND(1D0)), INTENT(IN)::MinQFMetab
169  REAL(KIND(1D0)), INTENT(IN)::NARP_EMIS_SNOW
170  REAL(KIND(1D0)), INTENT(IN)::NARP_TRANS_SITE
171  REAL(KIND(1D0)), INTENT(IN)::PipeCapacity
172  REAL(KIND(1D0)), INTENT(IN)::PopDensNighttime
173  REAL(KIND(1D0)), INTENT(IN)::PorMax_dec
174  REAL(KIND(1D0)), INTENT(IN)::PorMin_dec
175  REAL(KIND(1D0)), INTENT(IN)::Precip
176  REAL(KIND(1D0)), INTENT(IN)::PrecipLimit
177  REAL(KIND(1D0)), INTENT(IN)::PrecipLimitAlb
178  REAL(KIND(1D0)), INTENT(IN)::Press_hPa
179  REAL(KIND(1D0)), INTENT(IN)::qn1_obs
180  REAL(KIND(1D0)), INTENT(IN)::qs_obs
181  REAL(KIND(1D0)), INTENT(IN)::qf_obs
182  REAL(KIND(1D0)), INTENT(IN)::RadMeltFact
183  REAL(KIND(1D0)), INTENT(IN)::RAINCOVER
184  REAL(KIND(1D0)), INTENT(IN)::RainMaxRes
185  REAL(KIND(1D0)), INTENT(IN)::RunoffToWater
186  REAL(KIND(1D0)), INTENT(IN)::S1
187  REAL(KIND(1D0)), INTENT(IN)::S2
188  REAL(KIND(1D0)), INTENT(IN)::SnowAlbMax
189  REAL(KIND(1D0)), INTENT(IN)::SnowAlbMin
190  REAL(KIND(1D0)), INTENT(IN)::SnowDensMax
191  REAL(KIND(1D0)), INTENT(IN)::SnowDensMin
192  REAL(KIND(1D0)), INTENT(IN)::SnowLimBldg
193  REAL(KIND(1D0)), INTENT(IN)::SnowLimPaved
194  REAL(KIND(1D0)), INTENT(IN)::snowFrac_obs
195  REAL(KIND(1D0)), INTENT(IN)::SurfaceArea
196  REAL(KIND(1D0)), INTENT(IN)::tau_a
197  REAL(KIND(1D0)), INTENT(IN)::tau_f
198  REAL(KIND(1D0)), INTENT(IN)::tau_r
199  REAL(KIND(1D0)), INTENT(IN)::Temp_C
200  REAL(KIND(1D0)), INTENT(IN)::TempMeltFact
201  REAL(KIND(1D0)), INTENT(IN)::TH
202  REAL(KIND(1D0)), INTENT(IN)::timezone
203  REAL(KIND(1D0)), INTENT(IN)::TL
204  REAL(KIND(1D0)), INTENT(IN)::TrafficUnits
205  REAL(KIND(1D0)), INTENT(IN)::wu_m3
206  REAL(KIND(1D0)), INTENT(IN)::xsmd
207  REAL(KIND(1D0)), INTENT(IN)::Z
208  REAL(KIND(1D0)), INTENT(IN)::z0m_in
209  REAL(KIND(1D0)), INTENT(IN)::zdm_in
210 
211  INTEGER, DIMENSION(NVEGSURF), INTENT(IN)::LAIType
212 
213  REAL(KIND(1D0)), DIMENSION(2), INTENT(IN) ::AH_MIN
214  REAL(KIND(1D0)), DIMENSION(2), INTENT(IN) ::AH_SLOPE_Cooling
215  REAL(KIND(1D0)), DIMENSION(2), INTENT(IN) ::AH_SLOPE_Heating
216  REAL(KIND(1D0)), DIMENSION(2), INTENT(IN) ::FcEF_v_kgkm
217  REAL(KIND(1D0)), DIMENSION(2), INTENT(IN) ::QF0_BEU
218  REAL(KIND(1D0)), DIMENSION(2), INTENT(IN) ::Qf_A
219  REAL(KIND(1D0)), DIMENSION(2), INTENT(IN) ::Qf_B
220  REAL(KIND(1D0)), DIMENSION(2), INTENT(IN) ::Qf_C
221  REAL(KIND(1D0)), DIMENSION(2), INTENT(IN) ::PopDensDaytime
222  REAL(KIND(1D0)), DIMENSION(2), INTENT(IN) ::T_CRITIC_Cooling
223  REAL(KIND(1D0)), DIMENSION(2), INTENT(IN) ::T_CRITIC_Heating
224  REAL(KIND(1D0)), DIMENSION(2), INTENT(IN) ::TrafficRate
225  REAL(KIND(1D0)), DIMENSION(3), INTENT(IN) ::Ie_a
226  REAL(KIND(1D0)), DIMENSION(3), INTENT(IN) ::Ie_m
227  REAL(KIND(1D0)), DIMENSION(3), INTENT(IN) ::MaxConductance
228  REAL(KIND(1D0)), DIMENSION(7), INTENT(IN) ::DayWat
229  REAL(KIND(1D0)), DIMENSION(7), INTENT(IN) ::DayWatPer
230  REAL(KIND(1D0)), DIMENSION(nsurf + 1), INTENT(IN) ::OHM_threshSW
231  REAL(KIND(1D0)), DIMENSION(nsurf + 1), INTENT(IN) ::OHM_threshWD
232  REAL(KIND(1D0)), DIMENSION(NSURF), INTENT(IN) ::chAnOHM
233  REAL(KIND(1D0)), DIMENSION(NSURF), INTENT(IN) ::cpAnOHM
234  REAL(KIND(1D0)), DIMENSION(NSURF), INTENT(IN) ::emis
235  REAL(KIND(1D0)), DIMENSION(NSURF), INTENT(IN) ::kkAnOHM
236  REAL(KIND(1D0)), DIMENSION(NSURF), INTENT(IN) ::SatHydraulicConduct
237  REAL(KIND(1D0)), DIMENSION(NSURF), INTENT(IN) ::sfr
238  REAL(KIND(1D0)), DIMENSION(NSURF), INTENT(IN) ::SnowPackLimit
239  REAL(KIND(1D0)), DIMENSION(NSURF), INTENT(IN) ::SoilDepth
240  REAL(KIND(1D0)), DIMENSION(NSURF), INTENT(IN) ::SoilStoreCap
241  REAL(KIND(1D0)), DIMENSION(NSURF), INTENT(IN) ::StateLimit
242  REAL(KIND(1D0)), DIMENSION(NSURF), INTENT(IN) ::WetThresh
243  REAL(KIND(1D0)), DIMENSION(NVEGSURF), INTENT(IN) ::alpha_bioCO2
244  REAL(KIND(1D0)), DIMENSION(NVEGSURF), INTENT(IN) ::alpha_enh_bioCO2
245  REAL(KIND(1D0)), DIMENSION(NVEGSURF), INTENT(IN) ::BaseT
246  REAL(KIND(1D0)), DIMENSION(NVEGSURF), INTENT(IN) ::BaseTe
247  REAL(KIND(1D0)), DIMENSION(NVEGSURF), INTENT(IN) ::beta_bioCO2
248  REAL(KIND(1D0)), DIMENSION(NVEGSURF), INTENT(IN) ::beta_enh_bioCO2
249  REAL(KIND(1D0)), DIMENSION(NVEGSURF), INTENT(IN) ::GDDFull
250  REAL(KIND(1D0)), DIMENSION(NVEGSURF), INTENT(IN) ::LAIMax
251  REAL(KIND(1D0)), DIMENSION(NVEGSURF), INTENT(IN) ::LAIMin
252  REAL(KIND(1D0)), DIMENSION(NVEGSURF), INTENT(IN) ::min_res_bioCO2
253  REAL(KIND(1D0)), DIMENSION(NVEGSURF), INTENT(IN) ::resp_a
254  REAL(KIND(1D0)), DIMENSION(NVEGSURF), INTENT(IN) ::resp_b
255  REAL(KIND(1D0)), DIMENSION(NVEGSURF), INTENT(IN) ::SDDFull
256  REAL(KIND(1D0)), DIMENSION(0:23, 2), INTENT(IN) ::SnowProf_24hr
257  REAL(KIND(1D0)), DIMENSION(NVEGSURF), INTENT(IN) ::theta_bioCO2
258  REAL(KIND(1D0)), DIMENSION(4, NVEGSURF), INTENT(IN) ::LAIPower
259  REAL(KIND(1D0)), DIMENSION(nsurf + 1, 4, 3), INTENT(IN) ::OHM_coef
260  REAL(KIND(1D0)), DIMENSION(NSURF + 1, NSURF - 1), INTENT(IN) ::WaterDist
261  REAL(KIND(1d0)), DIMENSION(:), INTENT(IN) ::Ts5mindata_ir
262  REAL(KIND(1D0)), DIMENSION(:, :), INTENT(IN) ::MetForcingData_grid
263 
264  ! diurnal profile values for 24hr
265  REAL(KIND(1D0)), DIMENSION(0:23, 2), INTENT(IN) ::AHProf_24hr
266  REAL(KIND(1D0)), DIMENSION(0:23, 2), INTENT(IN) ::HumActivity_24hr
267  REAL(KIND(1D0)), DIMENSION(0:23, 2), INTENT(IN) ::PopProf_24hr
268  REAL(KIND(1D0)), DIMENSION(0:23, 2), INTENT(IN) ::TraffProf_24hr
269  REAL(KIND(1D0)), DIMENSION(0:23, 2), INTENT(IN) ::WUProfA_24hr
270  REAL(KIND(1D0)), DIMENSION(0:23, 2), INTENT(IN) ::WUProfM_24hr
271 
272  ! ########################################################################################
273 
274  ! ########################################################################################
275  ! inout variables
276  ! OHM related:
277  REAL(KIND(1d0)), INTENT(INOUT) ::qn1_av
278  REAL(KIND(1d0)), INTENT(INOUT) ::dqndt
279  REAL(KIND(1d0)), INTENT(INOUT) ::qn1_s_av
280  REAL(KIND(1d0)), INTENT(INOUT) ::dqnsdt
281 
282  ! snow related:
283  REAL(KIND(1D0)), INTENT(INOUT) ::SnowfallCum
284  REAL(KIND(1D0)), INTENT(INOUT) ::SnowAlb
285  REAL(KIND(1D0)), DIMENSION(NSURF), INTENT(INOUT)::IceFrac
286  REAL(KIND(1D0)), DIMENSION(NSURF), INTENT(INOUT)::SnowWater
287  REAL(KIND(1D0)), DIMENSION(NSURF), INTENT(INOUT)::SnowDens
288  REAL(KIND(1D0)), DIMENSION(NSURF), INTENT(INOUT)::SnowFrac
289  REAL(KIND(1D0)), DIMENSION(NSURF), INTENT(INOUT)::SnowPack
290 
291  ! water balance related:
292  REAL(KIND(1D0)), DIMENSION(NSURF), INTENT(INOUT) ::soilstore_id
293  REAL(KIND(1D0)), DIMENSION(NSURF), INTENT(INOUT) ::state_id
294  REAL(KIND(1D0)), DIMENSION(6, NSURF), INTENT(INOUT)::StoreDrainPrm
295 
296  ! phenology related:
297  REAL(KIND(1D0)), DIMENSION(NSURF), INTENT(INOUT) ::alb
298  REAL(KIND(1d0)), DIMENSION(nvegsurf), INTENT(INOUT)::GDD_id !Growing Degree Days (see SUEWS_DailyState.f95)
299  REAL(KIND(1d0)), DIMENSION(nvegsurf), INTENT(INout)::SDD_id !Senescence Degree Days (see SUEWS_DailyState.f95)
300  REAL(KIND(1d0)), DIMENSION(nvegsurf), INTENT(INOUT)::LAI_id !LAI for each veg surface [m2 m-2]
301  REAL(KIND(1d0)), INTENT(INout)::Tmin_id
302  REAL(KIND(1d0)), INTENT(INout)::Tmax_id
303  REAL(KIND(1d0)), INTENT(INout)::lenDay_id
304  REAL(KIND(1d0)), INTENT(INOUT)::DecidCap_id
305  REAL(KIND(1d0)), INTENT(INOUT)::albDecTr_id
306  REAL(KIND(1d0)), INTENT(INOUT)::albEveTr_id
307  REAL(KIND(1d0)), INTENT(INOUT)::albGrass_id
308  REAL(KIND(1d0)), INTENT(INOUT)::porosity_id
309 
310  ! anthropogenic heat related:
311  REAL(KIND(1d0)), DIMENSION(12), INTENT(INOUT) ::HDD_id !Heating Degree Days (see SUEWS_DailyState.f95)
312 
313  ! water use related:
314  REAL(KIND(1d0)), DIMENSION(9), INTENT(INOUT) ::WUDay_id
315 
316  ! ESTM related:
317  REAL(KIND(1d0)), INTENT(INOUT) ::Tair_av
318  ! ########################################################################################
319 
320  ! ########################################################################################
321  ! output variables
322  REAL(KIND(1D0)), DIMENSION(5), INTENT(OUT) ::datetimeLine
323  REAL(KIND(1D0)), DIMENSION(ncolumnsDataOutSUEWS - 5), INTENT(OUT) ::dataOutLineSUEWS
324  REAL(KIND(1D0)), DIMENSION(ncolumnsDataOutSnow - 5), INTENT(OUT) ::dataOutLineSnow
325  REAL(KIND(1d0)), DIMENSION(ncolumnsDataOutESTM - 5), INTENT(OUT) ::dataOutLineESTM
326  REAL(KIND(1d0)), DIMENSION(ncolumnsDataOutRSL - 5), INTENT(OUT) ::dataoutLineRSL ! RSL variable array
327  REAL(KIND(1d0)), DIMENSION(ncolumnsDataOutDailyState - 5), INTENT(OUT)::DailyStateLine
328  ! ########################################################################################
329 
330  ! ########################################################################################
331  ! local variables
332  REAL(KIND(1D0))::a1
333  REAL(KIND(1D0))::a2
334  REAL(KIND(1D0))::a3
335  REAL(KIND(1D0))::AdditionalWater
336  REAL(KIND(1D0))::U10_ms
337  REAL(KIND(1D0))::azimuth
338  REAL(KIND(1D0))::chSnow_per_interval
339 
340  REAL(KIND(1D0))::dens_dry
341  REAL(KIND(1d0))::deltaLAI
342  REAL(KIND(1D0))::drain_per_tstep
343  REAL(KIND(1D0))::Ea_hPa
344  REAL(KIND(1D0))::QE_LUMPS
345  REAL(KIND(1D0))::es_hPa
346  REAL(KIND(1D0))::ev_per_tstep
347  REAL(KIND(1D0))::ext_wu
348  REAL(KIND(1D0))::Fc
349  REAL(KIND(1D0))::Fc_anthro
350  REAL(KIND(1D0))::Fc_biogen
351  REAL(KIND(1D0))::Fc_build
352  REAL(KIND(1D0))::fcld
353  REAL(KIND(1D0))::Fc_metab
354  REAL(KIND(1D0))::Fc_photo
355  REAL(KIND(1D0))::Fc_point
356  REAL(KIND(1D0))::Fc_respi
357  REAL(KIND(1D0))::Fc_traff
358  REAL(KIND(1D0))::gfunc
359  REAL(KIND(1D0))::gsc
360  REAL(KIND(1D0))::QH_LUMPS
361  REAL(KIND(1D0))::int_wu
362  REAL(KIND(1D0))::kclear
363  REAL(KIND(1D0))::kup
364  REAL(KIND(1D0))::ldown
365  REAL(KIND(1D0))::lup
366  REAL(KIND(1D0))::L_mod
367  REAL(KIND(1D0))::mwh
368  REAL(KIND(1D0))::mwstore
369  REAL(KIND(1D0))::NWstate_per_tstep
370  REAL(KIND(1D0))::planF
371  REAL(KIND(1D0))::zL
372  REAL(KIND(1D0))::q2_gkg
373  REAL(KIND(1D0))::qe
374  REAL(KIND(1D0))::qf
375  REAL(KIND(1D0))::QF_SAHP
376  REAL(KIND(1D0))::qh
377  REAL(KIND(1D0))::qh_residual
378  REAL(KIND(1D0))::qh_resist
379  REAL(KIND(1D0))::Qm
380  REAL(KIND(1D0))::QmFreez
381  REAL(KIND(1D0))::QmRain
382  REAL(KIND(1D0))::qn1
383  REAL(KIND(1D0))::qn1_S
384  REAL(KIND(1D0))::qn1_snowfree
385  REAL(KIND(1D0))::qs
386  REAL(KIND(1D0))::RA
387  REAL(KIND(1D0))::ResistSurf
388  REAL(KIND(1D0))::RH2
389  REAL(KIND(1d0))::runoffAGveg
390  REAL(KIND(1d0))::runoffAGimpervious
391  REAL(KIND(1D0))::runoff_per_tstep
392  REAL(KIND(1D0))::runoffPipes
393  REAL(KIND(1D0))::runoffSoil_per_tstep
394  REAL(KIND(1D0))::runoffwaterbody
395  REAL(KIND(1D0))::smd
396  REAL(KIND(1D0))::SoilState
397  REAL(KIND(1D0))::state_per_tstep
398  REAL(KIND(1D0))::surf_chang_per_tstep
399  REAL(KIND(1D0))::swe
400  REAL(KIND(1D0))::t2_C
401  REAL(KIND(1D0))::TSfc_C
402  REAL(KIND(1D0))::TempVeg
403  REAL(KIND(1D0))::tot_chang_per_tstep
404  REAL(KIND(1D0))::TStar
405  REAL(KIND(1D0))::tsurf
406  REAL(KIND(1D0))::UStar
407  REAL(KIND(1D0))::VPD_Pa
408  REAL(KIND(1D0))::wu_DecTr
409  REAL(KIND(1D0))::wu_EveTr
410  REAL(KIND(1D0))::wu_Grass
411  REAL(KIND(1D0))::z0m
412  REAL(KIND(1D0))::zdm
413  REAL(KIND(1D0))::ZENITH_deg
414  REAL(KIND(1D0))::Zh
415 
416  REAL(KIND(1D0)), DIMENSION(2)::SnowRemoval
417  REAL(KIND(1D0)), DIMENSION(NSURF)::FreezMelt
418  REAL(KIND(1d0)), DIMENSION(nsurf)::kup_ind_snow
419  REAL(KIND(1D0)), DIMENSION(NSURF)::mw_ind
420  REAL(KIND(1D0)), DIMENSION(NSURF)::Qm_freezState
421  REAL(KIND(1D0)), DIMENSION(NSURF)::Qm_melt
422  REAL(KIND(1D0)), DIMENSION(NSURF)::Qm_rain
423  REAL(KIND(1D0)), DIMENSION(NSURF)::qn1_ind_snow
424  REAL(KIND(1D0)), DIMENSION(NSURF)::rainOnSnow
425  REAL(KIND(1D0)), DIMENSION(NSURF)::runoffSoil
426  REAL(KIND(1D0)), DIMENSION(NSURF)::smd_nsurf
427  REAL(KIND(1D0)), DIMENSION(NSURF)::snowDepth
428 
429  REAL(KIND(1d0)), DIMENSION(nsurf)::Tsurf_ind_snow
430 
431  INTEGER, DIMENSION(NSURF)::snowCalcSwitch
432  INTEGER, DIMENSION(3) ::dayofWeek_id
433  INTEGER::DLS
434 
435  ! REAL(KIND(1D0))::avcp
436  ! REAL(KIND(1D0))::avdens
437  ! REAL(KIND(1D0))::lv_J_kg
438  REAL(KIND(1D0))::dq
439  REAL(KIND(1D0))::lvS_J_kg
440  REAL(KIND(1D0))::psyc_hPa
441  REAL(KIND(1D0))::RAsnow
442  REAL(KIND(1D0))::rb
443  REAL(KIND(1D0))::runoff_per_interval
444  REAL(KIND(1D0))::s_hPa
445  REAL(KIND(1D0))::sIce_hpa
446  REAL(KIND(1D0))::SoilMoistCap
447  REAL(KIND(1D0))::veg_fr
448  REAL(KIND(1D0))::VegPhenLumps
449  REAL(KIND(1D0))::VPd_hpa
450  REAL(KIND(1D0))::vsmd
451  REAL(KIND(1D0))::ZZD
452 
453  REAL(KIND(1D0)), DIMENSION(NSURF)::deltaQi
454  REAL(KIND(1D0)), DIMENSION(NSURF)::drain
455  REAL(KIND(1D0)), DIMENSION(NSURF)::FreezState
456  REAL(KIND(1D0)), DIMENSION(NSURF)::FreezStateVol
457  REAL(KIND(1D0)), DIMENSION(NSURF)::soilstoreOld
458  REAL(KIND(1D0)), DIMENSION(NSURF)::stateOld
459  REAL(KIND(1D0)), DIMENSION(NSURF)::tsurf_ind
460 
461  ! TODO: TS 25 Oct 2017
462  ! the `add-*` variables are not used currently as grid-to-grid connection is NOT set up.
463  ! set these variables as zero.
464  REAL(KIND(1D0))::addImpervious = 0
465  REAL(KIND(1D0))::addPipes = 0
466  REAL(KIND(1D0))::addVeg = 0
467  REAL(KIND(1D0))::addWaterBody = 0
468  REAL(KIND(1D0)), DIMENSION(NSURF)::AddWater = 0
469  REAL(KIND(1D0)), DIMENSION(NSURF)::AddWaterRunoff = 0
470 
471  ! values that are derived from tstep
472  INTEGER::nsh ! number of timesteps per hour
473  REAL(KIND(1D0))::nsh_real ! nsh in type real
474  REAL(KIND(1D0))::tstep_real ! tstep in type real
475  REAL(KIND(1D0))::dectime
476 
477  ! values that are derived from sfr (surface fractions)
478  REAL(KIND(1D0))::VegFraction
479  REAL(KIND(1D0))::ImpervFraction
480  REAL(KIND(1D0))::PervFraction
481  REAL(KIND(1D0))::NonWaterFraction
482 
483  ! snow related temporary values
484  ! REAL(KIND(1D0))::alb0
485  REAL(KIND(1D0))::alb1
486 
487  ! ########################################################################################
488  ! TS 19 Sep 2019
489  ! temporary variables to save values for inout varialbes
490  ! suffixes `prev` and `next` denote values from last and to next tsteps, respectively
491  ! these variables are introduced to allow safe and robust iterations inccurred in this subroutine
492  ! so that these values won't updated in unexpectedly many times
493 
494  ! OHM related:
495  REAL(KIND(1D0))::qn1_av_prev, qn1_av_next
496  REAL(KIND(1D0))::dqndt_prev, dqndt_next
497  REAL(KIND(1D0))::qn1_s_av_prev, qn1_s_av_next
498  REAL(KIND(1D0))::dqnsdt_prev, dqnsdt_next
499 
500  ! snow related:
501  REAL(KIND(1D0))::SnowfallCum_prev, SnowfallCum_next
502  REAL(KIND(1D0))::SnowAlb_prev, SnowAlb_next
503 
504  REAL(KIND(1D0)), DIMENSION(NSURF)::IceFrac_prev, IceFrac_next
505  REAL(KIND(1D0)), DIMENSION(NSURF)::SnowWater_prev, SnowWater_next
506  REAL(KIND(1D0)), DIMENSION(NSURF)::SnowDens_prev, SnowDens_next
507  REAL(KIND(1D0)), DIMENSION(NSURF)::SnowFrac_prev, SnowFrac_next
508  REAL(KIND(1D0)), DIMENSION(NSURF)::SnowPack_prev, SnowPack_next
509 
510  ! water balance related:
511  REAL(KIND(1D0)), DIMENSION(NSURF) ::soilstore_id_prev, soilstore_id_next
512  REAL(KIND(1D0)), DIMENSION(NSURF) ::state_id_prev, state_id_next
513  REAL(KIND(1D0)), DIMENSION(6, NSURF)::StoreDrainPrm_prev, StoreDrainPrm_next
514 
515  ! phenology related:
516  REAL(KIND(1D0)), DIMENSION(NSURF) ::alb_prev, alb_next
517  REAL(KIND(1d0)), DIMENSION(nvegsurf)::GDD_id_prev, GDD_id_next
518  REAL(KIND(1d0)), DIMENSION(nvegsurf)::LAI_id_prev, LAI_id_next
519  REAL(KIND(1d0)), DIMENSION(nvegsurf)::SDD_id_prev, SDD_id_next
520 
521  REAL(KIND(1D0))::DecidCap_id_prev, DecidCap_id_next
522  REAL(KIND(1D0))::albDecTr_id_prev, albDecTr_id_next
523  REAL(KIND(1D0))::albEveTr_id_prev, albEveTr_id_next
524  REAL(KIND(1D0))::albGrass_id_prev, albGrass_id_next
525  REAL(KIND(1D0))::porosity_id_prev, porosity_id_next
526 
527  REAL(KIND(1d0))::Tmin_id_prev, Tmin_id_next
528  REAL(KIND(1d0))::Tmax_id_prev, Tmax_id_next
529  REAL(KIND(1d0))::lenDay_id_prev, lenDay_id_next
530 
531  ! anthropogenic heat related:
532  REAL(KIND(1d0)), DIMENSION(12)::HDD_id_prev, HDD_id_next
533 
534  ! water use related:
535  REAL(KIND(1d0)), DIMENSION(9)::WUDay_id_prev, WUDay_id_next
536 
537  REAL(KIND(1D0))::Tair_av_prev, Tair_av_next
538  ! ########################################################################################
539 
540  ! Related to RSL wind profiles
541  INTEGER, PARAMETER :: nz = 90 ! number of levels 10 levels in canopy plus 20 (3 x Zh) above the canopy
542 
543  ! flag for Tsurf convergence
544  logical:: flag_converge
545  REAL(KIND(1D0)):: Ts_iter
546  REAL(KIND(1D0)):: L_mod_iter
547  REAL(KIND(1D0)):: QH_Init
548  INTEGER::i_iter
549 
550  REAL(KIND(1d0)), DIMENSION(30):: psihatm_z
551  REAL(KIND(1d0)), DIMENSION(30):: psihath_z
552 
553  ! ########################################################################################
554  ! save initial values of inout variables
555  qn1_av_prev = qn1_av
556  dqndt_prev = dqndt
557  qn1_s_av_prev = qn1_s_av
558  dqnsdt_prev = dqnsdt
559  snowfallcum_prev = snowfallcum
560  snowalb_prev = snowalb
561  icefrac_prev = icefrac
562  snowwater_prev = snowwater
563  snowdens_prev = snowdens
564  snowfrac_prev = snowfrac
565  snowpack_prev = snowpack
566  soilstore_id_prev = soilstore_id
567  state_id_prev = state_id
568  tair_av_prev = tair_av
569  lai_id_prev = lai_id
570  gdd_id_prev = gdd_id
571  sdd_id_prev = sdd_id
572  tmin_id_prev = tmin_id
573  tmax_id_prev = tmax_id
574  lenday_id_prev = lenday_id
575  storedrainprm_prev = storedrainprm
576  decidcap_id_prev = decidcap_id
577  porosity_id_prev = porosity_id
578  alb_prev = alb
579  albdectr_id_prev = albdectr_id
580  albevetr_id_prev = albevetr_id
581  albgrass_id_prev = albgrass_id
582  hdd_id_prev = hdd_id
583  wuday_id_prev = wuday_id
584 
585  ! initialise `_next` variables
586  qn1_av_next = qn1_av
587  dqndt_next = dqndt
588  qn1_s_av_next = qn1_s_av
589  dqnsdt_next = dqnsdt
590  snowfallcum_next = snowfallcum
591  snowalb_next = snowalb
592  icefrac_next = icefrac
593  snowwater_next = snowwater
594  snowdens_next = snowdens
595  snowfrac_next = snowfrac
596  snowpack_next = snowpack
597  soilstore_id_next = soilstore_id
598  state_id_next = state_id
599  tair_av_next = tair_av
600  lai_id_next = lai_id
601  gdd_id_next = gdd_id
602  sdd_id_next = sdd_id
603  tmin_id_next = tmin_id
604  tmax_id_next = tmax_id
605  lenday_id_next = lenday_id
606  storedrainprm_next = storedrainprm
607  decidcap_id_next = decidcap_id
608  porosity_id_next = porosity_id
609  alb_next = alb
610  albdectr_id_next = albdectr_id
611  albevetr_id_next = albevetr_id
612  albgrass_id_next = albgrass_id
613  hdd_id_next = hdd_id
614  wuday_id_next = wuday_id
615 
616  !########################################################################################
617  ! main calculation starts here
618  !########################################################################################
619 
620  ! iteration is used below to get results converge
621  flag_converge = .false.
622  ts_iter = temp_c
623  l_mod_iter = 10
624  i_iter = 1
625  do while (.not. flag_converge)
626 
627  ! calculate dectime
628  CALL suews_cal_dectime( &
629  id, it, imin, isec, & ! input
630  dectime) ! output
631 
632  ! calculate tstep related VARIABLES
633  CALL suews_cal_tstep( &
634  tstep, & ! input
635  nsh, nsh_real, tstep_real) ! output
636 
637  ! calculate surface fraction related VARIABLES
638  CALL suews_cal_surf( &
639  sfr, & !input
640  vegfraction, impervfraction, pervfraction, nonwaterfraction) ! output
641 
642  ! calculate dayofweek information
643  CALL suews_cal_weekday( &
644  iy, id, lat, & !input
645  dayofweek_id) !output
646 
647  ! calculate dayofweek information
648  CALL suews_cal_dls( &
649  id, startdls, enddls, & !input
650  dls) !output
651 
652  ! calculate mean air temperature of past 24 hours
653  tair_av_next = cal_tair_av(tair_av_prev, dt_since_start, tstep, temp_c)
654 
655  !==============main calculation start=======================
656 
657  !==============surface roughness calculation=======================
658  IF (diagnose == 1) WRITE (*, *) 'Calling SUEWS_cal_RoughnessParameters...'
659  IF (diagnose == 1) print *, 'z0m_in =', z0m_in
661  roughlenmommethod, sfr, &!input
662  bldgh, evetreeh, dectreeh, &
663  porosity_id_prev, faibldg, faievetree, faidectree, &
664  z0m_in, zdm_in, z, &
665  planf, &!output
666  zh, z0m, zdm, zzd)
667 
668  !=================Calculate sun position=================
669  IF (diagnose == 1) WRITE (*, *) 'Calling NARP_cal_SunPosition...'
670  CALL narp_cal_sunposition( &
671  REAL(iy, KIND(1d0)), &!input:
672  dectime - tstep/2, &! sun position at middle of timestep before
673  timezone, lat, lng, alt, &
674  azimuth, zenith_deg)!output:
675 
676  !=================Call the SUEWS_cal_DailyState routine to get surface characteristics ready=================
677  IF (diagnose == 1) WRITE (*, *) 'Calling SUEWS_cal_DailyState...'
678  CALL suews_cal_dailystate( &
679  iy, id, it, imin, isec, tstep, tstep_prev, dt_since_start, dayofweek_id, &!input
680  tmin_id_prev, tmax_id_prev, lenday_id_prev, &
681  waterusemethod, ie_start, ie_end, &
682  laicalcyes, laitype, &
683  nsh_real, avkdn, temp_c, precip, basethdd, &
684  lat, faut, lai_obs, &
685  albmax_dectr, albmax_evetr, albmax_grass, &
686  albmin_dectr, albmin_evetr, albmin_grass, &
687  capmax_dec, capmin_dec, pormax_dec, pormin_dec, &
688  ie_a, ie_m, daywatper, daywat, &
689  baset, basete, gddfull, sddfull, laimin, laimax, laipower, &
690  decidcap_id_prev, storedrainprm_prev, lai_id_prev, gdd_id_prev, sdd_id_prev, &
691  albdectr_id_prev, albevetr_id_prev, albgrass_id_prev, porosity_id_prev, &!input
692  hdd_id_prev, &!input
693  hdd_id_next, &!output
694  tmin_id_next, tmax_id_next, lenday_id_next, &
695  albdectr_id_next, albevetr_id_next, albgrass_id_next, porosity_id_next, &!output
696  decidcap_id_next, storedrainprm_next, lai_id_next, gdd_id_next, sdd_id_next, deltalai, wuday_id_next)!output
697 
698  !=================Calculation of density and other water related parameters=================
699  IF (diagnose == 1) WRITE (*, *) 'Calling LUMPS_cal_AtmMoist...'
700  CALL cal_atmmoist( &
701  temp_c, press_hpa, avrh, dectime, &! input:
702  lv_j_kg, lvs_j_kg, &! output:
703  es_hpa, ea_hpa, vpd_hpa, vpd_pa, dq, dens_dry, avcp, avdens)
704 
705  !======== Calculate soil moisture =========
706  IF (diagnose == 1) WRITE (*, *) 'Calling SUEWS_update_SoilMoist...'
707  CALL suews_update_soilmoist( &
708  nonwaterfraction, &!input
709  soilstorecap, sfr, soilstore_id_prev, &
710  soilmoistcap, soilstate, &!output
711  vsmd, smd)
712 
713  IF (diagnose == 1) WRITE (*, *) 'Calling SUEWS_cal_WaterUse...'
714  !=================Gives the external and internal water uses per timestep=================
715  CALL suews_cal_wateruse( &
716  nsh_real, & ! input:
717  wu_m3, surfacearea, sfr, &
718  irrfracconif, irrfracdecid, irrfracgrass, &
719  dayofweek_id, wuprofa_24hr, wuprofm_24hr, &
720  internalwateruse_h, hdd_id_next, wuday_id_next, &
721  waterusemethod, nsh, it, imin, dls, &
722  wu_evetr, wu_dectr, wu_grass, int_wu, ext_wu)! output:
723 
724  ! ===================ANTHROPOGENIC HEAT AND CO2 FLUX======================
726  ah_min, ahprof_24hr, ah_slope_cooling, ah_slope_heating, co2pointsource, &! input:
727  dayofweek_id, diagnose, dls, ef_umolco2perj, emissionsmethod, enef_v_jkm, &
728  fcef_v_kgkm, frfossilfuel_heat, frfossilfuel_nonheat, hdd_id_next, humactivity_24hr, &
729  id, imin, it, maxfcmetab, maxqfmetab, minfcmetab, minqfmetab, nsh, &
730  popdensdaytime, popdensnighttime, popprof_24hr, qf, qf0_beu, qf_a, qf_b, qf_c, &
731  qf_obs, qf_sahp, sfr, snowfrac, surfacearea, t_critic_cooling, t_critic_heating, &
732  temp_c, trafficrate, trafficunits, traffprof_24hr, &
733  fc_anthro, fc_build, fc_metab, fc_point, fc_traff)! output:
734 
735  ! ========================================================================
736  ! N.B.: the following parts involves snow-related calculations.
737 
738  ! ===================NET ALLWAVE RADIATION================================
739  CALL suews_cal_qn( &
740  netradiationmethod, snowuse, &!input
741  tstep, snowpack_prev, tau_a, tau_f, snowalbmax, snowalbmin, &
742  diagnose, snowfrac_obs, ldown_obs, fcld_obs, &
743  dectime, zenith_deg, ts_iter, avkdn, temp_c, avrh, ea_hpa, qn1_obs, &
744  snowalb_prev, snowfrac_prev, diagqn, &
745  narp_trans_site, narp_emis_snow, icefrac_prev, sfr, emis, &
746  alb_prev, albdectr_id_next, albevetr_id_next, albgrass_id_next, &!input
747  alb_next, &!output
748  ldown, fcld, &!output
749  qn1, qn1_snowfree, qn1_s, kclear, kup, lup, tsurf, &
750  qn1_ind_snow, kup_ind_snow, tsurf_ind_snow, tsurf_ind, &
751  alb1, snowfrac_next, snowalb_next)
752 
753  ! =================STORAGE HEAT FLUX=======================================
754  CALL suews_cal_qs( &
755  storageheatmethod, qs_obs, ohmincqf, gridiv, &!input
756  id, tstep, dt_since_start, diagnose, sfr, &
757  ohm_coef, ohm_threshsw, ohm_threshwd, &
758  soilstore_id_prev, soilstorecap, state_id_prev, snowuse, snowfrac_next, diagqs, &
759  hdd_id_next, metforcingdata_grid, ts5mindata_ir, qf, qn1, &
760  avkdn, avu1, temp_c, zenith_deg, avrh, press_hpa, ldown, &
761  bldgh, alb_next, emis, cpanohm, kkanohm, chanohm, emissionsmethod, &
762  tair_av_next, qn1_av_prev, dqndt_prev, qn1_s_av_prev, dqnsdt_prev, &
763  storedrainprm_next, &
764  qn1_s, dataoutlineestm, qs, &!output
765  qn1_av_next, dqndt_next, qn1_s_av_next, dqnsdt_next, &
766  deltaqi, a1, a2, a3)
767 
768  !==================Energy related to snow melting/freezing processes=======
769  IF (diagnose == 1) WRITE (*, *) 'Calling MeltHeat'
770 
771  CALL snow_cal_meltheat( &
772  snowuse, &!input
773  tstep, tau_r, snowdensmax, &
774  lvs_j_kg, lv_j_kg, tstep_real, radmeltfact, tempmeltfact, snowalbmax, &
775  snowdensmin, temp_c, precip, preciplimit, preciplimitalb, &
776  nsh_real, sfr, tsurf_ind, tsurf_ind_snow, state_id_prev, qn1_ind_snow, &
777  kup_ind_snow, snowwater_prev, deltaqi, alb1, &
778  snowpack_prev, snowfrac_next, snowalb_next, snowdens_prev, snowfallcum_prev, &!input
779  snowpack_next, snowfrac_next, snowalb_next, snowdens_next, snowfallcum_next, &!output
780  mwh, qm, qmfreez, qmrain, &! output
781  veg_fr, snowcalcswitch, qm_melt, qm_freezstate, qm_rain, freezmelt, &
782  freezstate, freezstatevol, rainonsnow, snowdepth, mw_ind, &
783  dataoutlinesnow)!output
784 
785  !==========================Turbulent Fluxes================================
786  IF (diagnose == 1) WRITE (*, *) 'Calling LUMPS_cal_QHQE...'
787  if (i_iter == 1) then
788  !Calculate QH and QE from LUMPS in the first iteration of each time step
789  CALL lumps_cal_qhqe( &
790  veg_type, & !input
791  snowuse, qn1, qf, qs, qm, temp_c, veg_fr, avcp, press_hpa, lv_j_kg, &
792  tstep_real, drainrt, nsh_real, &
793  precip, rainmaxres, raincover, sfr, lai_id_next, laimax, laimin, &
794  qh_lumps, & !output
795  qe_lumps, psyc_hpa, s_hpa, sice_hpa, tempveg, vegphenlumps)
796 
797  ! use LUMPS QH to do stability correction
798  qh_init = qh_lumps
799  else
800  ! use SUEWS QH to do stability correction
801  qh_init = qh
802  end if
803 
804  !============= calculate water balance =============
805  CALL suews_cal_water( &
806  diagnose, &!input
807  snowuse, nonwaterfraction, addpipes, addimpervious, addveg, addwaterbody, &
808  state_id_prev, soilstore_id_prev, sfr, storedrainprm_next, waterdist, nsh_real, &
809  drain_per_tstep, & !output
810  drain, addwaterrunoff, &
811  additionalwater, runoffpipes, runoff_per_interval, &
812  addwater, stateold, soilstoreold)
813  !============= calculate water balance end =============
814 
815  !===============Resistance Calculations=======================
816  CALL suews_cal_resistance( &
817  stabilitymethod, &!input:
818  diagnose, aerodynamicresistancemethod, roughlenheatmethod, snowuse, &
819  id, it, gsmodel, smdmethod, &
820  avdens, avcp, qh_init, zzd, z0m, zdm, &
821  avu1, temp_c, vegfraction, avkdn, &
822  kmax, &
823  g1, g2, g3, g4, &
824  g5, g6, s1, s2, &
825  th, tl, &
826  dq, xsmd, vsmd, maxconductance, laimax, lai_id_next, snowfrac_next, sfr, &
827  ustar, tstar, l_mod, &!output
828  zl, gsc, resistsurf, ra, rasnow, rb)
829 
830  !======== Evaporation and surface state_id ========
831  CALL suews_cal_qe( &
832  diagnose, snowuse, &!input
833  tstep, imin, it, evapmethod, snowcalcswitch, dayofweek_id, crwmin, crwmax, &
834  dectime, lvs_j_kg, avrh, press_hpa, temp_c, &
835  rasnow, psyc_hpa, sice_hpa, &
836  pervfraction, vegfraction, addimpervious, qn1_snowfree, qf, qs, vpd_hpa, s_hpa, &
837  resistsurf, ra, rb, snowdensmin, precip, pipecapacity, runofftowater, &
838  nonwaterfraction, wu_evetr, wu_dectr, wu_grass, addveg, addwaterbody, snowlimpaved, snowlimbldg, &
839  surfacearea, flowchange, drain, wetthresh, stateold, mw_ind, soilstorecap, rainonsnow, &
840  freezmelt, freezstate, freezstatevol, qm_melt, qm_rain, tsurf_ind, sfr, &
841  statelimit, addwater, addwaterrunoff, storedrainprm_next, snowpacklimit, snowprof_24hr, &
842  snowpack_next, snowfrac_next, snowwater_prev, icefrac_prev, snowdens_next, &! input:
843  runoff_per_interval, state_id_prev, soilstore_id_prev, &! input:
844  state_id_next, soilstore_id_next, &! output:
845  snowpack_next, snowfrac_next, snowwater_next, icefrac_next, snowdens_next, &! output
846  runoffsoil, &! output:
847  snowremoval, &
848  state_per_tstep, nwstate_per_tstep, qe, &
849  swe, chsnow_per_interval, ev_per_tstep, runoff_per_tstep, &
850  surf_chang_per_tstep, runoffpipes, mwstore, runoffwaterbody, &
851  runoffagveg, runoffagimpervious)
852  !======== Evaporation and surface state_id end========
853 
854  !============ Sensible heat flux ===============
855  IF (diagnose == 1) WRITE (*, *) 'Calling SUEWS_cal_QH...'
856  CALL suews_cal_qh( &
857  1, &
858  qn1, qf, qmrain, qe, qs, qmfreez, qm, avdens, avcp, tsurf, temp_c, ra, &
859  qh, qh_residual, qh_resist)!output
860  !============ Sensible heat flux end===============
861 
862  ! N.B.: snow-related calculations end here.
863  !===================================================
864 
865  !=== Horizontal movement between soil stores ===
866  ! Now water is allowed to move horizontally between the soil stores
867  IF (diagnose == 1) WRITE (*, *) 'Calling SUEWS_cal_HorizontalSoilWater...'
869  sfr, &! input: ! surface fractions
870  soilstorecap, &!Capacity of soil store for each surface [mm]
871  soildepth, &!Depth of sub-surface soil store for each surface [mm]
872  sathydraulicconduct, &!Saturated hydraulic conductivity for each soil subsurface [mm s-1]
873  surfacearea, &!Surface area of the study area [m2]
874  nonwaterfraction, &! sum of surface cover fractions for all except water surfaces
875  tstep_real, & !tstep cast as a real for use in calculations
876  soilstore_id_next, &! inout:!Soil moisture of each surface type [mm]
877  runoffsoil, &!Soil runoff from each soil sub-surface [mm]
878  runoffsoil_per_tstep &! output:!Runoff to deep soil per timestep [mm] (for whole surface, excluding water body)
879  )
880 
881  !========== Calculate soil moisture ============
882  IF (diagnose == 1) WRITE (*, *) 'Calling SUEWS_cal_SoilState...'
883  CALL suews_cal_soilstate( &
884  smdmethod, xsmd, nonwaterfraction, soilmoistcap, &!input
885  soilstorecap, surf_chang_per_tstep, &
886  soilstore_id_next, soilstoreold, sfr, &
887  smd, smd_nsurf, tot_chang_per_tstep, soilstate)!output
888 
889  !============ surface-level diagonostics ===============
890  ! IF (Diagnose == 1) WRITE (*, *) 'Calling SUEWS_cal_Diagnostics...'
891  ! CALL SUEWS_cal_Diagnostics( &
892  ! dectime, &!input
893  ! avU1, Temp_C, avRH, Press_hPa, &
894  ! qh, qe, &
895  ! VegFraction, z, z0m, zdm, RA, avdens, avcp, lv_J_kg, tstep_real, &
896  ! RoughLenHeatMethod, StabilityMethod, &
897  ! U10_ms, t2_C, q2_gkg, tsfc_C, RH2)!output
898 
899  !============ roughness sub-layer diagonostics ===============
900  ! IF (Diagnose == 1) WRITE (*, *) 'Calling RSLProfile...'
901  ! CALL RSLProfile( &
902  ! Zh, z0m, zdm, &
903  ! UStar, L_MOD, sfr, planF, StabilityMethod, &
904  ! avcp, lv_J_kg, &
905  ! Temp_C, avRH, Press_hPa, z, qh, qe, & ! input
906  ! T2_C, q2_gkg, U10_ms, RH2, & !output
907  ! dataoutLineRSL) ! output
908 
909  !============ calculate surface temperature ===============
910  tsfc_c = cal_tsfc(qh, avdens, avcp, ra, temp_c)
911 
912  !============ surface-level diagonostics end ===============
913 
914  ! ============ BIOGENIC CO2 FLUX =======================
915  ! CALL SUEWS_cal_BiogenCO2( &
916  ! alpha_bioCO2, alpha_enh_bioCO2, avkdn, avRh, beta_bioCO2, beta_enh_bioCO2, BSoilSurf, &! input:
917  ! ConifSurf, DecidSurf, dectime, Diagnose, EmissionsMethod, Fc_anthro, G1, G2, G3, G4, &
918  ! G5, G6, gfunc, GrassSurf, gsmodel, id, it, ivConif, ivDecid, ivGrass, Kmax, LAI_id_next, LAIMin, &
919  ! LAIMax, MaxConductance, min_res_bioCO2, nsurf, NVegSurf, Press_hPa, resp_a, &
920  ! resp_b, S1, S2, sfr, SMDMethod, SnowFrac, t2_C, Temp_C, theta_bioCO2, TH, TL, vsmd, xsmd, &
921  ! Fc, Fc_biogen, Fc_photo, Fc_respi)! output:
922 
923  ! force quit do-while, i.e., skip iteration and use NARP for Tsurf calculation
924  ! if (NetRadiationMethod < 10 .or. NetRadiationMethod > 100) exit
925 
926  ! Test if sensible heat fluxes converge in iterations
927  i_iter = i_iter + 1
928  if (abs(qh - qh_init) > 0.1) then
929  flag_converge = .false.
930  else
931  flag_converge = .true.
932  end if
933 
934  ! force quit do-while loop if not convergent after 100 iterations
935  if (i_iter > 100) exit
936 
937  ts_iter = tsfc_c
938  l_mod_iter = l_mod
939 
940  !==============main calculation end=======================
941  ENDDO ! end iteration for tsurf calculations
942 
943  !==============================================================
944  ! Calculate diagnostics: these variables are decoupled from the main SUEWS calculation
945 
946  !============ roughness sub-layer diagonostics ===============
947  IF (diagnose == 1) WRITE (*, *) 'Calling RSLProfile...'
948  CALL rslprofile( &
949  zh, z0m, zdm, &
950  l_mod, sfr, planf, stabilitymethod, &
951  avcp, lv_j_kg, avdens, &
952  avu1, temp_c, avrh, press_hpa, z, qh, qe, & ! input
953  t2_c, q2_gkg, u10_ms, rh2, & !output
954  dataoutlinersl) ! output
955 
956  ! ============ BIOGENIC CO2 FLUX =======================
957  CALL suews_cal_biogenco2( &
958  alpha_bioco2, alpha_enh_bioco2, avkdn, avrh, beta_bioco2, beta_enh_bioco2, bsoilsurf, &! input:
959  conifsurf, decidsurf, dectime, diagnose, emissionsmethod, fc_anthro, g1, g2, g3, g4, &
960  g5, g6, gfunc, grasssurf, gsmodel, id, it, ivconif, ivdecid, ivgrass, kmax, lai_id_next, laimin, &
961  laimax, maxconductance, min_res_bioco2, nsurf, nvegsurf, press_hpa, resp_a, &
962  resp_b, s1, s2, sfr, smdmethod, snowfrac, t2_c, temp_c, theta_bioco2, th, tl, vsmd, xsmd, &
963  fc, fc_biogen, fc_photo, fc_respi)! output:
964 
965  ! calculations of diagnostics end
966  !==============================================================
967 
968  !==============================================================
969  ! update inout variables with new values
970  qn1_av = qn1_av_next
971  dqndt = dqndt_next
972  qn1_s_av = qn1_s_av_next
973  dqnsdt = dqnsdt_next
974  snowfallcum = snowfallcum_next
975  snowalb = snowalb_next
976  icefrac = icefrac_next
977  snowwater = snowwater_next
978  snowdens = snowdens_next
979  snowfrac = snowfrac_next
980  snowpack = snowpack_next
981 
982  soilstore_id = soilstore_id_next
983  state_id = state_id_next
984  alb = alb_next
985  gdd_id = gdd_id_next
986  sdd_id = sdd_id_next
987  lai_id = lai_id_next
988  decidcap_id = decidcap_id_next
989  albdectr_id = albdectr_id_next
990  albevetr_id = albevetr_id_next
991  albgrass_id = albgrass_id_next
992  porosity_id = porosity_id_next
993  storedrainprm = storedrainprm_next
994  tair_av = tair_av_next
995  tmin_id = tmin_id_next
996  tmax_id = tmax_id_next
997  lenday_id = lenday_id_next
998  hdd_id = hdd_id_next
999  wuday_id = wuday_id_next
1000 
1001  !==============translation of output variables into output array===========
1002  CALL suews_update_outputline( &
1003  additionalwater, alb, avkdn, u10_ms, azimuth, &!input
1004  chsnow_per_interval, dectime, &
1005  drain_per_tstep, qe_lumps, ev_per_tstep, ext_wu, fc, fc_build, fcld, &
1006  fc_metab, fc_photo, fc_respi, fc_point, fc_traff, flowchange, &
1007  qh_lumps, id, imin, int_wu, it, iy, &
1008  kup, lai_id, ldown, l_mod, lup, mwh, &
1009  mwstore, &
1010  nsh_real, nwstate_per_tstep, precip, q2_gkg, &
1011  qe, qf, qh, qh_resist, qm, qmfreez, &
1012  qmrain, qn1, qn1_s, qn1_snowfree, qs, ra, &
1013  resistsurf, rh2, runoffagimpervious, runoffagveg, &
1014  runoff_per_tstep, runoffpipes, runoffsoil_per_tstep, &
1015  runoffwaterbody, sfr, smd, smd_nsurf, snowalb, snowremoval, &
1016  state_id_next, state_per_tstep, surf_chang_per_tstep, swe, t2_c, tsfc_c, &
1017  tot_chang_per_tstep, tsurf, ustar, wu_dectr, &
1018  wu_evetr, wu_grass, z0m, zdm, zenith_deg, &
1019  datetimeline, dataoutlinesuews)!output
1020 
1021  ! model state_id:
1022 
1023  ! daily state_id:
1024  CALL update_dailystateline( &
1025  it, imin, nsh_real, &!input
1026  gdd_id, hdd_id, lai_id, &
1027  sdd_id, &
1028  tmin_id, tmax_id, lenday_id, &
1029  decidcap_id, &
1030  albdectr_id, &
1031  albevetr_id, &
1032  albgrass_id, &
1033  porosity_id, &
1034  wuday_id, &
1035  deltalai, vegphenlumps, &
1036  snowalb, snowdens, &
1037  a1, a2, a3, &
1038  dailystateline)!out
1039 
1040  !==============translation end ================
1041 
1042  END SUBROUTINE suews_cal_main
1043  ! ================================================================================
1044 
1045  ! ===================ANTHROPOGENIC HEAT + CO2 FLUX================================
1046  SUBROUTINE suews_cal_anthropogenicemission( &
1047  AH_MIN, AHProf_24hr, AH_SLOPE_Cooling, AH_SLOPE_Heating, CO2PointSource, &! input:
1048  dayofWeek_id, Diagnose, DLS, EF_umolCO2perJ, EmissionsMethod, EnEF_v_Jkm, &
1049  FcEF_v_kgkm, FrFossilFuel_Heat, FrFossilFuel_NonHeat, HDD_id, HumActivity_24hr, &
1050  id, imin, it, MaxFCMetab, MaxQFMetab, MinFCMetab, MinQFMetab, nsh, &
1051  PopDensDaytime, PopDensNighttime, PopProf_24hr, QF, QF0_BEU, Qf_A, Qf_B, Qf_C, &
1052  QF_obs, QF_SAHP, sfr, SnowFrac, SurfaceArea, T_CRITIC_Cooling, T_CRITIC_Heating, &
1053  Temp_C, TrafficRate, TrafficUnits, TraffProf_24hr, &
1054  Fc_anthro, Fc_build, Fc_metab, Fc_point, Fc_traff)! output:
1056  IMPLICIT NONE
1057 
1058  INTEGER, INTENT(in)::Diagnose
1059  INTEGER, INTENT(in)::DLS
1060  INTEGER, INTENT(in)::EmissionsMethod
1061  INTEGER, INTENT(in)::id
1062  INTEGER, INTENT(in)::it
1063  INTEGER, INTENT(in)::imin
1064  INTEGER, INTENT(in)::nsh
1065  INTEGER, DIMENSION(3), INTENT(in)::dayofWeek_id
1066 
1067  REAL(KIND(1d0)), DIMENSION(6, 2), INTENT(in)::HDD_id
1068 
1069  REAL(KIND(1d0)), DIMENSION(2), INTENT(in)::AH_MIN
1070  REAL(KIND(1d0)), DIMENSION(2), INTENT(in)::AH_SLOPE_Heating
1071  REAL(KIND(1d0)), DIMENSION(2), INTENT(in)::AH_SLOPE_Cooling
1072  REAL(KIND(1D0)), DIMENSION(2), INTENT(in)::FcEF_v_kgkm
1073  ! REAL(KIND(1d0)), DIMENSION(2), INTENT(in)::NumCapita
1074  REAL(KIND(1d0)), DIMENSION(2), INTENT(in)::PopDensDaytime
1075  REAL(KIND(1d0)), DIMENSION(2), INTENT(in)::QF0_BEU
1076  REAL(KIND(1d0)), DIMENSION(2), INTENT(in)::Qf_A
1077  REAL(KIND(1d0)), DIMENSION(2), INTENT(in)::Qf_B
1078  REAL(KIND(1d0)), DIMENSION(2), INTENT(in)::Qf_C
1079  REAL(KIND(1d0)), DIMENSION(2), INTENT(in)::T_CRITIC_Heating
1080  REAL(KIND(1d0)), DIMENSION(2), INTENT(in)::T_CRITIC_Cooling
1081  REAL(KIND(1d0)), DIMENSION(2), INTENT(in)::TrafficRate
1082 
1083  REAL(KIND(1d0)), DIMENSION(0:23, 2), INTENT(in)::AHProf_24hr
1084  REAL(KIND(1d0)), DIMENSION(0:23, 2), INTENT(in)::HumActivity_24hr
1085  REAL(KIND(1d0)), DIMENSION(0:23, 2), INTENT(in)::TraffProf_24hr
1086  REAL(KIND(1d0)), DIMENSION(0:23, 2), INTENT(in)::PopProf_24hr
1087 
1088  REAL(KIND(1D0)), INTENT(in)::CO2PointSource
1089  REAL(KIND(1D0)), INTENT(in)::EF_umolCO2perJ
1090  REAL(KIND(1D0)), INTENT(in)::EnEF_v_Jkm
1091  REAL(KIND(1D0)), INTENT(in)::FrFossilFuel_Heat
1092  REAL(KIND(1D0)), INTENT(in)::FrFossilFuel_NonHeat
1093  REAL(KIND(1D0)), INTENT(in)::MaxFCMetab
1094  REAL(KIND(1D0)), INTENT(in)::MaxQFMetab
1095  REAL(KIND(1D0)), INTENT(in)::MinFCMetab
1096  REAL(KIND(1D0)), INTENT(in)::MinQFMetab
1097  REAL(KIND(1D0)), INTENT(in)::PopDensNighttime
1098  REAL(KIND(1D0)), INTENT(in)::QF_obs
1099  REAL(KIND(1D0)), INTENT(in)::Temp_C
1100  REAL(KIND(1D0)), INTENT(in)::TrafficUnits
1101 
1102  REAL(KIND(1d0)), DIMENSION(nsurf), INTENT(in)::sfr
1103  REAL(KIND(1d0)), DIMENSION(nsurf), INTENT(in)::SnowFrac
1104  REAL(KIND(1D0)), INTENT(IN)::SurfaceArea
1105 
1106  REAL(KIND(1D0)), INTENT(out)::Fc_anthro
1107  REAL(KIND(1D0)), INTENT(out)::Fc_build
1108  REAL(KIND(1D0)), INTENT(out)::Fc_metab
1109  REAL(KIND(1D0)), INTENT(out)::Fc_point
1110  REAL(KIND(1D0)), INTENT(out)::Fc_traff
1111  REAL(KIND(1D0)), INTENT(out)::QF
1112  REAL(KIND(1D0)), INTENT(out)::QF_SAHP
1113 
1114  INTEGER, PARAMETER :: notUsedI = -999
1115  REAL(KIND(1D0)), PARAMETER::notUsed = -999
1116 
1117  IF (emissionsmethod == 0) THEN ! use observed qf
1118  qf = qf_obs
1119  ELSEIF ((emissionsmethod > 0 .AND. emissionsmethod <= 6) .OR. emissionsmethod >= 11) THEN
1120  CALL anthropogenicemissions( &
1121  co2pointsource, emissionsmethod, &
1122  id, it, imin, dls, nsh, dayofweek_id, &
1123  ef_umolco2perj, fcef_v_kgkm, enef_v_jkm, trafficunits, &
1124  frfossilfuel_heat, frfossilfuel_nonheat, &
1125  minfcmetab, maxfcmetab, minqfmetab, maxqfmetab, &
1126  popdensdaytime, popdensnighttime, &
1127  temp_c, hdd_id, qf_a, qf_b, qf_c, &
1128  ah_min, ah_slope_heating, ah_slope_cooling, &
1129  t_critic_heating, t_critic_cooling, &
1130  trafficrate, &
1131  qf0_beu, qf_sahp, &
1132  fc_anthro, fc_metab, fc_traff, fc_build, fc_point, &
1133  ahprof_24hr, humactivity_24hr, traffprof_24hr, popprof_24hr, surfacearea)
1134 
1135  ELSE
1136  CALL errorhint(73, 'RunControl.nml:EmissionsMethod unusable', notused, notused, emissionsmethod)
1137  ENDIF
1138 
1139  IF (emissionsmethod >= 1) qf = qf_sahp
1140 
1141  IF (emissionsmethod >= 0 .AND. emissionsmethod <= 6) THEN
1142  fc_anthro = 0
1143  fc_metab = 0
1144  fc_traff = 0
1145  fc_build = 0
1146  fc_point = 0
1147  ENDIF
1148 
1149  END SUBROUTINE suews_cal_anthropogenicemission
1150  ! ================================================================================
1151 
1152  !==============BIOGENIC CO2 flux==================================================
1153  SUBROUTINE suews_cal_biogenco2( &
1154  alpha_bioCO2, alpha_enh_bioCO2, avkdn, avRh, beta_bioCO2, beta_enh_bioCO2, BSoilSurf, &! input:
1155  ConifSurf, DecidSurf, dectime, Diagnose, EmissionsMethod, Fc_anthro, G1, G2, G3, G4, &
1156  G5, G6, gfunc, GrassSurf, gsmodel, id, it, ivConif, ivDecid, ivGrass, Kmax, LAI_id, LAIMin, &
1157  LAIMax, MaxConductance, min_res_bioCO2, nsurf, NVegSurf, Press_hPa, resp_a, &
1158  resp_b, S1, S2, sfr, SMDMethod, SnowFrac, t2_C, Temp_C, theta_bioCO2, TH, TL, vsmd, xsmd, &
1159  Fc, Fc_biogen, Fc_photo, Fc_respi)! output:
1161  IMPLICIT NONE
1162 
1163  REAL(KIND(1d0)), DIMENSION(nvegsurf), INTENT(in)::alpha_bioCO2
1164  REAL(KIND(1d0)), DIMENSION(nvegsurf), INTENT(in)::alpha_enh_bioCO2
1165  REAL(KIND(1d0)), DIMENSION(nvegsurf), INTENT(in)::beta_bioCO2
1166  REAL(KIND(1d0)), DIMENSION(nvegsurf), INTENT(in)::beta_enh_bioCO2
1167  REAL(KIND(1d0)), DIMENSION(nvegsurf), INTENT(in)::LAI_id
1168  REAL(KIND(1d0)), DIMENSION(nvegsurf), INTENT(in)::LAIMin
1169  REAL(KIND(1d0)), DIMENSION(nvegsurf), INTENT(in)::LAIMax
1170  REAL(KIND(1d0)), DIMENSION(nvegsurf), INTENT(in)::min_res_bioCO2
1171  REAL(KIND(1d0)), DIMENSION(nvegsurf), INTENT(in)::resp_a
1172  REAL(KIND(1d0)), DIMENSION(nvegsurf), INTENT(in)::resp_b
1173  REAL(KIND(1d0)), DIMENSION(nvegsurf), INTENT(in)::theta_bioCO2
1174 
1175  REAL(KIND(1d0)), DIMENSION(nsurf), INTENT(in)::sfr
1176  REAL(KIND(1d0)), DIMENSION(nsurf), INTENT(in)::SnowFrac
1177 
1178  REAL(KIND(1d0)), DIMENSION(3), INTENT(in)::MaxConductance
1179 
1180  INTEGER, INTENT(in)::BSoilSurf
1181  INTEGER, INTENT(in)::ConifSurf
1182  INTEGER, INTENT(in)::DecidSurf
1183  INTEGER, INTENT(in)::Diagnose
1184  INTEGER, INTENT(in)::EmissionsMethod
1185  INTEGER, INTENT(in)::GrassSurf
1186  INTEGER, INTENT(in)::gsmodel
1187  INTEGER, INTENT(in)::id
1188  INTEGER, INTENT(in)::it
1189  INTEGER, INTENT(in)::ivConif
1190  INTEGER, INTENT(in)::ivDecid
1191  INTEGER, INTENT(in)::ivGrass
1192  INTEGER, INTENT(in)::nsurf
1193  INTEGER, INTENT(in)::NVegSurf
1194  INTEGER, INTENT(in)::SMDMethod
1195 
1196  REAL(KIND(1d0)), INTENT(in)::avkdn
1197  REAL(KIND(1d0)), INTENT(in)::avRh
1198  REAL(KIND(1d0)), INTENT(in)::dectime
1199  REAL(KIND(1d0)), INTENT(in)::Fc_anthro
1200  REAL(KIND(1d0)), INTENT(in)::G1
1201  REAL(KIND(1d0)), INTENT(in)::G2
1202  REAL(KIND(1d0)), INTENT(in)::G3
1203  REAL(KIND(1d0)), INTENT(in)::G4
1204  REAL(KIND(1d0)), INTENT(in)::G5
1205  REAL(KIND(1d0)), INTENT(in)::G6
1206  REAL(KIND(1d0)), INTENT(in)::gfunc
1207  REAL(KIND(1d0)), INTENT(in)::Kmax
1208  REAL(KIND(1d0)), INTENT(in)::Press_hPa
1209  REAL(KIND(1d0)), INTENT(in)::S1
1210  REAL(KIND(1d0)), INTENT(in)::S2
1211  REAL(KIND(1d0)), INTENT(in)::t2_C
1212  REAL(KIND(1d0)), INTENT(in)::Temp_C
1213  REAL(KIND(1d0)), INTENT(in)::TH
1214  REAL(KIND(1d0)), INTENT(in)::TL
1215  REAL(KIND(1d0)), INTENT(in)::vsmd
1216  REAL(KIND(1d0)), INTENT(in)::xsmd
1217 
1218  REAL(KIND(1d0)), INTENT(out)::Fc_biogen
1219  REAL(KIND(1d0)), INTENT(out)::Fc_photo
1220  REAL(KIND(1d0)), INTENT(out)::Fc_respi
1221  REAL(KIND(1d0)), INTENT(out)::Fc
1222 
1223  REAL(KIND(1d0))::gfunc2
1224  REAL(KIND(1d0))::dq
1225  REAL(KIND(1d0))::t2
1226  REAL(KIND(1d0))::dummy1
1227  REAL(KIND(1d0))::dummy2
1228  REAL(KIND(1d0))::dummy3
1229  REAL(KIND(1d0))::dummy4
1230  REAL(KIND(1d0))::dummy5
1231  REAL(KIND(1d0))::dummy6
1232  REAL(KIND(1d0))::dummy7
1233  REAL(KIND(1d0))::dummy8
1234  REAL(KIND(1d0))::dummy9
1235  REAL(KIND(1d0))::dummy10
1236  REAL(KIND(1d0))::dummy11
1237 
1238  IF (emissionsmethod >= 11) THEN
1239 
1240  IF (gsmodel == 3 .OR. gsmodel == 4) THEN ! With modelled 2 meter temperature
1241  ! Call LUMPS_cal_AtmMoist for dq and SurfaceResistance for gfunc with 2 meter temperature
1242  ! If modelled 2 meter temperature is too different from measured air temperature then
1243  ! use temp_c
1244  IF (abs(temp_c - t2_c) > 5) THEN
1245  t2 = temp_c
1246  ELSE
1247  t2 = t2_c
1248  ENDIF
1249 
1250  CALL cal_atmmoist( &
1251  t2, press_hpa, avrh, dectime, &! input:
1252  dummy1, dummy2, &! output:
1253  dummy3, dummy4, dummy5, dummy6, dq, dummy7, dummy8, dummy9)
1254 
1255  CALL surfaceresistance( &
1256  id, it, &! input:
1257  smdmethod, snowfrac, sfr, avkdn, t2, dq, xsmd, vsmd, maxconductance, &
1258  laimax, lai_id, gsmodel, kmax, &
1259  g1, g2, g3, g4, g5, g6, th, tl, s1, s2, &
1260  gfunc2, dummy10, dummy11)! output:
1261  ENDIF
1262 
1263  ! Calculate CO2 fluxes from biogenic components
1264  IF (diagnose == 1) WRITE (*, *) 'Calling CO2_biogen...'
1265  CALL co2_biogen( &
1266  alpha_bioco2, alpha_enh_bioco2, avkdn, beta_bioco2, beta_enh_bioco2, bsoilsurf, &! input:
1267  conifsurf, decidsurf, dectime, emissionsmethod, gfunc, gfunc2, grasssurf, gsmodel, &
1268  id, it, ivconif, ivdecid, ivgrass, lai_id, laimin, laimax, min_res_bioco2, nsurf, &
1269  nvegsurf, resp_a, resp_b, sfr, snowfrac, t2, temp_c, theta_bioco2, &
1270  fc_biogen, fc_photo, fc_respi)! output:
1271  ENDIF
1272 
1273  IF (emissionsmethod >= 0 .AND. emissionsmethod <= 6) THEN
1274  fc_biogen = 0
1275  fc_photo = 0
1276  fc_respi = 0
1277  ENDIF
1278 
1279  fc = fc_anthro + fc_biogen
1280 
1281  END SUBROUTINE suews_cal_biogenco2
1282  !========================================================================
1283 
1284  !=============net all-wave radiation=====================================
1285  SUBROUTINE suews_cal_qn( &
1286  NetRadiationMethod, snowUse, &!input
1287  tstep, SnowPack_prev, tau_a, tau_f, SnowAlbMax, SnowAlbMin, &
1288  Diagnose, snowFrac_obs, ldown_obs, fcld_obs, &
1289  dectime, ZENITH_deg, Tsurf_0, avKdn, Temp_C, avRH, ea_hPa, qn1_obs, &
1290  SnowAlb_prev, snowFrac_prev, DiagQN, &
1291  NARP_TRANS_SITE, NARP_EMIS_SNOW, IceFrac, sfr, emis, &
1292  alb_prev, albDecTr_id, albEveTr_id, albGrass_id, &!input
1293  alb_next, &!output
1294  ldown, fcld, &!output
1295  qn1, qn1_snowfree, qn1_S, kclear, kup, lup, tsurf, &
1296  qn1_ind_snow, kup_ind_snow, Tsurf_ind_snow, Tsurf_ind, &
1297  alb1, snowFrac_next, SnowAlb_next)
1299 
1300  IMPLICIT NONE
1301  ! INTEGER,PARAMETER ::nsurf = 7 ! number of surface types
1302  ! INTEGER,PARAMETER ::ConifSurf = 3 !New surface classes: Grass = 5th/7 surfaces
1303  ! INTEGER,PARAMETER ::DecidSurf = 4 !New surface classes: Grass = 5th/7 surfaces
1304  ! INTEGER,PARAMETER ::GrassSurf = 5
1305 
1306  INTEGER, INTENT(in)::NetRadiationMethod
1307  INTEGER, INTENT(in)::snowUse
1308  INTEGER, INTENT(in)::Diagnose
1309  INTEGER, INTENT(in)::DiagQN
1310  INTEGER, INTENT(in)::tstep
1311 
1312  REAL(KIND(1d0)), INTENT(in)::snowFrac_obs
1313  REAL(KIND(1d0)), INTENT(in)::ldown_obs
1314  REAL(KIND(1d0)), INTENT(in)::fcld_obs
1315  REAL(KIND(1d0)), INTENT(in)::dectime
1316  REAL(KIND(1d0)), INTENT(in)::ZENITH_deg
1317  REAL(KIND(1d0)), INTENT(in)::Tsurf_0
1318  REAL(KIND(1d0)), INTENT(in)::avKdn
1319  REAL(KIND(1d0)), INTENT(in)::Temp_C
1320  REAL(KIND(1d0)), INTENT(in)::avRH
1321  REAL(KIND(1d0)), INTENT(in)::ea_hPa
1322  REAL(KIND(1d0)), INTENT(in)::qn1_obs
1323  REAL(KIND(1d0)), INTENT(in)::SnowAlb_prev
1324  REAL(KIND(1d0)), INTENT(in)::NARP_EMIS_SNOW
1325  REAL(KIND(1d0)), INTENT(in)::NARP_TRANS_SITE
1326  REAL(KIND(1d0)), INTENT(in)::tau_a, tau_f, SnowAlbMax, SnowAlbMin
1327 
1328  REAL(KIND(1d0)), DIMENSION(nsurf), INTENT(in):: IceFrac
1329  REAL(KIND(1d0)), DIMENSION(nsurf), INTENT(in):: sfr
1330  REAL(KIND(1d0)), DIMENSION(nsurf), INTENT(in):: emis
1331 
1332  REAL(KIND(1d0)), DIMENSION(nsurf) ::alb
1333  REAL(KIND(1d0)), DIMENSION(nsurf), INTENT(in) ::alb_prev
1334  REAL(KIND(1d0)), DIMENSION(nsurf), INTENT(out) ::alb_next
1335  REAL(KIND(1d0)), INTENT(in) ::albDecTr_id
1336  ! REAL(KIND(1d0)), INTENT(in) ::DecidCap_id
1337  REAL(KIND(1d0)), INTENT(in) ::albEveTr_id
1338  REAL(KIND(1d0)), INTENT(in) ::albGrass_id
1339 
1340  ! REAL(KIND(1d0)), DIMENSION(6, nsurf), INTENT(inout)::StoreDrainPrm
1341 
1342  REAL(KIND(1d0)), DIMENSION(nsurf), INTENT(in)::SnowPack_prev
1343  REAL(KIND(1d0)), DIMENSION(nsurf), INTENT(in)::snowFrac_prev
1344  REAL(KIND(1d0)), DIMENSION(nsurf), INTENT(out)::snowFrac_next
1345  REAL(KIND(1d0)), DIMENSION(nsurf)::SnowFrac
1346 
1347  REAL(KIND(1d0)), INTENT(out)::ldown
1348  REAL(KIND(1d0)), INTENT(out)::fcld
1349  REAL(KIND(1d0)), INTENT(out)::qn1
1350  REAL(KIND(1d0)), INTENT(out)::qn1_snowfree
1351  REAL(KIND(1d0)), INTENT(out)::qn1_S
1352  REAL(KIND(1d0)), INTENT(out)::kclear
1353  REAL(KIND(1d0)), INTENT(out)::kup
1354  REAL(KIND(1d0)), INTENT(out)::lup
1355  REAL(KIND(1d0)), INTENT(out)::tsurf
1356  REAL(KIND(1d0)), INTENT(out)::alb1
1357  REAL(KIND(1d0)), INTENT(out)::SnowAlb_next
1358  REAL(KIND(1d0))::alb0
1359  REAL(KIND(1d0))::SnowAlb
1360 
1361  REAL(KIND(1d0)), DIMENSION(nsurf), INTENT(out) ::qn1_ind_snow
1362  REAL(KIND(1d0)), DIMENSION(nsurf), INTENT(out) ::kup_ind_snow
1363  REAL(KIND(1d0)), DIMENSION(nsurf), INTENT(out) ::Tsurf_ind_snow
1364  REAL(KIND(1d0)), DIMENSION(nsurf), INTENT(out):: tsurf_ind
1365 
1366  REAL(KIND(1d0)), DIMENSION(nsurf):: lup_ind
1367  REAL(KIND(1d0)), DIMENSION(nsurf):: kup_ind
1368  REAL(KIND(1d0)), DIMENSION(nsurf):: qn1_ind
1369 
1370  REAL(KIND(1d0)), PARAMETER::NAN = -999
1371  INTEGER :: NetRadiationMethod_use
1372  INTEGER::AlbedoChoice, ldown_option
1373 
1374  ! translate values
1375  alb = alb_prev
1376 
1377  ! update snow albedo
1378  snowalb = update_snow_albedo( &
1379  tstep, snowpack_prev, snowalb_prev, temp_c, &
1380  tau_a, tau_f, snowalbmax, snowalbmin)
1381 
1382  CALL radmethod( &
1383  netradiationmethod, &!input
1384  snowuse, &!input
1385  netradiationmethod_use, albedochoice, ldown_option)!output
1386 
1387  snowfrac = snowfrac_prev
1388  IF (netradiationmethod_use > 0) THEN
1389 
1390  ! IF (snowUse==0) SnowFrac=snowFrac_obs
1391  IF (snowuse == 0) snowfrac = 0
1392 
1393  IF (ldown_option == 1) THEN !Observed ldown provided as forcing
1394  ldown = ldown_obs
1395  ELSE
1396  ldown = -9 !to be filled in NARP
1397  ENDIF
1398 
1399  IF (ldown_option == 2) THEN !observed cloud fraction provided as forcing
1400  fcld = fcld_obs
1401  ENDIF
1402 
1403  !write(*,*) DecidCap(id), id, it, imin, 'Calc - near start'
1404 
1405  ! Update variables that change daily and represent seasonal variability
1406  alb(decidsurf) = albdectr_id !Change deciduous albedo
1407  ! StoreDrainPrm(6, DecidSurf) = DecidCap_id !Change current storage capacity of deciduous trees
1408  ! Change EveTr and Grass albedo too
1409  alb(conifsurf) = albevetr_id
1410  alb(grasssurf) = albgrass_id
1411 
1412  IF (diagnose == 1) WRITE (*, *) 'Calling NARP...'
1413  IF (diagqn == 1) WRITE (*, *) 'NetRadiationMethodX:', netradiationmethod_use
1414  IF (diagqn == 1) WRITE (*, *) 'AlbedoChoice:', albedochoice
1415 
1416  CALL narp( &
1417  nsurf, sfr, snowfrac, alb, emis, icefrac, &! input:
1418  narp_trans_site, narp_emis_snow, &
1419  dectime, zenith_deg, tsurf_0, avkdn, temp_c, avrh, ea_hpa, qn1_obs, &
1420  snowalb, &
1421  albedochoice, ldown_option, netradiationmethod_use, diagqn, &
1422  qn1, qn1_snowfree, qn1_s, kclear, kup, ldown, lup, fcld, tsurf, &! output:
1423  qn1_ind_snow, kup_ind_snow, tsurf_ind_snow, tsurf_ind, alb0, alb1)
1424 
1425  ELSE ! NetRadiationMethod==0
1426  snowfrac = snowfrac_obs
1427  qn1 = qn1_obs
1428  qn1_snowfree = qn1_obs
1429  qn1_s = qn1_obs
1430  ldown = nan
1431  lup = nan
1432  kup = nan
1433  tsurf = nan
1434  lup_ind = nan
1435  kup_ind = nan
1436  tsurf_ind = nan
1437  qn1_ind = nan
1438  fcld = nan
1439  ENDIF
1440  snowfrac_next = snowfrac
1441 
1442  IF (ldown_option == 1) THEN
1443  fcld = nan
1444  ENDIF
1445 
1446  ! translate values
1447  alb_next = alb
1448  snowalb_next = snowalb
1449 
1450  END SUBROUTINE suews_cal_qn
1451  !========================================================================
1452 
1453  !=============storage heat flux=========================================
1454  SUBROUTINE suews_cal_qs( &
1455  StorageHeatMethod, qs_obs, OHMIncQF, Gridiv, &!input
1456  id, tstep, dt_since_start, Diagnose, sfr, &
1457  OHM_coef, OHM_threshSW, OHM_threshWD, &
1458  soilstore_id, SoilStoreCap, state_id, SnowUse, SnowFrac, DiagQS, &
1459  HDD_id, MetForcingData_grid, Ts5mindata_ir, qf, qn1, &
1460  avkdn, avu1, temp_c, zenith_deg, avrh, press_hpa, ldown, &
1461  bldgh, alb, emis, cpAnOHM, kkAnOHM, chAnOHM, EmissionsMethod, &
1462  Tair_av, qn1_av_prev, dqndt_prev, qn1_s_av_prev, dqnsdt_prev, &
1463  StoreDrainPrm, &
1464  qn1_S, dataOutLineESTM, qs, &!output
1465  qn1_av_next, dqndt_next, qn1_s_av_next, dqnsdt_next, &
1466  deltaQi, a1, a2, a3)
1468  IMPLICIT NONE
1469 
1470  INTEGER, INTENT(in) ::StorageHeatMethod
1471  INTEGER, INTENT(in) ::OHMIncQF
1472  INTEGER, INTENT(in) ::Gridiv
1473  INTEGER, INTENT(in) ::id
1474  INTEGER, INTENT(in) ::tstep ! time step [s]
1475  INTEGER, INTENT(in) ::dt_since_start ! time since simulation starts [s]
1476  INTEGER, INTENT(in) ::Diagnose
1477  ! INTEGER, INTENT(in) ::nsh ! number of timesteps in one hour
1478  INTEGER, INTENT(in) ::SnowUse ! option for snow related calculations
1479  INTEGER, INTENT(in) ::DiagQS ! diagnostic option
1480  INTEGER, INTENT(in) :: EmissionsMethod
1481 
1482  REAL(KIND(1d0)), INTENT(in)::OHM_coef(nsurf + 1, 4, 3) ! OHM coefficients
1483  REAL(KIND(1d0)), INTENT(in)::OHM_threshSW(nsurf + 1) ! OHM thresholds
1484  REAL(KIND(1d0)), INTENT(in)::OHM_threshWD(nsurf + 1) ! OHM thresholds
1485  REAL(KIND(1d0)), INTENT(in)::soilstore_id(nsurf) ! soil moisture
1486  REAL(KIND(1d0)), INTENT(in)::SoilStoreCap(nsurf) ! capacity of soil store
1487  REAL(KIND(1d0)), INTENT(in)::state_id(nsurf) ! wetness status
1488 
1489  REAL(KIND(1d0)), DIMENSION(12), INTENT(in)::HDD_id
1490  REAL(KIND(1d0)), INTENT(in)::qf
1491  REAL(KIND(1d0)), INTENT(in)::qn1
1492  REAL(KIND(1d0)), INTENT(in)::qs_obs
1493  REAL(KIND(1d0)), INTENT(in)::avkdn, avu1, temp_c, zenith_deg, avrh, press_hpa, ldown
1494  REAL(KIND(1d0)), INTENT(in)::bldgh
1495 
1496  REAL(KIND(1d0)), DIMENSION(nsurf), INTENT(in)::sfr
1497  REAL(KIND(1d0)), DIMENSION(nsurf), INTENT(in)::alb
1498  REAL(KIND(1d0)), DIMENSION(nsurf), INTENT(in)::emis
1499  REAL(KIND(1d0)), DIMENSION(nsurf), INTENT(in)::cpAnOHM
1500  REAL(KIND(1d0)), DIMENSION(nsurf), INTENT(in)::kkAnOHM
1501  REAL(KIND(1d0)), DIMENSION(nsurf), INTENT(in)::chAnOHM
1502  REAL(KIND(1d0)), DIMENSION(nsurf), INTENT(in)::SnowFrac
1503 
1504  REAL(KIND(1d0)), DIMENSION(:, :), INTENT(in)::MetForcingData_grid
1505 
1506  REAL(KIND(1d0)), DIMENSION(:), INTENT(in)::Ts5mindata_ir
1507 
1508  REAL(KIND(1d0)), INTENT(in)::Tair_av
1509  REAL(KIND(1d0)), INTENT(in)::qn1_av_prev
1510  REAL(KIND(1d0)), INTENT(out)::qn1_av_next
1511  REAL(KIND(1d0)), INTENT(in)::dqndt_prev ! Rate of change of net radiation [W m-2 h-1] at t-1
1512  REAL(KIND(1d0)), INTENT(out)::dqndt_next ! Rate of change of net radiation [W m-2 h-1] at t-1
1513  REAL(KIND(1d0)), INTENT(in)::qn1_s_av_prev
1514  REAL(KIND(1d0)), INTENT(out)::qn1_s_av_next
1515  REAL(KIND(1d0)), INTENT(in)::dqnsdt_prev ! Rate of change of net radiation [W m-2 h-1] at t-1
1516  REAL(KIND(1d0)), INTENT(out)::dqnsdt_next ! Rate of change of net radiation [W m-2 h-1] at t-1
1517  ! REAL(KIND(1d0)),DIMENSION(nsh),INTENT(inout) ::qn1_store_grid
1518  ! REAL(KIND(1d0)),DIMENSION(nsh),INTENT(inout) ::qn1_S_store_grid !< stored qn1 [W m-2]
1519 
1520  ! REAL(KIND(1d0)),DIMENSION(2*nsh+1),INTENT(inout)::qn1_av_store_grid
1521  ! REAL(KIND(1d0)),DIMENSION(2*nsh+1),INTENT(inout)::qn1_S_av_store_grid !< average net radiation over previous hour [W m-2]
1522  REAL(KIND(1d0)), DIMENSION(6, nsurf), INTENT(in)::StoreDrainPrm
1523 
1524  REAL(KIND(1d0)), DIMENSION(nsurf), INTENT(out)::deltaQi ! storage heat flux of snow surfaces
1525 
1526  REAL(KIND(1d0)), DIMENSION(27), INTENT(out):: dataOutLineESTM
1527  REAL(KIND(1d0)), INTENT(out)::qn1_S
1528  REAL(KIND(1d0)), INTENT(out):: qs ! storage heat flux
1529  REAL(KIND(1d0)), INTENT(out):: a1
1530  REAL(KIND(1d0)), INTENT(out):: a2
1531  REAL(KIND(1d0)), INTENT(out):: a3
1532 
1533  REAL(KIND(1d0))::Tair_mav_5d ! Tair_mav_5d=HDD(id-1,4) HDD at the begining of today (id-1)
1534  REAL(KIND(1d0))::qn1_use ! qn used in OHM calculations
1535 
1536  REAL(KIND(1d0)):: moist_surf(nsurf)
1537 
1538  ! initialise output variables
1539  !deltaQi = 0
1540  !SnowFrac = 0
1541  !qn1_S = 0
1542  dataoutlineestm = -999
1543  qs = -999
1544  a1 = -999
1545  a2 = -999
1546  a3 = -999
1547 
1548  ! calculate qn if qf should be included
1549  IF (ohmincqf == 1) THEN
1550  qn1_use = qf + qn1
1551  ELSEIF (ohmincqf == 0) THEN
1552  qn1_use = qn1
1553  ENDIF
1554 
1555  IF (storageheatmethod == 0) THEN !Use observed QS
1556  qs = qs_obs
1557 
1558  ELSEIF (storageheatmethod == 1) THEN !Use OHM to calculate QS
1559  tair_mav_5d = hdd_id(10)
1560  IF (diagnose == 1) WRITE (*, *) 'Calling OHM...'
1561  CALL ohm(qn1_use, qn1_av_prev, dqndt_prev, qn1_av_next, dqndt_next, &
1562  qn1_s, qn1_s_av_prev, dqnsdt_prev, qn1_s_av_next, dqnsdt_next, &
1563  tstep, dt_since_start, &
1564  sfr, nsurf, &
1565  tair_mav_5d, &
1566  ohm_coef, &
1567  ohm_threshsw, ohm_threshwd, &
1568  soilstore_id, soilstorecap, state_id, &
1569  bldgsurf, watersurf, &
1570  snowuse, snowfrac, &
1571  diagqs, &
1572  a1, a2, a3, qs, deltaqi)
1573 
1574  ! use AnOHM to calculate QS, TS 14 Mar 2016
1575  ELSEIF (storageheatmethod == 3) THEN
1576  IF (diagnose == 1) WRITE (*, *) 'Calling AnOHM...'
1577  ! CALL AnOHM(qn1_use,qn1_store_grid,qn1_av_store_grid,qf,&
1578  ! MetForcingData_grid,state_id/StoreDrainPrm(6,:),&
1579  ! alb, emis, cpAnOHM, kkAnOHM, chAnOHM,&
1580  ! sfr,nsurf,nsh,EmissionsMethod,id,Gridiv,&
1581  ! a1,a2,a3,qs,deltaQi)
1582  moist_surf = state_id/storedrainprm(6, :)
1583  CALL anohm( &
1584  tstep, dt_since_start, &
1585  qn1_use, qn1_av_prev, dqndt_prev, qf, &
1586  metforcingdata_grid, moist_surf, &
1587  alb, emis, cpanohm, kkanohm, chanohm, &! input
1588  sfr, nsurf, emissionsmethod, id, gridiv, &
1589  qn1_av_next, dqndt_next, &
1590  a1, a2, a3, qs, deltaqi)! output
1591 
1592  ! !Calculate QS using ESTM
1593  ELSEIF (storageheatmethod == 4 .OR. storageheatmethod == 14) THEN
1594  ! !CALL ESTM(QSestm,iMB)
1595  IF (diagnose == 1) WRITE (*, *) 'Calling ESTM...'
1596  CALL estm( &
1597  gridiv, &!input
1598  tstep, &
1599  avkdn, avu1, temp_c, zenith_deg, avrh, press_hpa, ldown, &
1600  bldgh, ts5mindata_ir, &
1601  tair_av, &
1602  dataoutlineestm, qs)!output
1603  ! CALL ESTM(QSestm,Gridiv,ir) ! iMB corrected to Gridiv, TS 09 Jun 2016
1604  ! QS=QSestm ! Use ESTM qs
1605  ENDIF
1606 
1607  END SUBROUTINE suews_cal_qs
1608  !=======================================================================
1609 
1610  !==========================water balance================================
1611  SUBROUTINE suews_cal_water( &
1612  Diagnose, &!input
1613  snowUse, NonWaterFraction, addPipes, addImpervious, addVeg, addWaterBody, &
1614  state_id, soilstore_id, sfr, StoreDrainPrm, WaterDist, nsh_real, &
1615  drain_per_tstep, & !output
1616  drain, AddWaterRunoff, &
1617  AdditionalWater, runoffPipes, runoff_per_interval, &
1618  AddWater, stateOld, soilstoreOld)
1620  IMPLICIT NONE
1621  ! INTEGER,PARAMETER :: nsurf=7! number of surface types
1622  ! INTEGER,PARAMETER ::WaterSurf = 7
1623  INTEGER, INTENT(in) ::Diagnose
1624  INTEGER, INTENT(in) ::snowUse
1625 
1626  REAL(KIND(1d0)), INTENT(in)::NonWaterFraction
1627  REAL(KIND(1d0)), INTENT(in)::addPipes
1628  REAL(KIND(1d0)), INTENT(in)::addImpervious
1629  REAL(KIND(1d0)), INTENT(in)::addVeg
1630  REAL(KIND(1d0)), INTENT(in)::addWaterBody
1631  REAL(KIND(1d0)), INTENT(in)::nsh_real !nsh cast as a real for use in calculations
1632 
1633  REAL(KIND(1d0)), DIMENSION(nsurf), INTENT(in) ::state_id
1634  REAL(KIND(1d0)), DIMENSION(nsurf), INTENT(in) ::soilstore_id
1635  REAL(KIND(1d0)), DIMENSION(nsurf), INTENT(in) ::sfr
1636  REAL(KIND(1d0)), DIMENSION(6, nsurf), INTENT(in) ::StoreDrainPrm
1637  REAL(KIND(1d0)), DIMENSION(nsurf + 1, nsurf - 1), INTENT(in)::WaterDist
1638 
1639  REAL(KIND(1d0)), DIMENSION(nsurf), INTENT(out):: drain !Drainage of surface type "is" [mm]
1640  REAL(KIND(1d0)), DIMENSION(nsurf), INTENT(out):: AddWaterRunoff!Fraction of water going to runoff/sub-surface soil (WGWaterDist) [-]
1641  REAL(KIND(1d0)), DIMENSION(nsurf), INTENT(out):: AddWater !water from other surfaces (WGWaterDist in SUEWS_ReDistributeWater.f95) [mm]
1642  REAL(KIND(1d0)), DIMENSION(nsurf), INTENT(out):: stateOld
1643  REAL(KIND(1d0)), DIMENSION(nsurf), INTENT(out):: soilstoreOld
1644 
1645  REAL(KIND(1d0)), INTENT(out)::drain_per_tstep
1646  REAL(KIND(1d0)), INTENT(out)::AdditionalWater
1647  REAL(KIND(1d0)), INTENT(out)::runoffPipes
1648  REAL(KIND(1d0)), INTENT(out)::runoff_per_interval
1649  INTEGER:: is
1650 
1651  ! Retain previous surface state_id and soil moisture state_id
1652  stateold = state_id !state_id of each surface [mm] for the previous timestep
1653  soilstoreold = soilstore_id !Soil moisture of each surface [mm] for the previous timestep
1654 
1655  !============= Grid-to-grid runoff =============
1656  ! Calculate additional water coming from other grids
1657  ! i.e. the variables addImpervious, addVeg, addWaterBody, addPipes
1658  !call RunoffFromGrid(GridFromFrac) !!Need to code between-grid water transfer
1659 
1660  ! Sum water coming from other grids (these are expressed as depths over the whole surface)
1661  additionalwater = addpipes + addimpervious + addveg + addwaterbody ![mm]
1662 
1663  ! Initialise runoff in pipes
1664  runoffpipes = addpipes !Water flowing in pipes from other grids. QUESTION: No need for scaling?
1665  !! CHECK p_i
1666  runoff_per_interval = addpipes !pipe plor added to total runoff.
1667 
1668  !================== Drainage ===================
1669  ! Calculate drainage for each soil subsurface (excluding water body)
1670  IF (diagnose == 1) WRITE (*, *) 'Calling Drainage...'
1671 
1672  IF (nonwaterfraction /= 0) THEN !Soil states only calculated if soil exists. LJ June 2017
1673  DO is = 1, nsurf - 1
1674 
1675  CALL drainage( &
1676  is, &! input:
1677  state_id(is), &
1678  storedrainprm(6, is), &
1679  storedrainprm(2, is), &
1680  storedrainprm(3, is), &
1681  storedrainprm(4, is), &
1682  nsh_real, &
1683  drain(is))! output
1684 
1685  ! !HCW added and changed to StoreDrainPrm(6,is) here 20 Feb 2015
1686  ! drain_per_tstep=drain_per_tstep+(drain(is)*sfr(is)/NonWaterFraction) !No water body included
1687  ENDDO
1688  drain_per_tstep = dot_product(drain(1:nsurf - 1), sfr(1:nsurf - 1))/nonwaterfraction !No water body included
1689  ELSE
1690  drain(1:nsurf - 1) = 0
1691  drain_per_tstep = 0
1692  ENDIF
1693 
1694  drain(watersurf) = 0 ! Set drainage from water body to zero
1695 
1696  ! Distribute water within grid, according to WithinGridWaterDist matrix (Cols 1-7)
1697  IF (diagnose == 1) WRITE (*, *) 'Calling ReDistributeWater...'
1698  ! CALL ReDistributeWater
1699  !Calculates AddWater(is)
1700  CALL redistributewater( &
1701  snowuse, waterdist, sfr, drain, &! input:
1702  addwaterrunoff, addwater)! output
1703 
1704  END SUBROUTINE suews_cal_water
1705  !=======================================================================
1706 
1707  !===============initialize sensible heat flux============================
1708  SUBROUTINE suews_init_qh( &
1709  avdens, avcp, h_mod, qn1, dectime, &!input
1710  H_init)!output
1712  IMPLICIT NONE
1713  ! REAL(KIND(1d0)), INTENT(in)::qh_obs
1714  REAL(KIND(1d0)), INTENT(in)::avdens
1715  REAL(KIND(1d0)), INTENT(in)::avcp
1716  REAL(KIND(1d0)), INTENT(in)::h_mod
1717  REAL(KIND(1d0)), INTENT(in)::qn1
1718  REAL(KIND(1d0)), INTENT(in)::dectime
1719  REAL(KIND(1d0)), INTENT(out)::H_init
1720 
1721  REAL(KIND(1d0)), PARAMETER::NAN = -999
1722  INTEGER, PARAMETER::notUsedI = -999
1723 
1724  ! Calculate kinematic heat flux (w'T') from sensible heat flux [W m-2] from observed data (if available) or LUMPS
1725  ! IF (qh_obs /= NAN) THEN !if(qh_obs/=NAN) qh=qh_obs !Commented out by HCW 04 Mar 2015
1726  ! H_init = qh_obs/(avdens*avcp) !Use observed value
1727  ! ELSE
1728  IF (h_mod /= nan) THEN
1729  h_init = h_mod/(avdens*avcp) !Use LUMPS value
1730  ELSE
1731  h_init = (qn1*0.2)/(avdens*avcp) !If LUMPS has had a problem, we still need a value
1732  CALL errorhint(38, 'LUMPS unable to calculate realistic value for H_mod.', h_mod, dectime, notusedi)
1733  ENDIF
1734  ! ENDIF
1735 
1736  END SUBROUTINE suews_init_qh
1737  !========================================================================
1738 
1739  !================latent heat flux and surface wetness===================
1740  ! TODO: optimise the structure of this function
1741  SUBROUTINE suews_cal_qe( &
1742  Diagnose, snowuse, &!input
1743  tstep, imin, it, EvapMethod, snowCalcSwitch, dayofWeek_id, CRWmin, CRWmax, &
1744  dectime, lvS_J_kg, avRh, Press_hPa, Temp_C, &
1745  RAsnow, psyc_hPa, sIce_hPa, &
1746  PervFraction, vegfraction, addimpervious, qn1_snowfree, qf, qs, vpd_hPa, s_hPa, &
1747  ResistSurf, RA, rb, snowdensmin, precip, PipeCapacity, RunoffToWater, &
1748  NonWaterFraction, wu_EveTr, wu_DecTr, wu_Grass, addVeg, addWaterBody, SnowLimPaved, SnowLimBldg, &
1749  SurfaceArea, FlowChange, drain, WetThresh, stateOld, mw_ind, SoilStoreCap, rainonsnow, &
1750  freezmelt, freezstate, freezstatevol, Qm_Melt, Qm_rain, Tsurf_ind, sfr, &
1751  StateLimit, AddWater, addwaterrunoff, StoreDrainPrm, SnowPackLimit, SnowProf_24hr, &
1752  SnowPack_in, SnowFrac_in, SnowWater_in, iceFrac_in, SnowDens_in, &! input:
1753  runoff_per_interval_in, state_id_in, soilstore_id_in, &! input:
1754  state_id_out, soilstore_id_out, &! output:
1755  SnowPack_out, SnowFrac_out, SnowWater_out, iceFrac_out, SnowDens_out, &! output
1756  runoffSoil, &! output:
1757  SnowRemoval, &
1758  state_per_tstep, NWstate_per_tstep, qe, &
1759  swe, chSnow_per_interval, ev_per_tstep, runoff_per_tstep, &
1760  surf_chang_per_tstep, runoffPipes, mwstore, runoffwaterbody, &
1761  runoffAGveg, runoffAGimpervious)
1763  IMPLICIT NONE
1764 
1765  INTEGER, INTENT(in) ::Diagnose
1766  INTEGER, INTENT(in) ::snowuse
1767  INTEGER, INTENT(in) ::tstep
1768  INTEGER, INTENT(in) ::imin
1769  INTEGER, INTENT(in) ::it
1770  INTEGER, INTENT(in) ::EvapMethod !Evaporation calculated according to Rutter (1) or Shuttleworth (2)
1771 
1772  INTEGER, DIMENSION(nsurf), INTENT(in)::snowCalcSwitch
1773  INTEGER, DIMENSION(3), INTENT(in)::dayofWeek_id
1774 
1775  REAL(KIND(1d0)), INTENT(in)::CRWmin
1776  REAL(KIND(1d0)), INTENT(in)::CRWmax
1777  REAL(KIND(1d0)), INTENT(in)::dectime
1778  REAL(KIND(1d0)), INTENT(in)::lvS_J_kg
1779  ! REAL(KIND(1d0)), INTENT(in)::lv_j_kg
1780  ! REAL(KIND(1d0)), INTENT(in)::avdens
1781  REAL(KIND(1d0)), INTENT(in)::avRh
1782  REAL(KIND(1d0)), INTENT(in)::Press_hPa
1783  REAL(KIND(1d0)), INTENT(in)::Temp_C
1784  REAL(KIND(1d0)), INTENT(in)::RAsnow
1785  REAL(KIND(1d0)), INTENT(in)::psyc_hPa
1786  ! REAL(KIND(1d0)), INTENT(in)::avcp
1787  REAL(KIND(1d0)), INTENT(in)::sIce_hPa
1788  REAL(KIND(1d0)), INTENT(in)::PervFraction
1789  REAL(KIND(1d0)), INTENT(in)::vegfraction
1790  REAL(KIND(1d0)), INTENT(in)::addimpervious
1791  REAL(KIND(1d0)), INTENT(in)::qn1_snowfree
1792  REAL(KIND(1d0)), INTENT(in)::qf
1793  REAL(KIND(1d0)), INTENT(in)::qs
1794  REAL(KIND(1d0)), INTENT(in)::vpd_hPa
1795  REAL(KIND(1d0)), INTENT(in)::s_hPa
1796  REAL(KIND(1d0)), INTENT(in)::ResistSurf
1797  REAL(KIND(1d0)), INTENT(in)::RA
1798  REAL(KIND(1d0)), INTENT(in)::rb
1799  REAL(KIND(1d0)), INTENT(in)::snowdensmin
1800  REAL(KIND(1d0)), INTENT(in)::precip
1801  REAL(KIND(1d0)), INTENT(in)::PipeCapacity
1802  REAL(KIND(1d0)), INTENT(in)::RunoffToWater
1803  REAL(KIND(1d0)), INTENT(in)::NonWaterFraction
1804  REAL(KIND(1d0)), INTENT(in)::wu_EveTr!Water use for evergreen trees/shrubs [mm]
1805  REAL(KIND(1d0)), INTENT(in)::wu_DecTr!Water use for deciduous trees/shrubs [mm]
1806  REAL(KIND(1d0)), INTENT(in)::wu_Grass!Water use for grass [mm]
1807  REAL(KIND(1d0)), INTENT(in)::addVeg!Water from vegetated surfaces of other grids [mm] for whole surface area
1808  REAL(KIND(1d0)), INTENT(in)::addWaterBody!Water from water surface of other grids [mm] for whole surface area
1809  REAL(KIND(1d0)), INTENT(in)::SnowLimPaved
1810  REAL(KIND(1d0)), INTENT(in)::SnowLimBldg
1811  REAL(KIND(1d0)), INTENT(in)::SurfaceArea
1812  REAL(KIND(1d0)), INTENT(in)::FlowChange
1813 
1814  REAL(KIND(1d0)), DIMENSION(nsurf), INTENT(in)::drain
1815  REAL(KIND(1d0)), DIMENSION(nsurf), INTENT(in)::WetThresh
1816  REAL(KIND(1d0)), DIMENSION(nsurf), INTENT(in)::stateOld
1817  REAL(KIND(1d0)), DIMENSION(nsurf), INTENT(in)::mw_ind
1818  REAL(KIND(1d0)), DIMENSION(nsurf), INTENT(in)::SoilStoreCap
1819  REAL(KIND(1d0)), DIMENSION(nsurf), INTENT(in)::rainonsnow
1820  REAL(KIND(1d0)), DIMENSION(nsurf), INTENT(in)::freezmelt
1821  REAL(KIND(1d0)), DIMENSION(nsurf), INTENT(in)::freezstate
1822  REAL(KIND(1d0)), DIMENSION(nsurf), INTENT(in)::freezstatevol
1823  REAL(KIND(1d0)), DIMENSION(nsurf), INTENT(in)::Qm_Melt
1824  REAL(KIND(1d0)), DIMENSION(nsurf), INTENT(in)::Qm_rain
1825  REAL(KIND(1d0)), DIMENSION(nsurf), INTENT(in)::Tsurf_ind
1826  REAL(KIND(1d0)), DIMENSION(nsurf), INTENT(in)::sfr
1827  REAL(KIND(1d0)), DIMENSION(nsurf), INTENT(in)::SnowPackLimit
1828  REAL(KIND(1d0)), DIMENSION(nsurf), INTENT(in)::StateLimit !Limit for state_id of each surface type [mm] (specified in input files)
1829  REAL(KIND(1d0)), DIMENSION(nsurf), INTENT(in)::AddWater
1830  REAL(KIND(1d0)), DIMENSION(nsurf), INTENT(in)::addwaterrunoff
1831  REAL(KIND(1d0)), DIMENSION(6, nsurf), INTENT(in)::StoreDrainPrm
1832  REAL(KIND(1d0)), DIMENSION(0:23, 2), INTENT(in):: SnowProf_24hr
1833 
1834  ! Total water transported to each grid for grid-to-grid connectivity
1835  REAL(KIND(1d0)), INTENT(in)::runoff_per_interval_in
1836  REAL(KIND(1d0)), DIMENSION(nsurf), INTENT(in)::state_id_in
1837  REAL(KIND(1d0)), DIMENSION(nsurf), INTENT(in)::soilstore_id_in
1838  REAL(KIND(1d0)), DIMENSION(nsurf), INTENT(in)::SnowPack_in
1839  REAL(KIND(1d0)), DIMENSION(nsurf), INTENT(in)::SnowFrac_in
1840  REAL(KIND(1d0)), DIMENSION(nsurf), INTENT(in)::SnowWater_in
1841  REAL(KIND(1d0)), DIMENSION(nsurf), INTENT(in)::iceFrac_in
1842  REAL(KIND(1d0)), DIMENSION(nsurf), INTENT(in)::SnowDens_in
1843 
1844  ! output:
1845  REAL(KIND(1d0))::runoff_per_interval_out
1846  REAL(KIND(1d0)), DIMENSION(nsurf), INTENT(out)::state_id_out
1847  REAL(KIND(1d0)), DIMENSION(nsurf), INTENT(out)::soilstore_id_out
1848  REAL(KIND(1d0)), DIMENSION(nsurf), INTENT(out)::SnowPack_out
1849  REAL(KIND(1d0)), DIMENSION(nsurf), INTENT(out)::SnowFrac_out
1850  REAL(KIND(1d0)), DIMENSION(nsurf), INTENT(out)::SnowWater_out
1851  REAL(KIND(1d0)), DIMENSION(nsurf), INTENT(out)::iceFrac_out
1852  REAL(KIND(1d0)), DIMENSION(nsurf), INTENT(out)::SnowDens_out
1853 
1854  REAL(KIND(1d0)), DIMENSION(nsurf)::runoffSnow !Initialize for runoff caused by snowmelting
1855  REAL(KIND(1d0)), DIMENSION(nsurf)::runoff
1856  REAL(KIND(1d0)), DIMENSION(nsurf), INTENT(out)::runoffSoil
1857  REAL(KIND(1d0)), DIMENSION(nsurf)::chang
1858  REAL(KIND(1d0)), DIMENSION(nsurf)::changSnow
1859  REAL(KIND(1d0)), DIMENSION(nsurf)::snowDepth
1860  REAL(KIND(1d0)), DIMENSION(nsurf)::SnowToSurf
1861  REAL(KIND(1d0)), DIMENSION(nsurf)::ev_snow
1862  REAL(KIND(1d0)), DIMENSION(2), INTENT(out)::SnowRemoval
1863  REAL(KIND(1d0)), DIMENSION(nsurf)::evap
1864  REAL(KIND(1d0)), DIMENSION(nsurf)::rss_nsurf
1865 
1866  REAL(KIND(1d0))::p_mm!Inputs to surface water balance
1867  ! REAL(KIND(1d0)),INTENT(out)::rss
1868  REAL(KIND(1d0))::qe_surf ! latent heat flux [W m-2]
1869  REAL(KIND(1d0)), INTENT(out)::state_per_tstep
1870  REAL(KIND(1d0)), INTENT(out)::NWstate_per_tstep
1871  REAL(KIND(1d0)), INTENT(out)::qe
1872  REAL(KIND(1d0)), INTENT(out)::swe
1873  REAL(KIND(1d0))::ev
1874  REAL(KIND(1d0)), INTENT(out)::chSnow_per_interval
1875  REAL(KIND(1d0)), INTENT(out)::ev_per_tstep
1876  REAL(KIND(1d0))::qe_per_tstep
1877  REAL(KIND(1d0)), INTENT(out)::runoff_per_tstep
1878  REAL(KIND(1d0)), INTENT(out)::surf_chang_per_tstep
1879  REAL(KIND(1d0)), INTENT(out)::runoffPipes
1880  REAL(KIND(1d0)), INTENT(out)::mwstore
1881  REAL(KIND(1d0)), INTENT(out)::runoffwaterbody
1882  REAL(KIND(1d0))::runoffWaterBody_m3
1883  REAL(KIND(1d0))::runoffPipes_m3
1884  REAL(KIND(1d0)), INTENT(out)::runoffAGveg
1885  REAL(KIND(1d0)), INTENT(out)::runoffAGimpervious
1886 
1887  ! local:
1888  INTEGER:: is
1889 
1890  REAL(KIND(1d0))::runoff_per_interval
1891  REAL(KIND(1d0)), DIMENSION(nsurf)::state_id
1892  REAL(KIND(1d0)), DIMENSION(nsurf)::soilstore_id
1893  REAL(KIND(1d0)), DIMENSION(nsurf)::SnowPack
1894  REAL(KIND(1d0)), DIMENSION(nsurf)::SnowFrac
1895  REAL(KIND(1d0)), DIMENSION(nsurf)::SnowWater
1896  REAL(KIND(1d0)), DIMENSION(nsurf)::iceFrac
1897  REAL(KIND(1d0)), DIMENSION(nsurf)::SnowDens
1898 
1899  REAL(KIND(1d0)), DIMENSION(2) ::SurplusEvap !Surplus for evaporation in 5 min timestep
1900  REAL(KIND(1d0))::surplusWaterBody
1901  REAL(KIND(1d0))::pin!Rain per time interval
1902  ! REAL(KIND(1d0))::sae
1903  ! REAL(KIND(1d0))::vdrc
1904  ! REAL(KIND(1d0))::sp
1905  ! REAL(KIND(1d0))::numPM
1906  REAL(KIND(1d0))::qn_e
1907  REAL(KIND(1d0))::tlv
1908  REAL(KIND(1d0))::runoffAGimpervious_m3
1909  REAL(KIND(1d0))::runoffAGveg_m3
1910  REAL(KIND(1d0))::nsh_real
1911  REAL(KIND(1d0))::tstep_real
1912  REAL(KIND(1d0))::ev_tot
1913  REAL(KIND(1d0))::qe_tot
1914  REAL(KIND(1d0))::surf_chang_tot
1915  REAL(KIND(1d0))::runoff_tot
1916  REAL(KIND(1d0))::chSnow_tot
1917 
1918  REAL(KIND(1d0)), DIMENSION(7)::capStore ! current storage capacity [mm]
1919 
1920  runoff_per_interval = runoff_per_interval_in
1921  state_id = state_id_in
1922  soilstore_id = soilstore_id_in
1923  snowpack = snowpack_in
1924  snowfrac = snowfrac_in
1925  snowwater = snowwater_in
1926  icefrac = icefrac_in
1927  snowdens = snowdens_in
1928 
1929  tstep_real = tstep*1.d0
1930  nsh_real = 3600/tstep_real
1931 
1932  capstore = 0 !initialise capStore
1933 
1934  tlv = lv_j_kg/tstep_real !Latent heat of vapourisation per timestep
1935 
1936  pin = max(0., precip)!Initiate rain data [mm]
1937 
1938  ! Initialize the output variables
1939  qe_surf = 0
1940  ev = 0
1941  qe_tot = 0
1942  ev_tot = 0
1943  swe = 0
1944  ev_snow = 0
1945  ev_per_tstep = 0
1946  qe_per_tstep = 0
1947  surf_chang_per_tstep = 0
1948  runoff_per_tstep = 0
1949  state_per_tstep = 0
1950  nwstate_per_tstep = 0
1951  qe = 0
1952  runoffwaterbody = 0
1953  chsnow_per_interval = 0
1954  mwstore = 0
1955  runoffagveg = 0
1956  runoffagimpervious = 0
1957  surpluswaterbody = 0
1958  runoffsoil = 0
1959  runoff = 0
1960  chang = 0
1961  surplusevap = 0
1962  snowremoval = 0
1963 
1964  ! net available energy for evaporation
1965  qn_e = qn1_snowfree + qf - qs ! qn1 changed to qn1_snowfree, lj in May 2013
1966 
1967  IF (diagnose == 1) WRITE (*, *) 'Calling evap_SUEWS and SoilStore...'
1968  DO is = 1, nsurf !For each surface in turn
1969  IF (snowuse == 1 .and. snowcalcswitch(is) == 1) THEN ! snow calculation
1970  IF (sfr(is) /= 0) THEN
1971  ! IF (Diagnose == 1) WRITE (*, *) 'Calling SnowCalc...'
1972  CALL snowcalc( &
1973  tstep, imin, it, dectime, is, &!input
1974  evapmethod, crwmin, crwmax, nsh_real, lvs_j_kg, lv_j_kg, avdens, &
1975  avrh, press_hpa, temp_c, rasnow, psyc_hpa, avcp, sice_hpa, &
1976  pervfraction, vegfraction, addimpervious, &
1977  vpd_hpa, qn_e, s_hpa, resistsurf, ra, rb, tlv, snowdensmin, snowprof_24hr, precip, &
1978  pipecapacity, runofftowater, &
1979  addveg, snowlimpaved, snowlimbldg, flowchange, drain, &
1980  wetthresh, stateold, mw_ind, soilstorecap, rainonsnow, &
1981  freezmelt, freezstate, freezstatevol, &
1982  qm_melt, qm_rain, tsurf_ind, sfr, dayofweek_id, storedrainprm, snowpacklimit, &
1983  addwater, addwaterrunoff, &
1984  soilstore_id, snowpack, surplusevap, &!inout
1985  snowfrac, snowwater, icefrac, snowdens, &
1986  runoffagimpervious, runoffagveg, surpluswaterbody, &
1987  rss_nsurf, runoffsnow, & ! output
1988  runoff, runoffsoil, chang, changsnow, snowtosurf, state_id, ev_snow, &
1989  snowdepth, snowremoval, swe, ev, chsnow_tot, &
1990  ev_tot, qe_tot, runoff_tot, surf_chang_tot, &
1991  runoffpipes, mwstore, runoffwaterbody)
1992 
1993  !Actual updates here as xx_tstep variables not taken as input to snowcalc
1994  ev_per_tstep = ev_per_tstep + ev_tot
1995  qe_per_tstep = qe_per_tstep + qe_tot
1996  runoff_per_tstep = runoff_per_tstep + runoff_tot
1997  surf_chang_per_tstep = surf_chang_per_tstep + surf_chang_tot
1998  chsnow_per_interval = chsnow_per_interval + chsnow_tot
1999  ELSE
2000  snowfrac(is) = 0
2001  snowdens(is) = 0
2002  snowpack(is) = 0
2003  ENDIF
2004 
2005  !Store ev_tot for each surface
2006  evap(is) = ev_tot
2007  ELSE ! snow-free calculation
2008 
2009  capstore(is) = storedrainprm(6, is)
2010  !Calculates ev [mm]
2011  CALL evap_suews( &
2012  evapmethod, state_id(is), wetthresh(is), capstore(is), &!input
2013  vpd_hpa, avdens, avcp, qn_e, s_hpa, psyc_hpa, resistsurf, ra, rb, tlv, &
2014  rss_nsurf(is), ev, qe_surf) !output
2015 
2016  !Surface water balance and soil store updates (can modify ev, updates state_id)
2017  CALL soilstore( &
2018  is, sfr, pipecapacity, runofftowater, pin, & ! input:
2019  wu_evetr, wu_dectr, wu_grass, drain, addwater, addimpervious, nsh_real, stateold, addwaterrunoff, &
2020  pervfraction, addveg, soilstorecap, addwaterbody, flowchange, statelimit, runoffagimpervious, surpluswaterbody, &
2021  runoffagveg, runoffpipes, ev, soilstore_id, surplusevap, runoffwaterbody, &
2022  p_mm, chang, runoff, state_id)!output:
2023 
2024  !Store ev for each surface
2025  evap(is) = ev
2026 
2027  ! Sum evaporation from different surfaces to find total evaporation [mm]
2028  ev_per_tstep = ev_per_tstep + evap(is)*sfr(is)
2029 
2030  ! Sum latent heat flux from different surfaces to find total latent heat flux
2031  qe_per_tstep = qe_per_tstep + qe_surf*sfr(is)
2032 
2033  ! Sum change from different surfaces to find total change to surface state_id
2034  surf_chang_per_tstep = surf_chang_per_tstep + (state_id(is) - stateold(is))*sfr(is)
2035 
2036  ! Sum runoff from different surfaces to find total runoff
2037  runoff_per_tstep = runoff_per_tstep + runoff(is)*sfr(is)
2038 
2039  ! Calculate total state_id (including water body)
2040  state_per_tstep = state_per_tstep + state_id(is)*sfr(is)
2041 
2042  ! sum The total runoff from the area !!Check (HCW)
2043  runoff_per_interval = runoff_per_interval + (runoff(is)*sfr(is))
2044 
2045  IF (nonwaterfraction /= 0 .AND. is /= watersurf) THEN
2046  nwstate_per_tstep = nwstate_per_tstep + (state_id(is)*sfr(is)/nonwaterfraction)
2047  ENDIF
2048 
2049  changsnow(is) = 0
2050  runoffsnow(is) = 0
2051 
2052  ENDIF
2053  ENDDO !end loop over surfaces
2054 
2055  qe = qe_per_tstep
2056 
2057  ! Calculate volume of water that will move between grids
2058  ! Volume [m3] = Depth relative to whole area [mm] / 1000 [mm m-1] * SurfaceArea [m2]
2059  ! Need to use these volumes when converting back to addImpervious, AddVeg and AddWater
2060  runoffagimpervious_m3 = runoffagimpervious/1000*surfacearea
2061  runoffagveg_m3 = runoffagveg/1000*surfacearea
2062  runoffwaterbody_m3 = runoffwaterbody/1000*surfacearea
2063  runoffpipes_m3 = runoffpipes/1000*surfacearea
2064 
2065  runoff_per_interval_out = runoff_per_interval
2066  state_id_out = state_id
2067  soilstore_id_out = soilstore_id
2068  snowpack_out = snowpack
2069  snowfrac_out = snowfrac
2070  snowwater_out = snowwater
2071  icefrac_out = icefrac
2072  snowdens_out = snowdens
2073 
2074  END SUBROUTINE suews_cal_qe
2075  !========================================================================
2076 
2077  !===============sensible heat flux======================================
2078  SUBROUTINE suews_cal_qh( &
2079  QHMethod, &!input
2080  qn1, qf, QmRain, qeOut, qs, QmFreez, qm, avdens, avcp, tsurf, Temp_C, RA, &
2081  qh, qh_residual, qh_resist)!output
2082  IMPLICIT NONE
2083 
2084  INTEGER, INTENT(in) :: QHMethod ! option for QH calculation: 1, residual; 2, resistance-based
2085 
2086  REAL(KIND(1d0)), INTENT(in)::qn1
2087  REAL(KIND(1d0)), INTENT(in)::qf
2088  REAL(KIND(1d0)), INTENT(in)::QmRain
2089  REAL(KIND(1d0)), INTENT(in)::qeOut
2090  REAL(KIND(1d0)), INTENT(in)::qs
2091  REAL(KIND(1d0)), INTENT(in)::QmFreez
2092  REAL(KIND(1d0)), INTENT(in)::qm
2093  REAL(KIND(1d0)), INTENT(in)::avdens
2094  REAL(KIND(1d0)), INTENT(in)::avcp
2095  REAL(KIND(1d0)), INTENT(in)::tsurf
2096  REAL(KIND(1d0)), INTENT(in)::Temp_C
2097  REAL(KIND(1d0)), INTENT(in)::RA
2098 
2099  REAL(KIND(1d0)), INTENT(out)::qh
2100  REAL(KIND(1d0)), INTENT(out)::qh_resist
2101  REAL(KIND(1d0)), INTENT(out)::qh_residual
2102 
2103  REAL(KIND(1d0)), PARAMETER::NAN = -999
2104 
2105  ! Calculate sensible heat flux as a residual (Modified by LJ in Nov 2012)
2106  qh_residual = (qn1 + qf + qmrain) - (qeout + qs + qm + qmfreez) !qh=(qn1+qf+QmRain+QmFreez)-(qeOut+qs+Qm)
2107 
2108  ! ! Calculate QH using resistance method (for testing HCW 06 Jul 2016)
2109  ! Aerodynamic-Resistance-based method
2110  IF (ra /= 0) THEN
2111  qh_resist = avdens*avcp*(tsurf - temp_c)/ra
2112  ELSE
2113  qh_resist = nan
2114  ENDIF
2115 
2116  ! choose output QH
2117  SELECT CASE (qhmethod)
2118  CASE (1)
2119  qh = qh_residual
2120  CASE (2)
2121  qh = qh_resist
2122  END SELECT
2123 
2124  END SUBROUTINE suews_cal_qh
2125  !========================================================================
2126 
2127  !===============Resistance Calculations=======================
2128  SUBROUTINE suews_cal_resistance( &
2129  StabilityMethod, &!input:
2130  Diagnose, AerodynamicResistanceMethod, RoughLenHeatMethod, snowUse, &
2131  id, it, gsModel, SMDMethod, &
2132  avdens, avcp, QH_init, zzd, z0m, zdm, &
2133  avU1, Temp_C, VegFraction, &
2134  avkdn, Kmax, G1, G2, G3, G4, G5, G6, S1, S2, TH, TL, dq, &
2135  xsmd, vsmd, MaxConductance, LAIMax, LAI_id, SnowFrac, sfr, &
2136  UStar, TStar, L_mod, &!output
2137  zL, gsc, ResistSurf, RA, RAsnow, rb)
2139  IMPLICIT NONE
2140 
2141  INTEGER, INTENT(in)::StabilityMethod
2142  INTEGER, INTENT(in)::Diagnose
2143  INTEGER, INTENT(in)::AerodynamicResistanceMethod
2144  INTEGER, INTENT(in)::RoughLenHeatMethod
2145  INTEGER, INTENT(in)::snowUse
2146  INTEGER, INTENT(in)::id
2147  INTEGER, INTENT(in)::it !time: day of year and hour
2148  INTEGER, INTENT(in)::gsModel !Choice of gs parameterisation (1 = Ja11, 2 = Wa16)
2149  INTEGER, INTENT(in)::SMDMethod!Method of measured soil moisture
2150 
2151  ! REAL(KIND(1d0)), INTENT(in)::qh_obs
2152  REAL(KIND(1d0)), INTENT(in)::avdens
2153  REAL(KIND(1d0)), INTENT(in)::avcp
2154  REAL(KIND(1d0)), INTENT(in)::QH_init
2155  REAL(KIND(1d0)), INTENT(in)::zzd !Active measurement height (meas. height-displac. height)
2156  REAL(KIND(1d0)), INTENT(in)::z0m !Aerodynamic roughness length
2157  REAL(KIND(1d0)), INTENT(in)::zdm !Displacement height
2158  REAL(KIND(1d0)), INTENT(in)::avU1 !Average wind speed
2159  REAL(KIND(1d0)), INTENT(in)::Temp_C !Air temperature
2160  REAL(KIND(1d0)), INTENT(in)::VegFraction!Fraction of vegetation
2161  REAL(KIND(1d0)), INTENT(in)::avkdn !Average downwelling shortwave radiation
2162  REAL(KIND(1d0)), INTENT(in)::Kmax !Annual maximum hourly solar radiation
2163  REAL(KIND(1d0)), INTENT(in)::G1 !Fitted parameters related to surface res. calculations
2164  REAL(KIND(1d0)), INTENT(in)::G2 !Fitted parameters related to surface res. calculations
2165  REAL(KIND(1d0)), INTENT(in)::G3 !Fitted parameters related to surface res. calculations
2166  REAL(KIND(1d0)), INTENT(in)::G4 !Fitted parameters related to surface res. calculations
2167  REAL(KIND(1d0)), INTENT(in)::G5 !Fitted parameters related to surface res. calculations
2168  REAL(KIND(1d0)), INTENT(in)::G6 !Fitted parameters related to surface res. calculations
2169  REAL(KIND(1d0)), INTENT(in)::S1 !Fitted parameters related to surface res. calculations
2170  REAL(KIND(1d0)), INTENT(in)::S2 !Fitted parameters related to surface res. calculations
2171  REAL(KIND(1d0)), INTENT(in)::TH !Maximum temperature limit
2172  REAL(KIND(1d0)), INTENT(in)::TL !Minimum temperature limit
2173  REAL(KIND(1d0)), INTENT(in)::dq !Specific humidity deficit
2174  REAL(KIND(1d0)), INTENT(in)::xsmd !Measured soil moisture deficit
2175  REAL(KIND(1d0)), INTENT(in)::vsmd !Soil moisture deficit for vegetated surfaces only (QUESTION: what about BSoil?)
2176 
2177  REAL(KIND(1d0)), DIMENSION(3), INTENT(in) ::MaxConductance!Max conductance [mm s-1]
2178  REAL(KIND(1d0)), DIMENSION(3), INTENT(in) ::LAIMax !Max LAI [m2 m-2]
2179  REAL(KIND(1d0)), DIMENSION(3), INTENT(in) ::LAI_id !=LAI_id(id-1,:), LAI for each veg surface [m2 m-2]
2180 
2181  REAL(KIND(1d0)), DIMENSION(nsurf), INTENT(in)::SnowFrac !Surface fraction of snow cover
2182  REAL(KIND(1d0)), DIMENSION(nsurf), INTENT(in)::sfr !Surface fractions [-]
2183 
2184  REAL(KIND(1d0)), INTENT(out)::TStar !T*
2185  REAL(KIND(1d0)), INTENT(out) ::UStar !Friction velocity
2186  REAL(KIND(1d0)), INTENT(out) ::zL !
2187  REAL(KIND(1d0)), INTENT(out) ::gsc !Surface Layer Conductance
2188  REAL(KIND(1d0)), INTENT(out) ::ResistSurf!Surface resistance
2189  REAL(KIND(1d0)), INTENT(out) ::RA !Aerodynamic resistance [s m^-1]
2190  REAL(KIND(1d0)), INTENT(out) ::RAsnow !Aerodynamic resistance for snow [s m^-1]
2191  REAL(KIND(1d0)), INTENT(out) ::rb !boundary layer resistance shuttleworth
2192  REAL(KIND(1d0)), INTENT(out) ::L_mod !Obukhov length
2193  REAL(KIND(1d0)) ::gfunc !gdq*gtemp*gs*gq for photosynthesis calculations
2194  ! REAL(KIND(1d0)) ::H_init !Kinematic sensible heat flux [K m s-1] used to calculate friction velocity
2195 
2196  ! Get first estimate of sensible heat flux. Modified by HCW 26 Feb 2015
2197  ! CALL SUEWS_init_QH( &
2198  ! avdens, avcp, QH_init, qn1, dectime, &
2199  ! H_init)
2200 
2201  IF (diagnose == 1) WRITE (*, *) 'Calling STAB_lumps...'
2202  !u* and Obukhov length out
2203  CALL cal_stab( &
2204  stabilitymethod, & ! input
2205  zzd, & !Active measurement height (meas. height-displac. height)
2206  z0m, & !Aerodynamic roughness length
2207  zdm, & !zero-plane displacement
2208  avu1, & !Average wind speed
2209  temp_c, & !Air temperature
2210  qh_init, & !sensible heat flux
2211  avdens, & ! air density
2212  avcp, & ! heat capacity of air
2213  l_mod, &! output: !Obukhov length
2214  tstar, & !T*, temperature scale
2215  ustar, & !Friction velocity
2216  zl)!Stability scale
2217 
2218  IF (diagnose == 1) WRITE (*, *) 'Calling AerodynamicResistance...'
2219  CALL aerodynamicresistance( &
2220  zzd, &! input:
2221  z0m, &
2222  avu1, &
2223  l_mod, &
2224  ustar, &
2225  vegfraction, &
2226  aerodynamicresistancemethod, &
2227  stabilitymethod, &
2228  roughlenheatmethod, &
2229  ra) ! output:
2230 
2231  IF (snowuse == 1) THEN
2232  IF (diagnose == 1) WRITE (*, *) 'Calling AerodynamicResistance for snow...'
2233  CALL aerodynamicresistance( &
2234  zzd, &! input:
2235  z0m, &
2236  avu1, &
2237  l_mod, &
2238  ustar, &
2239  vegfraction, &
2240  aerodynamicresistancemethod, &
2241  stabilitymethod, &
2242  3, &
2243  rasnow) ! output:
2244  ENDIF
2245 
2246  IF (diagnose == 1) WRITE (*, *) 'Calling SurfaceResistance...'
2247  ! CALL SurfaceResistance(id,it) !qsc and surface resistance out
2248  CALL surfaceresistance( &
2249  id, it, &! input:
2250  smdmethod, snowfrac, sfr, avkdn, temp_c, dq, xsmd, vsmd, maxconductance, &
2251  laimax, lai_id, gsmodel, kmax, &
2252  g1, g2, g3, g4, g5, g6, th, tl, s1, s2, &
2253  gfunc, gsc, resistsurf)! output:
2254 
2255  IF (diagnose == 1) WRITE (*, *) 'Calling BoundaryLayerResistance...'
2256  CALL boundarylayerresistance( &
2257  zzd, &! input: !Active measurement height (meas. height- zero-plane displacement)
2258  z0m, & !Aerodynamic roughness length
2259  avu1, & !Average wind speed
2260  ustar, & ! input/output:
2261  rb) ! output:
2262 
2263  END SUBROUTINE suews_cal_resistance
2264  !========================================================================
2265 
2266  !==============Update output arrays=========================
2267  SUBROUTINE suews_update_outputline( &
2268  AdditionalWater, alb, avkdn, avU10_ms, azimuth, &!input
2269  chSnow_per_interval, dectime, &
2270  drain_per_tstep, E_mod, ev_per_tstep, ext_wu, Fc, Fc_build, fcld, &
2271  Fc_metab, Fc_photo, Fc_respi, Fc_point, Fc_traff, FlowChange, &
2272  h_mod, id, imin, int_wu, it, iy, &
2273  kup, LAI_id, ldown, l_mod, lup, mwh, &
2274  MwStore, &
2275  nsh_real, NWstate_per_tstep, Precip, q2_gkg, &
2276  qeOut, qf, qh, qh_resist, Qm, QmFreez, &
2277  QmRain, qn1, qn1_S, qn1_snowfree, qs, RA, &
2278  resistsurf, RH2, runoffAGimpervious, runoffAGveg, &
2279  runoff_per_tstep, runoffPipes, runoffSoil_per_tstep, &
2280  runoffWaterBody, sfr, smd, smd_nsurf, SnowAlb, SnowRemoval, &
2281  state_id, state_per_tstep, surf_chang_per_tstep, swe, t2_C, tskin_C, &
2282  tot_chang_per_tstep, tsurf, UStar, wu_DecTr, &
2283  wu_EveTr, wu_Grass, z0m, zdm, zenith_deg, &
2284  datetimeLine, dataOutLineSUEWS)!output
2285  IMPLICIT NONE
2286 
2287  REAL(KIND(1d0)), PARAMETER :: NAN = -999
2288  INTEGER, INTENT(in) :: iy
2289  ! INTEGER,INTENT(in) :: iy_prev_t
2290  INTEGER, INTENT(in) :: id
2291  ! INTEGER,INTENT(in) :: id_prev_t
2292  INTEGER, INTENT(in) :: it
2293  INTEGER, INTENT(in) :: imin
2294  ! INTEGER, INTENT(in) :: Gridiv
2295  REAL(KIND(1d0)), INTENT(in) :: AdditionalWater
2296  REAL(KIND(1d0)), INTENT(in) :: alb(nsurf)
2297  REAL(KIND(1d0)), INTENT(in) :: avkdn
2298  REAL(KIND(1d0)), INTENT(in) :: avU10_ms
2299  REAL(KIND(1d0)), INTENT(in) :: azimuth
2300  REAL(KIND(1d0)), INTENT(in) :: chSnow_per_interval
2301  REAL(KIND(1d0)), INTENT(in) :: dectime
2302  REAL(KIND(1d0)), INTENT(in) :: drain_per_tstep
2303  REAL(KIND(1d0)), INTENT(in) :: E_mod
2304  REAL(KIND(1d0)), INTENT(in) :: ev_per_tstep
2305  REAL(KIND(1d0)), INTENT(in) :: ext_wu
2306  REAL(KIND(1d0)), INTENT(in) :: Fc
2307  REAL(KIND(1d0)), INTENT(in) :: Fc_build
2308  REAL(KIND(1d0)), INTENT(in) :: Fc_metab
2309  REAL(KIND(1d0)), INTENT(in) :: Fc_photo
2310  REAL(KIND(1d0)), INTENT(in) :: Fc_respi
2311  REAL(KIND(1d0)), INTENT(in) :: Fc_point
2312  REAL(KIND(1d0)), INTENT(in) :: Fc_traff
2313  REAL(KIND(1d0)), INTENT(in) :: fcld
2314  REAL(KIND(1d0)), INTENT(in) :: FlowChange
2315  REAL(KIND(1d0)), INTENT(in) :: h_mod
2316  REAL(KIND(1d0)), INTENT(in) :: int_wu
2317  REAL(KIND(1d0)), INTENT(in) :: kup
2318  REAL(KIND(1d0)), INTENT(in) :: l_mod
2319  REAL(KIND(1d0)), INTENT(in) :: LAI_id(nvegsurf)
2320  REAL(KIND(1d0)), INTENT(in) :: ldown
2321  REAL(KIND(1d0)), INTENT(in) :: lup
2322  REAL(KIND(1d0)), INTENT(in) :: mwh
2323  REAL(KIND(1d0)), INTENT(in) :: MwStore
2324  REAL(KIND(1d0)), INTENT(in) :: nsh_real
2325  REAL(KIND(1d0)), INTENT(in) :: NWstate_per_tstep
2326  REAL(KIND(1d0)), INTENT(in) :: Precip
2327  REAL(KIND(1d0)), INTENT(in) :: q2_gkg
2328  REAL(KIND(1d0)), INTENT(in) :: qeOut
2329  REAL(KIND(1d0)), INTENT(in) :: qf
2330  REAL(KIND(1d0)), INTENT(in) :: qh
2331  REAL(KIND(1d0)), INTENT(in) :: qh_resist
2332  REAL(KIND(1d0)), INTENT(in) :: Qm
2333  REAL(KIND(1d0)), INTENT(in) :: QmFreez
2334  REAL(KIND(1d0)), INTENT(in) :: QmRain
2335  REAL(KIND(1d0)), INTENT(in) :: qn1
2336  REAL(KIND(1d0)), INTENT(in) :: qn1_S
2337  REAL(KIND(1d0)), INTENT(in) :: qn1_snowfree
2338  REAL(KIND(1d0)), INTENT(in) :: qs
2339  REAL(KIND(1d0)), INTENT(in) :: RA
2340  REAL(KIND(1d0)), INTENT(in) :: resistsurf
2341  REAL(KIND(1d0)), INTENT(in) :: RH2
2342  REAL(KIND(1d0)), INTENT(in) :: runoff_per_tstep
2343  REAL(KIND(1d0)), INTENT(in) :: runoffAGimpervious
2344  REAL(KIND(1d0)), INTENT(in) :: runoffAGveg
2345  REAL(KIND(1d0)), INTENT(in) :: runoffPipes
2346  REAL(KIND(1d0)), INTENT(in) :: runoffSoil_per_tstep
2347  REAL(KIND(1d0)), INTENT(in) :: runoffWaterBody
2348  REAL(KIND(1d0)), INTENT(in) :: sfr(nsurf)
2349  REAL(KIND(1d0)), INTENT(in) :: smd
2350  REAL(KIND(1d0)), INTENT(in) :: smd_nsurf(nsurf)
2351  REAL(KIND(1d0)), INTENT(in) :: SnowAlb
2352  REAL(KIND(1d0)), INTENT(in) :: SnowRemoval(2)
2353  REAL(KIND(1d0)), INTENT(in) :: state_id(nsurf)
2354  REAL(KIND(1d0)), INTENT(in) :: state_per_tstep
2355  REAL(KIND(1d0)), INTENT(in) :: surf_chang_per_tstep
2356  REAL(KIND(1d0)), INTENT(in) :: swe
2357  REAL(KIND(1d0)), INTENT(in) :: t2_C
2358  REAL(KIND(1d0)), INTENT(in) :: tskin_C
2359  REAL(KIND(1d0)), INTENT(in) :: tot_chang_per_tstep
2360  REAL(KIND(1d0)), INTENT(in) :: tsurf
2361  REAL(KIND(1d0)), INTENT(in) :: UStar
2362  REAL(KIND(1d0)), INTENT(in) :: wu_DecTr
2363  REAL(KIND(1d0)), INTENT(in) :: wu_EveTr
2364  REAL(KIND(1d0)), INTENT(in) :: wu_Grass
2365  REAL(KIND(1d0)), INTENT(in) :: z0m
2366  REAL(KIND(1d0)), INTENT(in) :: zdm
2367  REAL(KIND(1d0)), INTENT(in) :: zenith_deg
2368 
2369  REAL(KIND(1D0)), DIMENSION(5), INTENT(OUT)::datetimeLine
2370  REAL(KIND(1d0)), DIMENSION(ncolumnsDataOutSUEWS - 5), INTENT(out) :: dataOutLineSUEWS
2371  ! REAL(KIND(1d0)),DIMENSION(ncolumnsDataOutSnow-5),INTENT(out) :: dataOutLineSnow
2372  ! REAL(KIND(1d0)),DIMENSION(ncolumnsDataOutESTM-5),INTENT(out) :: dataOutLineESTM
2373  ! INTEGER:: is
2374  REAL(KIND(1d0)):: LAI_wt
2375  REAL(KIND(1d0)):: RH2_pct ! RH2 in percentage
2376 
2377  ! the variables below with '_x' endings stand for 'exported' values
2378  REAL(KIND(1d0))::ResistSurf_x
2379  REAL(KIND(1d0))::l_mod_x
2380  REAL(KIND(1d0))::bulkalbedo
2381  REAL(KIND(1d0))::smd_nsurf_x(nsurf)
2382  REAL(KIND(1d0))::state_x(nsurf)
2383 
2384  !=====================================================================
2385  !====================== Prepare data for output ======================
2386  ! values outside of reasonable range are set as NAN-like numbers. TS 10 Jun 2018
2387 
2388  ! Remove non-existing surface type from surface and soil outputs ! Added back in with NANs by HCW 24 Aug 2016
2389  state_x = unpack(spread(nan, dim=1, ncopies=SIZE(sfr)), mask=(sfr < 0.00001), field=state_id)
2390  smd_nsurf_x = unpack(spread(nan, dim=1, ncopies=SIZE(sfr)), mask=(sfr < 0.00001), field=smd_nsurf)
2391 
2392  resistsurf_x = min(9999., resistsurf)
2393 
2394  l_mod_x = max(min(9999., l_mod), -9999.)
2395 
2396  ! Calculate areally-weighted LAI
2397  ! IF(iy == (iy_prev_t +1) .AND. (id-1) == 0) THEN !Check for start of next year and avoid using LAI(id-1) as this is at the start of the year
2398  ! LAI_wt=DOT_PRODUCT(LAI(id_prev_t,:),sfr(1+2:nvegsurf+2))
2399  ! ELSE
2400  ! LAI_wt=DOT_PRODUCT(LAI(id-1,:),sfr(1+2:nvegsurf+2))
2401  ! ENDIF
2402 
2403  lai_wt = dot_product(lai_id(:), sfr(1 + 2:nvegsurf + 2))
2404 
2405  ! Calculate areally-weighted albedo
2406  bulkalbedo = dot_product(alb, sfr)
2407 
2408  ! convert RH2 to a percentage form
2409  rh2_pct = rh2*100.0
2410 
2411  ! NB: this part needs to be reconsidered for calculation logic. TS, 27 Sep 2018
2412  ! TODO: this part should be reconnected to an improved CBL interface. TS 10 Jun 2018
2413  ! uncomment the CBL.HWJ,21 Aug 2019
2414  ! Save qh and qe for CBL in next iteration
2415  ! IF(Qh_choice==1) THEN !use QH and QE from SUEWS
2416  ! qhforCBL(Gridiv) = qh
2417  ! qeforCBL(Gridiv) = qeOut
2418  ! ELSEIF(Qh_choice==2)THEN !use QH and QE from LUMPS
2419  ! qhforCBL(Gridiv) = h_mod
2420  ! qeforCBL(Gridiv) = e_mod
2421  ! ELSEIF(qh_choice==3)THEN !use QH and QE from OBS
2422  ! qhforCBL(Gridiv) = qh_obs
2423  ! qeforCBL(Gridiv) = qe_obs
2424  ! IF(qh_obs<-900.OR.qe_obs<-900)THEN ! observed data has a problem
2425  ! CALL ErrorHint(22,'Unrealistic observed qh or qe_value.',qh_obs,qe_obs,qh_choice)
2426  ! ENDIF
2427  ! ENDIF
2428 
2429  !====================== update output line ==============================
2430  ! date & time:
2431  datetimeline = [ &
2432  REAL(iy, KIND(1D0)), REAL(id, KIND(1D0)), &
2433  REAL(it, KIND(1D0)), REAL(imin, KIND(1D0)), dectime]
2434  !Define the overall output matrix to be printed out step by step
2435  dataoutlinesuews = [ &
2436  avkdn, kup, ldown, lup, tsurf, &
2437  qn1, qf, qs, qh, qeout, &
2438  h_mod, e_mod, qh_resist, &
2439  precip, ext_wu, ev_per_tstep, runoff_per_tstep, tot_chang_per_tstep, &
2440  surf_chang_per_tstep, state_per_tstep, nwstate_per_tstep, drain_per_tstep, smd, &
2441  flowchange/nsh_real, additionalwater, &
2442  runoffsoil_per_tstep, runoffpipes, runoffagimpervious, runoffagveg, runoffwaterbody, &
2443  int_wu, wu_evetr, wu_dectr, wu_grass, &
2444  smd_nsurf_x(1:nsurf - 1), &
2445  state_x(1:nsurf), &
2446  zenith_deg, azimuth, bulkalbedo, fcld, &
2447  lai_wt, z0m, zdm, &
2448  ustar, l_mod, ra, resistsurf, &
2449  fc, &
2450  fc_photo, fc_respi, fc_metab, fc_traff, fc_build, fc_point, &
2451  qn1_snowfree, qn1_s, snowalb, &
2452  qm, qmfreez, qmrain, swe, mwh, mwstore, chsnow_per_interval, &
2453  snowremoval(1:2), &
2454  tskin_c, t2_c, q2_gkg, avu10_ms, rh2_pct & ! surface-level diagonostics
2455  ]
2456  ! set invalid values to NAN
2457  ! dataOutLineSUEWS = set_nan(dataOutLineSUEWS)
2458 
2459  !====================update output line end==============================
2460 
2461  END SUBROUTINE suews_update_outputline
2462  !========================================================================
2463 
2464  !==============Update output arrays=========================
2465  SUBROUTINE suews_update_output( &
2466  SnowUse, storageheatmethod, &!input
2467  ReadLinesMetdata, NumberOfGrids, &
2468  ir, gridiv, &
2469  datetimeLine, dataOutLineSUEWS, dataOutLineSnow, dataOutLineESTM, dataoutLineURSL, &!input
2470  dataOutSUEWS, dataOutSnow, dataOutESTM, dataOutRSL)!inout
2471  IMPLICIT NONE
2472 
2473  INTEGER, INTENT(in) ::ReadLinesMetdata
2474  INTEGER, INTENT(in) ::NumberOfGrids
2475  INTEGER, INTENT(in) ::Gridiv
2476  INTEGER, INTENT(in) ::SnowUse
2477  INTEGER, INTENT(in) ::storageheatmethod
2478  INTEGER, INTENT(in) ::ir
2479 
2480  REAL(KIND(1d0)), DIMENSION(5), INTENT(in) :: datetimeLine
2481  REAL(KIND(1d0)), DIMENSION(ncolumnsDataOutSUEWS - 5), INTENT(in) :: dataOutLineSUEWS
2482  REAL(KIND(1d0)), DIMENSION(ncolumnsDataOutESTM - 5), INTENT(in) :: dataOutLineESTM
2483  REAL(KIND(1d0)), DIMENSION(ncolumnsDataOutSnow - 5), INTENT(in) :: dataOutLineSnow
2484  REAL(KIND(1d0)), DIMENSION(ncolumnsDataOutRSL - 5), INTENT(in) :: dataoutLineURSL
2485 
2486  REAL(KIND(1d0)), INTENT(inout) :: dataOutSUEWS(readlinesmetdata, ncolumnsdataoutsuews, numberofgrids)
2487  REAL(KIND(1d0)), INTENT(inout) :: dataOutSnow(readlinesmetdata, ncolumnsdataoutsnow, numberofgrids)
2488  REAL(KIND(1d0)), INTENT(inout) :: dataOutESTM(readlinesmetdata, ncolumnsdataoutestm, numberofgrids)
2489  REAL(KIND(1d0)), INTENT(inout) :: dataOutRSL(readlinesmetdata, ncolumnsdataoutrsl, numberofgrids)
2490 
2491  !====================== update output arrays ==============================
2492  !Define the overall output matrix to be printed out step by step
2493  dataoutsuews(ir, 1:ncolumnsdataoutsuews, gridiv) = [datetimeline, set_nan(dataoutlinesuews)]
2494  dataoutrsl(ir, 1:ncolumnsdataoutrsl, gridiv) = [datetimeline, set_nan(dataoutlineursl)]
2495  ! ! set invalid values to NAN
2496  ! dataOutSUEWS(ir,6:ncolumnsDataOutSUEWS,Gridiv)=set_nan(dataOutSUEWS(ir,6:ncolumnsDataOutSUEWS,Gridiv))
2497 
2498  IF (snowuse == 1) THEN
2499  dataoutsnow(ir, 1:ncolumnsdataoutsnow, gridiv) = [datetimeline, set_nan(dataoutlinesnow)]
2500  END IF
2501 
2502  IF (storageheatmethod == 4) THEN
2503  dataoutestm(ir, 1:ncolumnsdataoutestm, gridiv) = [datetimeline, set_nan(dataoutlineestm)]
2504  END IF
2505 
2506  !====================update output arrays end==============================
2507 
2508  END SUBROUTINE suews_update_output
2509 
2510  ! !========================================================================
2511  ! SUBROUTINE SUEWS_cal_Diagnostics( &
2512  ! dectime, &!input
2513  ! avU1, Temp_C, avRH, Press_hPa, &
2514  ! qh, qe, &
2515  ! VegFraction, zMeas, z0m, zdm, RA, avdens, avcp, lv_J_kg, tstep_real, &
2516  ! RoughLenHeatMethod, StabilityMethod, &
2517  ! avU10_ms, t2_C, q2_gkg, tskin_C, RH2)!output
2518  ! ! TS 03 Aug 2018: added limit on q2 by restricting RH2_max to 100%
2519  ! ! TS 31 Jul 2018: removed dependence on surface variables (Tsurf, qsat)
2520  ! ! TS 26 Jul 2018: improved the calculation logic
2521  ! ! TS 05 Sep 2017: improved interface
2522  ! ! TS 20 May 2017: calculate surface-level diagonostics
2523  ! IMPLICIT NONE
2524  ! REAL(KIND(1d0)), INTENT(in) ::dectime
2525  ! REAL(KIND(1d0)), INTENT(in) ::avU1, Temp_C, avRH
2526  ! REAL(KIND(1d0)), INTENT(in) ::qh
2527  ! REAL(KIND(1d0)), INTENT(in) ::Press_hPa, qe
2528  ! REAL(KIND(1d0)), INTENT(in) :: VegFraction, z0m, RA, avdens, avcp, lv_J_kg, tstep_real
2529  ! REAL(KIND(1d0)), INTENT(in) :: zMeas! height for measurement
2530  ! REAL(KIND(1d0)), INTENT(in) :: zdm ! displacement height
2531 
2532  ! ! INTEGER,INTENT(in) :: opt ! 0 for momentum, 1 for temperature, 2 for humidity
2533  ! INTEGER, INTENT(in) :: RoughLenHeatMethod, StabilityMethod
2534 
2535  ! REAL(KIND(1d0)), INTENT(out):: avU10_ms, t2_C, q2_gkg, tskin_C, RH2
2536  ! REAL(KIND(1d0))::qa_gkg
2537  ! REAL(KIND(1d0)), PARAMETER::k = 0.4
2538 
2539  ! ! wind speed:
2540  ! CALL diagSfc( &
2541  ! 0, &
2542  ! zMeas, avU1, 0d0, 10d0, avU10_ms, &
2543  ! VegFraction, &
2544  ! z0m, zdm, avdens, avcp, lv_J_kg, &
2545  ! avU1, Temp_C, qh, &
2546  ! RoughLenHeatMethod, StabilityMethod, tstep_real, dectime)
2547 
2548  ! ! temperature at 2 m agl:
2549  ! CALL diagSfc( &
2550  ! 1, &
2551  ! zMeas, Temp_C, qh, 2d0, t2_C, &
2552  ! VegFraction, &
2553  ! z0m, zdm, avdens, avcp, lv_J_kg, &
2554  ! avU1, Temp_C, qh, &
2555  ! RoughLenHeatMethod, StabilityMethod, tstep_real, dectime)
2556 
2557  ! ! skin temperature:
2558  ! tskin_C = qh/(avdens*avcp)*RA + temp_C
2559 
2560  ! ! humidity:
2561  ! qa_gkg = RH2qa(avRH/100, Press_hPa, Temp_c)
2562  ! CALL diagSfc( &
2563  ! 2, &
2564  ! zMeas, qa_gkg, qe, 2d0, q2_gkg, &
2565  ! VegFraction, &
2566  ! z0m, zdm, avdens, avcp, lv_J_kg, &
2567  ! avU1, Temp_C, qh, &
2568  ! RoughLenHeatMethod, StabilityMethod, tstep_real, dectime)
2569  ! ! re-examine if the diagnostic RH2 > 100% ?
2570  ! RH2 = qa2RH(q2_gkg, Press_hPa, Temp_c)
2571  ! IF (RH2 > 1) THEN
2572  ! ! if so, limit RH2 to 100%
2573  ! RH2 = 1d0
2574  ! ! and adjust the diagnostic q2_gkg
2575  ! q2_gkg = RH2qa(RH2, Press_hPa, Temp_c)
2576  ! END IF
2577 
2578  ! END SUBROUTINE SUEWS_cal_Diagnostics
2579 
2580  ! calculate several surface fraction related parameters
2581  SUBROUTINE suews_cal_surf( &
2582  sfr, & !input
2583  vegfraction, ImpervFraction, PervFraction, NonWaterFraction) ! output
2584  IMPLICIT NONE
2585 
2586  REAL(KIND(1D0)), DIMENSION(NSURF), INTENT(IN)::sfr
2587  REAL(KIND(1D0)), INTENT(OUT)::VegFraction
2588  REAL(KIND(1D0)), INTENT(OUT)::ImpervFraction
2589  REAL(KIND(1D0)), INTENT(OUT)::PervFraction
2590  REAL(KIND(1D0)), INTENT(OUT)::NonWaterFraction
2591 
2592  vegfraction = sfr(conifsurf) + sfr(decidsurf) + sfr(grasssurf)
2593  impervfraction = sfr(pavsurf) + sfr(bldgsurf)
2594  pervfraction = 1 - impervfraction
2595  nonwaterfraction = 1 - sfr(watersurf)
2596 
2597  END SUBROUTINE suews_cal_surf
2598 
2599  ! SUBROUTINE diagSfc( &
2600  ! opt, &
2601  ! zMeas, xMeas, xFlux, zDiag, xDiag, &
2602  ! VegFraction, &
2603  ! z0m, zd, avdens, avcp, lv_J_kg, &
2604  ! avU1, Temp_C, qh, &
2605  ! RoughLenHeatMethod, StabilityMethod, tstep_real, dectime)
2606  ! ! TS 31 Jul 2018: removed dependence on surface variables (Tsurf, qsat)
2607  ! ! TS 26 Jul 2018: improved the calculation logic
2608  ! ! TS 05 Sep 2017: improved interface
2609  ! ! TS 20 May 2017: calculate surface-level diagonostics
2610 
2611  ! IMPLICIT NONE
2612  ! REAL(KIND(1d0)), INTENT(in) :: dectime
2613  ! REAL(KIND(1d0)), INTENT(in) :: qh ! sensible heat flux
2614  ! REAL(KIND(1d0)), INTENT(in) :: z0m, avdens, avcp, lv_J_kg, tstep_real
2615  ! REAL(KIND(1d0)), INTENT(in) :: avU1, Temp_C ! atmospheric level variables
2616  ! REAL(KIND(1d0)), INTENT(in) :: zDiag ! height for diagonostics
2617  ! REAL(KIND(1d0)), INTENT(in) :: zMeas! height for measurement
2618  ! REAL(KIND(1d0)), INTENT(in) :: zd ! displacement height
2619  ! REAL(KIND(1d0)), INTENT(in) :: xMeas ! measurement at height
2620  ! REAL(KIND(1d0)), INTENT(in) :: xFlux!
2621  ! REAL(KIND(1d0)), INTENT(in) :: VegFraction ! vegetation fraction
2622 
2623  ! INTEGER, INTENT(in) :: opt ! 0 for momentum, 1 for temperature, 2 for humidity
2624  ! INTEGER, INTENT(in) :: RoughLenHeatMethod, StabilityMethod
2625 
2626  ! REAL(KIND(1d0)), INTENT(out):: xDiag
2627 
2628  ! REAL(KIND(1d0)) :: L_mod
2629  ! REAL(KIND(1d0)) :: psimz0, psihzDiag, psihzMeas, psihz0, psimzDiag ! stability correction functions
2630  ! REAL(KIND(1d0)) :: z0h ! Roughness length for heat
2631  ! REAL(KIND(1d0)) :: zDiagzd! height for diagnositcs
2632  ! REAL(KIND(1d0)) :: zMeaszd
2633  ! REAL(KIND(1d0)) :: tlv, H_kms, TStar, zL, UStar
2634  ! REAL(KIND(1d0)), PARAMETER :: muu = 1.46e-5 !molecular viscosity
2635  ! REAL(KIND(1d0)), PARAMETER :: nan = -999
2636  ! REAL(KIND(1d0)), PARAMETER :: zdm = 0 ! assuming Displacement height is ZERO
2637  ! REAL(KIND(1d0)), PARAMETER::k = 0.4
2638 
2639  ! tlv = lv_J_kg/tstep_real !Latent heat of vapourisation per timestep
2640  ! zDiagzd = zDiag + z0m ! height at hgtX assuming Displacement height is ZERO; set lower limit as z0 to prevent arithmetic error, zd=0
2641 
2642  ! ! get !Kinematic sensible heat flux [K m s-1] used to calculate friction velocity
2643  ! CALL SUEWS_init_QH( &
2644  ! avdens, avcp, qh, 0d0, dectime, & ! use qh as qh_obs to initialise H_init
2645  ! H_kms)!output
2646 
2647  ! ! redo the calculation for stability correction
2648  ! CALL cal_Stab( &
2649  ! ! input
2650  ! StabilityMethod, &
2651  ! dectime, & !Decimal time
2652  ! zDiagzd, & !Active measurement height (meas. height-displac. height)
2653  ! z0m, & !Aerodynamic roughness length
2654  ! zdm, & !Displacement height
2655  ! avU1, & !Average wind speed
2656  ! Temp_C, & !Air temperature
2657  ! H_kms, & !Kinematic sensible heat flux [K m s-1] used to calculate friction velocity
2658  ! ! output:
2659  ! L_MOD, & !Obukhov length
2660  ! TStar, & !T*
2661  ! UStar, & !Friction velocity
2662  ! zL)!Stability scale
2663 
2664  ! !***************************************************************
2665  ! ! log-law based stability corrections:
2666  ! ! Roughness length for heat
2667  ! z0h = cal_z0V(RoughLenHeatMethod, z0m, VegFraction, UStar)
2668 
2669  ! ! stability correction functions
2670  ! ! momentum:
2671  ! psimzDiag = stab_psi_mom(StabilityMethod, zDiagzd/L_mod)
2672  ! ! psimz2=stab_fn_mom(StabilityMethod,z2zd/L_mod,z2zd/L_mod)
2673  ! psimz0 = stab_psi_mom(StabilityMethod, z0m/L_mod)
2674 
2675  ! ! heat and vapor: assuming both are the same
2676  ! ! psihz2=stab_fn_heat(StabilityMethod,z2zd/L_mod,z2zd/L_mod)
2677  ! psihz0 = stab_psi_heat(StabilityMethod, z0h/L_mod)
2678 
2679  ! !***************************************************************
2680  ! SELECT CASE (opt)
2681  ! CASE (0) ! wind (momentum) at hgtX=10 m
2682  ! zDiagzd = zDiag + z0m! set lower limit as z0h to prevent arithmetic error, zd=0
2683 
2684  ! ! stability correction functions
2685  ! ! momentum:
2686  ! psimzDiag = stab_psi_mom(StabilityMethod, zDiagzd/L_mod)
2687  ! psimz0 = stab_psi_mom(StabilityMethod, z0m/L_mod)
2688  ! xDiag = UStar/k*(LOG(zDiagzd/z0m) - psimzDiag + psimz0) ! Brutsaert (2005), p51, eq.2.54
2689 
2690  ! CASE (1) ! temperature at hgtX=2 m
2691  ! zMeaszd = zMeas - zd
2692  ! zDiagzd = zDiag + z0h! set lower limit as z0h to prevent arithmetic error, zd=0
2693 
2694  ! ! heat and vapor: assuming both are the same
2695  ! psihzMeas = stab_psi_heat(StabilityMethod, zMeaszd/L_mod)
2696  ! psihzDiag = stab_psi_heat(StabilityMethod, zDiagzd/L_mod)
2697  ! ! psihz0=stab_fn_heat(StabilityMethod,z0h/L_mod,z0h/L_mod)
2698  ! xDiag = xMeas + xFlux/(k*UStar*avdens*avcp)*(LOG(zMeaszd/zDiagzd) - (psihzMeas - psihzDiag)) ! Brutsaert (2005), p51, eq.2.55
2699  ! ! IF ( ABS((LOG(z2zd/z0h)-psihz2+psihz0))>10 ) THEN
2700  ! ! PRINT*, '#####################################'
2701  ! ! PRINT*, 'xSurf',xSurf
2702  ! ! PRINT*, 'xFlux',xFlux
2703  ! ! PRINT*, 'k*us*avdens*avcp',k*us*avdens*avcp
2704  ! ! PRINT*, 'k',k
2705  ! ! PRINT*, 'us',us
2706  ! ! PRINT*, 'avdens',avdens
2707  ! ! PRINT*, 'avcp',avcp
2708  ! ! PRINT*, 'xFlux/X',xFlux/(k*us*avdens*avcp)
2709  ! ! PRINT*, 'stab',(LOG(z2zd/z0h)-psihz2+psihz0)
2710  ! ! PRINT*, 'LOG(z2zd/z0h)',LOG(z2zd/z0h)
2711  ! ! PRINT*, 'z2zd',z2zd,'L_mod',L_mod,'z0h',z0h
2712  ! ! PRINT*, 'z2zd/L_mod',z2zd/L_mod
2713  ! ! PRINT*, 'psihz2',psihz2
2714  ! ! PRINT*, 'psihz0',psihz0
2715  ! ! PRINT*, 'psihz2-psihz0',psihz2-psihz0
2716  ! ! PRINT*, 'xDiag',xDiag
2717  ! ! PRINT*, '*************************************'
2718  ! ! END IF
2719 
2720  ! CASE (2) ! humidity at hgtX=2 m
2721  ! zMeaszd = zMeas - zd
2722  ! zDiagzd = zDiag + z0h! set lower limit as z0h to prevent arithmetic error, zd=0
2723 
2724  ! ! heat and vapor: assuming both are the same
2725  ! psihzMeas = stab_psi_heat(StabilityMethod, zMeaszd/L_mod)
2726  ! psihzDiag = stab_psi_heat(StabilityMethod, zDiagzd/L_mod)
2727  ! ! psihz0=stab_fn_heat(StabilityMethod,z0h/L_mod,z0h/L_mod)
2728 
2729  ! xDiag = xMeas + xFlux/(k*UStar*avdens*tlv)*(LOG(zMeaszd/zDiagzd) - (psihzMeas - psihzDiag)) ! Brutsaert (2005), p51, eq.2.56
2730 
2731  ! END SELECT
2732 
2733  ! END SUBROUTINE diagSfc
2734 
2735  !===============set variable of invalid value to NAN=====================
2736  ELEMENTAL FUNCTION set_nan(x) RESULT(xx)
2737  IMPLICIT NONE
2738  REAL(KIND(1d0)), PARAMETER::pNAN = 30000 ! 30000 to prevent water_state being filtered out as it can be large
2739  REAL(KIND(1d0)), PARAMETER::NAN = -999
2740  REAL(KIND(1d0)), INTENT(in)::x
2741  REAL(KIND(1d0))::xx
2742 
2743  IF (abs(x) > pnan) THEN
2744  xx = nan
2745  ELSE
2746  xx = x
2747  ENDIF
2748 
2749  END FUNCTION set_nan
2750  !========================================================================
2751 
2752  !===============the functions below are only for test in f2py conversion===
2753  FUNCTION square(x) RESULT(xx)
2754  IMPLICIT NONE
2755  REAL(KIND(1d0)), PARAMETER::pNAN = 9999
2756  REAL(KIND(1d0)), PARAMETER::NAN = -999
2757  REAL(KIND(1d0)), INTENT(in)::x
2758  REAL(KIND(1d0))::xx
2759 
2760  xx = x**2 + nan/pnan
2761  xx = x**2
2762 
2763  END FUNCTION square
2764 
2765  FUNCTION square_real(x) RESULT(xx)
2766  IMPLICIT NONE
2767  REAL, PARAMETER::pNAN = 9999
2768  REAL, PARAMETER::NAN = -999
2769  REAL, INTENT(in)::x
2770  REAL::xx
2771 
2772  xx = x**2 + nan/pnan
2773  xx = x**2
2774 
2775  END FUNCTION square_real
2776 
2777  SUBROUTINE output_name_n(i, name, group, aggreg)
2778  ! used by f2py module `SuPy` to handle output names
2779  IMPLICIT NONE
2780  ! the dimension is potentially incorrect,
2781  ! which should be consistent with that in output module
2782  INTEGER, INTENT(in) :: i
2783  CHARACTER(len=15), INTENT(out) :: name, group, aggreg
2784 
2785  INTEGER :: nVar
2786  nvar = SIZE(varlistall, dim=1)
2787  IF (i < nvar .AND. i > 0) THEN
2788  name = trim(varlistall(i)%header)
2789  group = trim(varlistall(i)%group)
2790  aggreg = trim(varlistall(i)%aggreg)
2791  ELSE
2792  name = ''
2793  group = ''
2794  aggreg = ''
2795  END IF
2796 
2797  END SUBROUTINE output_name_n
2798 
2799  SUBROUTINE output_size(nVar)
2800  ! used by f2py module `SuPy` to get size of the output list
2801  IMPLICIT NONE
2802  ! the dimension is potentially incorrect,
2803  ! which should be consistent with that in output module
2804  INTEGER, INTENT(out) :: nVar
2805 
2806  nvar = SIZE(varlistall, dim=1)
2807 
2808  END SUBROUTINE output_size
2809 
2810  SUBROUTINE suews_cal_multitsteps( &
2811  MetForcingBlock, len_sim, &
2812  AerodynamicResistanceMethod, AH_MIN, AHProf_24hr, AH_SLOPE_Cooling, & ! input&inout in alphabetical order
2813  AH_SLOPE_Heating, &
2814  alb, AlbMax_DecTr, AlbMax_EveTr, AlbMax_Grass, &
2815  AlbMin_DecTr, AlbMin_EveTr, AlbMin_Grass, &
2816  alpha_bioCO2, alpha_enh_bioCO2, alt, BaseT, BaseTe, &
2817  BaseTHDD, beta_bioCO2, beta_enh_bioCO2, bldgH, CapMax_dec, CapMin_dec, &
2818  chAnOHM, CO2PointSource, cpAnOHM, CRWmax, CRWmin, DayWat, DayWatPer, &
2819  DecTreeH, Diagnose, DiagQN, DiagQS, DRAINRT, &
2820  dt_since_start, dqndt, qn1_av, dqnsdt, qn1_s_av, &
2821  EF_umolCO2perJ, emis, EmissionsMethod, EnEF_v_Jkm, endDLS, EveTreeH, FAIBldg, &
2822  FAIDecTree, FAIEveTree, Faut, FcEF_v_kgkm, FlowChange, &
2823  FrFossilFuel_Heat, FrFossilFuel_NonHeat, G1, G2, G3, G4, G5, G6, GDD_id, &
2824  GDDFull, Gridiv, gsModel, HDD_id, HumActivity_24hr, &
2825  IceFrac, Ie_a, Ie_end, Ie_m, Ie_start, &
2826  InternalWaterUse_h, IrrFracConif, IrrFracDecid, IrrFracGrass, EvapMethod, &
2827  kkAnOHM, Kmax, LAI_id, LAICalcYes, LAIMax, LAIMin, &
2828  LAIPower, LAIType, lat, lng, MaxConductance, MaxFCMetab, MaxQFMetab, &
2829  SnowWater, MinFCMetab, MinQFMetab, min_res_bioCO2, &
2830  NARP_EMIS_SNOW, NARP_TRANS_SITE, NetRadiationMethod, &
2831  OHM_coef, OHMIncQF, OHM_threshSW, &
2832  OHM_threshWD, PipeCapacity, PopDensDaytime, &
2833  PopDensNighttime, PopProf_24hr, PorMax_dec, PorMin_dec, &
2834  PrecipLimit, PrecipLimitAlb, &
2835  QF0_BEU, Qf_A, Qf_B, Qf_C, &
2836  RadMeltFact, RAINCOVER, RainMaxRes, resp_a, resp_b, &
2837  RoughLenHeatMethod, RoughLenMomMethod, RunoffToWater, S1, S2, &
2838  SatHydraulicConduct, SDDFull, SDD_id, sfr, SMDMethod, SnowAlb, SnowAlbMax, &
2839  SnowAlbMin, SnowPackLimit, SnowDens, SnowDensMax, SnowDensMin, SnowfallCum, SnowFrac, &
2840  SnowLimBldg, SnowLimPaved, SnowPack, SnowProf_24hr, snowUse, SoilDepth, &
2841  soilstore_id, SoilStoreCap, StabilityMethod, startDLS, state_id, StateLimit, &
2842  StorageHeatMethod, StoreDrainPrm, SurfaceArea, Tair_av, tau_a, tau_f, tau_r, &
2843  T_CRITIC_Cooling, T_CRITIC_Heating, TempMeltFact, TH, &
2844  theta_bioCO2, timezone, TL, TrafficRate, TrafficUnits, &
2845  Tmin_id, Tmax_id, lenday_id, &
2846  TraffProf_24hr, Ts5mindata_ir, tstep, tstep_prev, veg_type, &
2847  WaterDist, WaterUseMethod, WetThresh, &
2848  WUDay_id, DecidCap_id, albDecTr_id, albEveTr_id, albGrass_id, porosity_id, &
2849  WUProfA_24hr, WUProfM_24hr, Z, z0m_in, zdm_in, &
2850  dataOutBlockSUEWS, dataOutBlockSnow, dataOutBlockESTM, dataOutBlockRSL, &!output
2851  DailyStateBlock)
2853  IMPLICIT NONE
2854  ! input:
2855  ! met forcing block
2856  REAL(KIND(1D0)), DIMENSION(len_sim, 24), INTENT(IN) ::MetForcingBlock
2857  INTEGER, INTENT(IN) :: len_sim
2858  ! input variables
2859  INTEGER, INTENT(IN)::AerodynamicResistanceMethod
2860  INTEGER, INTENT(IN)::Diagnose
2861  INTEGER, INTENT(IN)::DiagQN
2862  INTEGER, INTENT(IN)::DiagQS
2863  INTEGER, INTENT(IN)::startDLS
2864  INTEGER, INTENT(IN)::endDLS
2865  INTEGER, INTENT(IN)::EmissionsMethod
2866  INTEGER, INTENT(IN)::Gridiv
2867  INTEGER, INTENT(IN)::gsModel
2868  INTEGER, INTENT(IN)::Ie_end
2869  INTEGER, INTENT(IN)::Ie_start
2870  INTEGER, INTENT(IN)::EvapMethod
2871  INTEGER, INTENT(IN)::LAICalcYes
2872  INTEGER, INTENT(IN)::NetRadiationMethod
2873  INTEGER, INTENT(IN)::OHMIncQF
2874  INTEGER, INTENT(IN)::RoughLenHeatMethod
2875  INTEGER, INTENT(IN)::RoughLenMomMethod
2876  INTEGER, INTENT(IN)::SMDMethod
2877  INTEGER, INTENT(IN)::snowUse
2878  INTEGER, INTENT(IN)::StabilityMethod
2879  INTEGER, INTENT(IN)::StorageHeatMethod
2880  INTEGER, INTENT(IN)::tstep
2881  INTEGER, INTENT(IN)::tstep_prev ! tstep size of the previous step
2882  ! dt_since_start is intentionally made as inout to keep naming consistency with the embedded subroutine
2883  INTEGER, INTENT(inout)::dt_since_start ! time since simulation starts [s]
2884  INTEGER, INTENT(IN)::veg_type
2885  INTEGER, INTENT(IN)::WaterUseMethod
2886 
2887  INTEGER, DIMENSION(NVEGSURF), INTENT(IN)::LAIType
2888 
2889  REAL(KIND(1D0)), INTENT(IN)::AlbMax_DecTr
2890  REAL(KIND(1D0)), INTENT(IN)::AlbMax_EveTr
2891  REAL(KIND(1D0)), INTENT(IN)::AlbMax_Grass
2892  REAL(KIND(1D0)), INTENT(IN)::AlbMin_DecTr
2893  REAL(KIND(1D0)), INTENT(IN)::AlbMin_EveTr
2894  REAL(KIND(1D0)), INTENT(IN)::AlbMin_Grass
2895  REAL(KIND(1D0)), INTENT(IN)::alt
2896  ! REAL(KIND(1D0)),INTENT(IN)::avkdn
2897  ! REAL(KIND(1D0)),INTENT(IN)::avRh
2898  ! REAL(KIND(1D0)),INTENT(IN)::avU1
2899  REAL(KIND(1D0)), INTENT(IN)::BaseTHDD
2900  REAL(KIND(1D0)), INTENT(IN)::bldgH
2901  REAL(KIND(1D0)), INTENT(IN)::CapMax_dec
2902  REAL(KIND(1D0)), INTENT(IN)::CapMin_dec
2903  REAL(KIND(1D0)), INTENT(IN)::CO2PointSource
2904  REAL(KIND(1D0)), INTENT(IN)::CRWmax
2905  REAL(KIND(1D0)), INTENT(IN)::CRWmin
2906  REAL(KIND(1D0)), INTENT(IN)::DecTreeH
2907  REAL(KIND(1D0)), INTENT(IN)::DRAINRT
2908  REAL(KIND(1D0)), INTENT(IN)::EF_umolCO2perJ
2909  REAL(KIND(1D0)), INTENT(IN)::EnEF_v_Jkm
2910  REAL(KIND(1D0)), INTENT(IN)::EveTreeH
2911  REAL(KIND(1D0)), INTENT(IN)::FAIBldg
2912  REAL(KIND(1D0)), INTENT(IN)::FAIDecTree
2913  REAL(KIND(1D0)), INTENT(IN)::FAIEveTree
2914  REAL(KIND(1D0)), INTENT(IN)::Faut
2915  ! REAL(KIND(1D0)),INTENT(IN)::fcld_obs
2916  REAL(KIND(1D0)), INTENT(IN)::FlowChange
2917  REAL(KIND(1D0)), INTENT(IN)::FrFossilFuel_Heat
2918  REAL(KIND(1D0)), INTENT(IN)::FrFossilFuel_NonHeat
2919  REAL(KIND(1D0)), INTENT(IN)::G1
2920  REAL(KIND(1D0)), INTENT(IN)::G2
2921  REAL(KIND(1D0)), INTENT(IN)::G3
2922  REAL(KIND(1D0)), INTENT(IN)::G4
2923  REAL(KIND(1D0)), INTENT(IN)::G5
2924  REAL(KIND(1D0)), INTENT(IN)::G6
2925  REAL(KIND(1D0)), INTENT(IN)::InternalWaterUse_h
2926  REAL(KIND(1D0)), INTENT(IN)::IrrFracConif
2927  REAL(KIND(1D0)), INTENT(IN)::IrrFracDecid
2928  REAL(KIND(1D0)), INTENT(IN)::IrrFracGrass
2929  REAL(KIND(1D0)), INTENT(IN)::Kmax
2930  ! REAL(KIND(1D0)),INTENT(IN)::LAI_obs
2931  REAL(KIND(1D0)), INTENT(IN)::lat
2932  ! REAL(KIND(1D0)),INTENT(IN)::ldown_obs
2933  REAL(KIND(1D0)), INTENT(IN)::lng
2934  REAL(KIND(1D0)), INTENT(IN)::MaxFCMetab
2935  REAL(KIND(1D0)), INTENT(IN)::MaxQFMetab
2936  REAL(KIND(1D0)), INTENT(IN)::MinFCMetab
2937  REAL(KIND(1D0)), INTENT(IN)::MinQFMetab
2938  REAL(KIND(1D0)), INTENT(IN)::NARP_EMIS_SNOW
2939  REAL(KIND(1D0)), INTENT(IN)::NARP_TRANS_SITE
2940  REAL(KIND(1D0)), INTENT(IN)::PipeCapacity
2941  REAL(KIND(1D0)), INTENT(IN)::PopDensNighttime
2942  REAL(KIND(1D0)), INTENT(IN)::PorMax_dec
2943  REAL(KIND(1D0)), INTENT(IN)::PorMin_dec
2944  ! REAL(KIND(1D0)),INTENT(IN)::Precip
2945  REAL(KIND(1D0)), INTENT(IN)::PrecipLimit
2946  REAL(KIND(1D0)), INTENT(IN)::PrecipLimitAlb
2947  ! REAL(KIND(1D0)),INTENT(IN)::Press_hPa
2948  ! REAL(KIND(1D0)),INTENT(IN)::qh_obs
2949  ! REAL(KIND(1D0)),INTENT(IN)::qn1_obs
2950  ! REAL(KIND(1D0)),INTENT(IN)::qs_obs
2951  ! REAL(KIND(1D0)),INTENT(IN)::qf_obs
2952  REAL(KIND(1D0)), INTENT(IN)::RadMeltFact
2953  REAL(KIND(1D0)), INTENT(IN)::RAINCOVER
2954  REAL(KIND(1D0)), INTENT(IN)::RainMaxRes
2955  REAL(KIND(1D0)), INTENT(IN)::RunoffToWater
2956  REAL(KIND(1D0)), INTENT(IN)::S1
2957  REAL(KIND(1D0)), INTENT(IN)::S2
2958  REAL(KIND(1D0)), INTENT(IN)::SnowAlbMax
2959  REAL(KIND(1D0)), INTENT(IN)::SnowAlbMin
2960  REAL(KIND(1D0)), INTENT(IN)::SnowDensMax
2961  REAL(KIND(1D0)), INTENT(IN)::SnowDensMin
2962  REAL(KIND(1D0)), INTENT(IN)::SnowLimBldg
2963  REAL(KIND(1D0)), INTENT(IN)::SnowLimPaved
2964  ! REAL(KIND(1D0)),INTENT(IN)::snowFrac_obs
2965  REAL(KIND(1D0)), INTENT(IN)::SurfaceArea
2966  REAL(KIND(1D0)), INTENT(IN)::tau_a
2967  REAL(KIND(1D0)), INTENT(IN)::tau_f
2968  REAL(KIND(1D0)), INTENT(IN)::tau_r
2969  ! REAL(KIND(1D0)),INTENT(IN)::Temp_C
2970  REAL(KIND(1D0)), INTENT(IN)::TempMeltFact
2971  REAL(KIND(1D0)), INTENT(IN)::TH
2972  REAL(KIND(1D0)), INTENT(IN)::timezone
2973  REAL(KIND(1D0)), INTENT(IN)::TL
2974  REAL(KIND(1D0)), INTENT(IN)::TrafficUnits
2975  ! REAL(KIND(1D0)),INTENT(IN)::xsmd
2976  REAL(KIND(1D0)), INTENT(IN)::Z
2977  REAL(KIND(1D0)), INTENT(IN)::z0m_in
2978  REAL(KIND(1D0)), INTENT(IN)::zdm_in
2979 
2980  REAL(KIND(1D0)), DIMENSION(2), INTENT(IN) ::AH_MIN
2981  REAL(KIND(1D0)), DIMENSION(2), INTENT(IN) ::AH_SLOPE_Cooling
2982  REAL(KIND(1D0)), DIMENSION(2), INTENT(IN) ::AH_SLOPE_Heating
2983  REAL(KIND(1D0)), DIMENSION(2), INTENT(IN) ::FcEF_v_kgkm
2984  REAL(KIND(1D0)), DIMENSION(2), INTENT(IN) ::QF0_BEU
2985  REAL(KIND(1D0)), DIMENSION(2), INTENT(IN) ::Qf_A
2986  REAL(KIND(1D0)), DIMENSION(2), INTENT(IN) ::Qf_B
2987  REAL(KIND(1D0)), DIMENSION(2), INTENT(IN) ::Qf_C
2988  ! REAL(KIND(1D0)), DIMENSION(2), INTENT(IN) ::Numcapita
2989  REAL(KIND(1D0)), DIMENSION(2), INTENT(IN) ::PopDensDaytime
2990  REAL(KIND(1D0)), DIMENSION(2), INTENT(IN) ::T_CRITIC_Cooling
2991  REAL(KIND(1D0)), DIMENSION(2), INTENT(IN) ::T_CRITIC_Heating
2992  REAL(KIND(1D0)), DIMENSION(2), INTENT(IN) ::TrafficRate
2993  REAL(KIND(1D0)), DIMENSION(3), INTENT(IN) ::Ie_a
2994  REAL(KIND(1D0)), DIMENSION(3), INTENT(IN) ::Ie_m
2995  REAL(KIND(1D0)), DIMENSION(3), INTENT(IN) ::MaxConductance
2996  REAL(KIND(1D0)), DIMENSION(7), INTENT(IN) ::DayWat
2997  REAL(KIND(1D0)), DIMENSION(7), INTENT(IN) ::DayWatPer
2998  REAL(KIND(1D0)), DIMENSION(nsurf + 1), INTENT(IN)::OHM_threshSW
2999  REAL(KIND(1D0)), DIMENSION(nsurf + 1), INTENT(IN)::OHM_threshWD
3000  REAL(KIND(1D0)), DIMENSION(NSURF), INTENT(IN) ::chAnOHM
3001  REAL(KIND(1D0)), DIMENSION(NSURF), INTENT(IN) ::cpAnOHM
3002  REAL(KIND(1D0)), DIMENSION(NSURF), INTENT(IN) ::emis
3003  REAL(KIND(1D0)), DIMENSION(NSURF), INTENT(IN) ::kkAnOHM
3004  REAL(KIND(1D0)), DIMENSION(NSURF), INTENT(IN) ::SatHydraulicConduct
3005  REAL(KIND(1D0)), DIMENSION(NSURF), INTENT(IN) ::sfr
3006  REAL(KIND(1D0)), DIMENSION(NSURF), INTENT(IN) ::SnowPackLimit
3007  REAL(KIND(1D0)), DIMENSION(NSURF), INTENT(IN) ::SoilDepth
3008  REAL(KIND(1D0)), DIMENSION(NSURF), INTENT(IN) ::SoilStoreCap
3009  REAL(KIND(1D0)), DIMENSION(NSURF), INTENT(IN) ::StateLimit
3010  REAL(KIND(1D0)), DIMENSION(NSURF), INTENT(IN) ::WetThresh
3011  REAL(KIND(1D0)), DIMENSION(NVEGSURF), INTENT(IN) ::alpha_bioCO2
3012  REAL(KIND(1D0)), DIMENSION(NVEGSURF), INTENT(IN) ::alpha_enh_bioCO2
3013  REAL(KIND(1D0)), DIMENSION(NVEGSURF), INTENT(IN) ::BaseT
3014  REAL(KIND(1D0)), DIMENSION(NVEGSURF), INTENT(IN) ::BaseTe
3015  REAL(KIND(1D0)), DIMENSION(NVEGSURF), INTENT(IN) ::beta_bioCO2
3016  REAL(KIND(1D0)), DIMENSION(NVEGSURF), INTENT(IN) ::beta_enh_bioCO2
3017  REAL(KIND(1D0)), DIMENSION(NVEGSURF), INTENT(IN) ::GDDFull
3018  REAL(KIND(1D0)), DIMENSION(NVEGSURF), INTENT(IN) ::LAIMax
3019  REAL(KIND(1D0)), DIMENSION(NVEGSURF), INTENT(IN) ::LAIMin
3020  REAL(KIND(1D0)), DIMENSION(NVEGSURF), INTENT(IN) ::min_res_bioCO2
3021  REAL(KIND(1D0)), DIMENSION(NVEGSURF), INTENT(IN) ::resp_a
3022  REAL(KIND(1D0)), DIMENSION(NVEGSURF), INTENT(IN) ::resp_b
3023  REAL(KIND(1D0)), DIMENSION(NVEGSURF), INTENT(IN) ::SDDFull
3024  REAL(KIND(1D0)), DIMENSION(0:23, 2), INTENT(IN) ::SnowProf_24hr
3025  REAL(KIND(1D0)), DIMENSION(NVEGSURF), INTENT(IN) ::theta_bioCO2
3026  REAL(KIND(1D0)), DIMENSION(4, NVEGSURF), INTENT(IN) ::LAIPower
3027  REAL(KIND(1D0)), DIMENSION(nsurf + 1, 4, 3), INTENT(IN) ::OHM_coef
3028  REAL(KIND(1D0)), DIMENSION(NSURF + 1, NSURF - 1), INTENT(IN) ::WaterDist
3029  REAL(KIND(1d0)), DIMENSION(:), INTENT(IN) ::Ts5mindata_ir
3030 
3031  ! diurnal profile values for 24hr
3032  REAL(KIND(1D0)), DIMENSION(0:23, 2), INTENT(IN) ::AHProf_24hr
3033  REAL(KIND(1D0)), DIMENSION(0:23, 2), INTENT(IN) ::HumActivity_24hr
3034  REAL(KIND(1D0)), DIMENSION(0:23, 2), INTENT(IN) ::PopProf_24hr
3035  REAL(KIND(1D0)), DIMENSION(0:23, 2), INTENT(IN) ::TraffProf_24hr
3036  REAL(KIND(1D0)), DIMENSION(0:23, 2), INTENT(IN) ::WUProfA_24hr
3037  REAL(KIND(1D0)), DIMENSION(0:23, 2), INTENT(IN) ::WUProfM_24hr
3038  ! ########################################################################################
3039 
3040  ! ########################################################################################
3041  ! inout variables
3042  ! OHM related:
3043  REAL(KIND(1d0)), INTENT(INOUT) ::qn1_av
3044  REAL(KIND(1d0)), INTENT(INOUT) ::dqndt
3045  REAL(KIND(1d0)), INTENT(INOUT) ::qn1_s_av
3046  REAL(KIND(1d0)), INTENT(INOUT) ::dqnsdt
3047 
3048  ! snow related:
3049  REAL(KIND(1D0)), INTENT(INOUT) ::SnowfallCum
3050  REAL(KIND(1D0)), INTENT(INOUT) ::SnowAlb
3051  REAL(KIND(1D0)), DIMENSION(NSURF), INTENT(INOUT) ::IceFrac
3052  REAL(KIND(1D0)), DIMENSION(NSURF), INTENT(INOUT) ::SnowWater
3053  REAL(KIND(1D0)), DIMENSION(NSURF), INTENT(INOUT) ::SnowDens
3054  REAL(KIND(1D0)), DIMENSION(NSURF), INTENT(INOUT) ::SnowFrac
3055  REAL(KIND(1D0)), DIMENSION(NSURF), INTENT(INOUT) ::SnowPack
3056 
3057  ! water balance related:
3058  REAL(KIND(1D0)), DIMENSION(NSURF), INTENT(INOUT) ::soilstore_id
3059  REAL(KIND(1D0)), DIMENSION(NSURF), INTENT(INOUT) ::state_id
3060  REAL(KIND(1D0)), DIMENSION(6, NSURF), INTENT(INOUT) ::StoreDrainPrm
3061 
3062  ! phenology related:
3063  REAL(KIND(1D0)), DIMENSION(NSURF), INTENT(INOUT) ::alb
3064  REAL(KIND(1d0)), DIMENSION(nvegsurf), INTENT(INOUT) ::GDD_id !Growing Degree Days (see SUEWS_DailyState.f95)
3065  REAL(KIND(1d0)), DIMENSION(nvegsurf), INTENT(INOUT) ::SDD_id !Senescence Degree Days (see SUEWS_DailyState.f95)
3066  REAL(KIND(1d0)), DIMENSION(nvegsurf), INTENT(INOUT)::LAI_id !LAI for each veg surface [m2 m-2]
3067  REAL(KIND(1d0)), INTENT(INOUT) :: DecidCap_id
3068  REAL(KIND(1d0)), INTENT(INOUT) :: albDecTr_id
3069  REAL(KIND(1d0)), INTENT(INOUT) :: albEveTr_id
3070  REAL(KIND(1d0)), INTENT(INOUT) :: albGrass_id
3071  REAL(KIND(1d0)), INTENT(INOUT) :: porosity_id
3072  REAL(KIND(1d0)), INTENT(INOUT) :: Tmin_id
3073  REAL(KIND(1d0)), INTENT(INOUT) :: Tmax_id
3074  REAL(KIND(1d0)), INTENT(INOUT) :: lenday_id
3075 
3076  ! anthropogenic heat related:
3077  REAL(KIND(1d0)), DIMENSION(12), INTENT(INOUT) ::HDD_id !Heating Degree Days (see SUEWS_DailyState.f95)
3078 
3079  ! water use related:
3080  REAL(KIND(1d0)), DIMENSION(9), INTENT(INOUT) ::WUDay_id
3081 
3082  ! ESTM related:
3083  REAL(KIND(1d0)), INTENT(INOUT) ::Tair_av
3084  ! ########################################################################################
3085 
3086  ! ########################################################################################
3087  ! output variables
3088  ! REAL(KIND(1D0)),DIMENSION(:,:,:),ALLOCATABLE,INTENT(OUT) ::datetimeBlock
3089  REAL(KIND(1D0)), DIMENSION(len_sim, ncolumnsDataOutSUEWS), INTENT(OUT) ::dataOutBlockSUEWS
3090  REAL(KIND(1D0)), DIMENSION(len_sim, ncolumnsDataOutSnow), INTENT(OUT) ::dataOutBlockSnow
3091  REAL(KIND(1d0)), DIMENSION(len_sim, ncolumnsDataOutESTM), INTENT(OUT) ::dataOutBlockESTM
3092  REAL(KIND(1d0)), DIMENSION(len_sim, ncolumnsDataOutRSL), INTENT(OUT) ::dataOutBlockRSL
3093  REAL(KIND(1d0)), DIMENSION(len_sim, ncolumnsDataOutDailyState), INTENT(OUT) ::DailyStateBlock
3094  ! ########################################################################################
3095 
3096  ! internal temporal iteration related variables
3097  ! INTEGER::dt_since_start ! time since simulation starts [s]
3098 
3099  ! model output blocks of the same size as met forcing block
3100 
3101  ! local variables
3102  ! length of met forcing block
3103  INTEGER :: ir
3104  ! met forcing variables
3105  INTEGER :: iy
3106  INTEGER :: id
3107  INTEGER :: it
3108  INTEGER :: imin
3109  INTEGER :: isec
3110  INTEGER, PARAMETER :: gridiv_x = 1 ! a dummy gridiv as this routine is only one grid
3111  REAL(KIND(1D0))::qn1_obs
3112  REAL(KIND(1D0))::qh_obs
3113  REAL(KIND(1D0))::qe_obs
3114  REAL(KIND(1D0))::qs_obs
3115  REAL(KIND(1D0))::qf_obs
3116  REAL(KIND(1D0))::avu1
3117  REAL(KIND(1D0))::avrh
3118  REAL(KIND(1D0))::Temp_C
3119  REAL(KIND(1D0))::Press_hPa
3120  REAL(KIND(1D0))::Precip
3121  REAL(KIND(1D0))::avkdn
3122  REAL(KIND(1D0))::snowFrac_obs
3123  REAL(KIND(1D0))::ldown_obs
3124  REAL(KIND(1D0))::fcld_obs
3125  REAL(KIND(1D0))::wu_m3
3126  REAL(KIND(1D0))::xsmd
3127  REAL(KIND(1D0))::LAI_obs
3128  REAL(KIND(1D0))::kdiff
3129  REAL(KIND(1D0))::kdir
3130  REAL(KIND(1D0))::wdir
3131 
3132  REAL(KIND(1D0)), DIMENSION(5)::datetimeLine
3133  REAL(KIND(1D0)), DIMENSION(ncolumnsDataOutSUEWS - 5)::dataOutLineSUEWS
3134  REAL(KIND(1D0)), DIMENSION(ncolumnsDataOutSnow - 5)::dataOutLineSnow
3135  REAL(KIND(1d0)), DIMENSION(ncolumnsDataOutESTM - 5)::dataOutLineESTM
3136  REAL(KIND(1d0)), DIMENSION(ncolumnsDataOutRSL - 5)::dataOutLineRSL
3137  REAL(KIND(1d0)), DIMENSION(ncolumnsDataOutDailyState - 5)::DailyStateLine
3138 
3139  REAL(KIND(1D0)), DIMENSION(len_sim, ncolumnsDataOutSUEWS, 1) ::dataOutBlockSUEWS_X
3140  REAL(KIND(1D0)), DIMENSION(len_sim, ncolumnsDataOutSnow, 1) ::dataOutBlockSnow_X
3141  REAL(KIND(1d0)), DIMENSION(len_sim, ncolumnsDataOutESTM, 1) ::dataOutBlockESTM_X
3142  REAL(KIND(1d0)), DIMENSION(len_sim, ncolumnsDataOutRSL, 1) ::dataOutBlockRSL_X
3143  ! REAL(KIND(1d0)),DIMENSION(len_sim,ncolumnsDataOutDailyState,1) ::DailyStateBlock_X
3144 
3145  REAL(KIND(1D0)), DIMENSION(10, 10) ::MetForcingData_grid ! fake array as a placeholder
3146 
3147  CHARACTER(len=150):: FileStateInit
3148  CHARACTER(len=4):: year_txt
3149  CHARACTER(len=3):: id_text
3150  CHARACTER(len=2):: it_text, imin_text
3151 
3152  ! get initial dt_since_start_x from dt_since_start, dt_since_start_x is used for Qn averaging. TS 28 Nov 2018
3153  ! dt_since_start = dt_since_start
3154 
3155  DO ir = 1, len_sim, 1
3156  ! =============================================================================
3157  ! === Translate met data from MetForcingBlock to variable names used in model ==
3158  ! =============================================================================
3159  iy = int(metforcingblock(ir, 1)) !Integer variables
3160  id = int(metforcingblock(ir, 2))
3161  it = int(metforcingblock(ir, 3))
3162  imin = int(metforcingblock(ir, 4))
3163  isec = 0 ! NOT used by SUEWS but by WRF-SUEWS via the cal_main interface
3164  qn1_obs = metforcingblock(ir, 5) !Real values (kind(1d0))
3165  qh_obs = metforcingblock(ir, 6)
3166  qe_obs = metforcingblock(ir, 7)
3167  qs_obs = metforcingblock(ir, 8)
3168  qf_obs = metforcingblock(ir, 9)
3169  avu1 = metforcingblock(ir, 10)
3170  avrh = metforcingblock(ir, 11)
3171  temp_c = metforcingblock(ir, 12)
3172  press_hpa = metforcingblock(ir, 13)
3173  precip = metforcingblock(ir, 14)
3174  avkdn = metforcingblock(ir, 15)
3175  snowfrac_obs = metforcingblock(ir, 16)
3176  ldown_obs = metforcingblock(ir, 17)
3177  fcld_obs = metforcingblock(ir, 18)
3178  wu_m3 = metforcingblock(ir, 19)
3179  xsmd = metforcingblock(ir, 20)
3180  lai_obs = metforcingblock(ir, 21)
3181  kdiff = metforcingblock(ir, 22)
3182  kdir = metforcingblock(ir, 23)
3183  wdir = metforcingblock(ir, 24)
3184 
3185  ! !================================================
3186  ! ! below is for debugging
3187  ! WRITE (year_txt, '(I4)') INT(iy)
3188  ! WRITE (id_text, '(I3)') INT(id)
3189  ! WRITE (it_text, '(I4)') INT(it)
3190  ! WRITE (imin_text, '(I4)') INT(imin)
3191 
3192  ! FileStateInit = './'//TRIM(ADJUSTL(year_txt))//'_'&
3193  ! //TRIM(ADJUSTL(id_text))//'_'&
3194  ! //TRIM(ADJUSTL(it_text))//'_'&
3195  ! //TRIM(ADJUSTL(imin_text))//'_'&
3196  ! //'state_init.nml'
3197 
3198  ! OPEN (12, file=FileStateInit, position='rewind')
3199 
3200  ! write (12, *) '&state_init'
3201  ! write (12, *) 'aerodynamicresistancemethod=', aerodynamicresistancemethod
3202  ! write (12, *) 'ah_min=', ah_min
3203  ! write (12, *) 'ahprof_24hr=', ahprof_24hr
3204  ! write (12, *) 'ah_slope_cooling=', ah_slope_cooling
3205  ! write (12, *) 'ah_slope_heating=', ah_slope_heating
3206  ! write (12, *) 'alb=', alb
3207  ! write (12, *) 'albmax_dectr=', albmax_dectr
3208  ! write (12, *) 'albmax_evetr=', albmax_evetr
3209  ! write (12, *) 'albmax_grass=', albmax_grass
3210  ! write (12, *) 'albmin_dectr=', albmin_dectr
3211  ! write (12, *) 'albmin_evetr=', albmin_evetr
3212  ! write (12, *) 'albmin_grass=', albmin_grass
3213  ! write (12, *) 'alpha_bioco2=', alpha_bioco2
3214  ! write (12, *) 'alpha_enh_bioco2=', alpha_enh_bioco2
3215  ! write (12, *) 'alt=', alt
3216  ! write (12, *) 'avkdn=', avkdn
3217  ! write (12, *) 'avrh=', avrh
3218  ! write (12, *) 'avu1=', avu1
3219  ! write (12, *) 'baset=', baset
3220  ! write (12, *) 'basete=', basete
3221  ! write (12, *) 'basethdd=', basethdd
3222  ! write (12, *) 'beta_bioco2=', beta_bioco2
3223  ! write (12, *) 'beta_enh_bioco2=', beta_enh_bioco2
3224  ! write (12, *) 'bldgh=', bldgh
3225  ! write (12, *) 'capmax_dec=', capmax_dec
3226  ! write (12, *) 'capmin_dec=', capmin_dec
3227  ! write (12, *) 'chanohm=', chanohm
3228  ! write (12, *) 'co2pointsource=', co2pointsource
3229  ! write (12, *) 'cpanohm=', cpanohm
3230  ! write (12, *) 'crwmax=', crwmax
3231  ! write (12, *) 'crwmin=', crwmin
3232  ! write (12, *) 'daywat=', daywat
3233  ! write (12, *) 'daywatper=', daywatper
3234  ! write (12, *) 'dectreeh=', dectreeh
3235  ! write (12, *) 'diagnose=', diagnose
3236  ! write (12, *) 'diagqn=', diagqn
3237  ! write (12, *) 'diagqs=', diagqs
3238  ! write (12, *) 'drainrt=', drainrt
3239  ! write (12, *) 'dt_since_start=', dt_since_start
3240  ! write (12, *) 'dqndt=', dqndt
3241  ! write (12, *) 'qn1_av=', qn1_av
3242  ! write (12, *) 'dqnsdt=', dqnsdt
3243  ! write (12, *) 'qn1_s_av=', qn1_s_av
3244  ! write (12, *) 'ef_umolco2perj=', ef_umolco2perj
3245  ! write (12, *) 'emis=', emis
3246  ! write (12, *) 'emissionsmethod=', emissionsmethod
3247  ! write (12, *) 'enef_v_jkm=', enef_v_jkm
3248  ! write (12, *) 'enddls=', enddls
3249  ! write (12, *) 'evetreeh=', evetreeh
3250  ! write (12, *) 'faibldg=', faibldg
3251  ! write (12, *) 'faidectree=', faidectree
3252  ! write (12, *) 'faievetree=', faievetree
3253  ! write (12, *) 'faut=', faut
3254  ! write (12, *) 'fcef_v_kgkm=', fcef_v_kgkm
3255  ! write (12, *) 'fcld_obs=', fcld_obs
3256  ! write (12, *) 'flowchange=', flowchange
3257  ! write (12, *) 'frfossilfuel_heat=', frfossilfuel_heat
3258  ! write (12, *) 'frfossilfuel_nonheat=', frfossilfuel_nonheat
3259  ! write (12, *) 'g1=', g1
3260  ! write (12, *) 'g2=', g2
3261  ! write (12, *) 'g3=', g3
3262  ! write (12, *) 'g4=', g4
3263  ! write (12, *) 'g5=', g5
3264  ! write (12, *) 'g6=', g6
3265  ! write (12, *) 'gdd_id=', gdd_id
3266  ! write (12, *) 'gddfull=', gddfull
3267  ! write (12, *) 'gridiv=', gridiv
3268  ! write (12, *) 'gsmodel=', gsmodel
3269  ! write (12, *) 'hdd_id=', hdd_id
3270  ! write (12, *) 'humactivity_24hr=', humactivity_24hr
3271  ! write (12, *) 'icefrac=', icefrac
3272  ! write (12, *) 'id=', id
3273  ! write (12, *) 'ie_a=', ie_a
3274  ! write (12, *) 'ie_end=', ie_end
3275  ! write (12, *) 'ie_m=', ie_m
3276  ! write (12, *) 'ie_start=', ie_start
3277  ! write (12, *) 'imin=', imin
3278  ! write (12, *) 'internalwateruse_h=', internalwateruse_h
3279  ! write (12, *) 'irrfracconif=', irrfracconif
3280  ! write (12, *) 'irrfracdecid=', irrfracdecid
3281  ! write (12, *) 'irrfracgrass=', irrfracgrass
3282  ! write (12, *) 'isec=', isec
3283  ! write (12, *) 'it=', it
3284  ! write (12, *) 'evapmethod=', evapmethod
3285  ! write (12, *) 'iy=', iy
3286  ! write (12, *) 'kkanohm=', kkanohm
3287  ! write (12, *) 'kmax=', kmax
3288  ! write (12, *) 'lai_id=', lai_id
3289  ! write (12, *) 'laicalcyes=', laicalcyes
3290  ! write (12, *) 'laimax=', laimax
3291  ! write (12, *) 'laimin=', laimin
3292  ! write (12, *) 'lai_obs=', lai_obs
3293  ! write (12, *) 'laipower=', laipower
3294  ! write (12, *) 'laitype=', laitype
3295  ! write (12, *) 'lat=', lat
3296  ! write (12, *) 'lenday_id=', lenday_id
3297  ! write (12, *) 'ldown_obs=', ldown_obs
3298  ! write (12, *) 'lng=', lng
3299  ! write (12, *) 'maxconductance=', maxconductance
3300  ! write (12, *) 'maxfcmetab=', maxfcmetab
3301  ! write (12, *) 'maxqfmetab=', maxqfmetab
3302  ! write (12, *) 'snowwater=', snowwater
3303  ! ! write (12, *) 'metforcingdata_grid=', metforcingdata_grid
3304  ! write (12, *) 'minfcmetab=', minfcmetab
3305  ! write (12, *) 'minqfmetab=', minqfmetab
3306  ! write (12, *) 'min_res_bioco2=', min_res_bioco2
3307  ! write (12, *) 'narp_emis_snow=', narp_emis_snow
3308  ! write (12, *) 'narp_trans_site=', narp_trans_site
3309  ! write (12, *) 'netradiationmethod=', netradiationmethod
3310  ! write (12, *) 'ohm_coef=', ohm_coef
3311  ! write (12, *) 'ohmincqf=', ohmincqf
3312  ! write (12, *) 'ohm_threshsw=', ohm_threshsw
3313  ! write (12, *) 'ohm_threshwd=', ohm_threshwd
3314  ! write (12, *) 'pipecapacity=', pipecapacity
3315  ! write (12, *) 'popdensdaytime=', popdensdaytime
3316  ! write (12, *) 'popdensnighttime=', popdensnighttime
3317  ! write (12, *) 'popprof_24hr=', popprof_24hr
3318  ! write (12, *) 'pormax_dec=', pormax_dec
3319  ! write (12, *) 'pormin_dec=', pormin_dec
3320  ! write (12, *) 'precip=', precip
3321  ! write (12, *) 'preciplimit=', preciplimit
3322  ! write (12, *) 'preciplimitalb=', preciplimitalb
3323  ! write (12, *) 'press_hpa=', press_hpa
3324  ! write (12, *) 'qf0_beu=', qf0_beu
3325  ! write (12, *) 'qf_a=', qf_a
3326  ! write (12, *) 'qf_b=', qf_b
3327  ! write (12, *) 'qf_c=', qf_c
3328  ! write (12, *) 'qn1_obs=', qn1_obs
3329  ! write (12, *) 'qh_obs=', qh_obs
3330  ! write (12, *) 'qs_obs=', qs_obs
3331  ! write (12, *) 'qf_obs=', qf_obs
3332  ! write (12, *) 'radmeltfact=', radmeltfact
3333  ! write (12, *) 'raincover=', raincover
3334  ! write (12, *) 'rainmaxres=', rainmaxres
3335  ! write (12, *) 'resp_a=', resp_a
3336  ! write (12, *) 'resp_b=', resp_b
3337  ! write (12, *) 'roughlenheatmethod=', roughlenheatmethod
3338  ! write (12, *) 'roughlenmommethod=', roughlenmommethod
3339  ! write (12, *) 'runofftowater=', runofftowater
3340  ! write (12, *) 's1=', s1
3341  ! write (12, *) 's2=', s2
3342  ! write (12, *) 'sathydraulicconduct=', sathydraulicconduct
3343  ! write (12, *) 'sddfull=', sddfull
3344  ! write (12, *) 'sdd_id=', sdd_id
3345  ! write (12, *) 'sfr=', sfr
3346  ! write (12, *) 'smdmethod=', smdmethod
3347  ! write (12, *) 'snowalb=', snowalb
3348  ! write (12, *) 'snowalbmax=', snowalbmax
3349  ! write (12, *) 'snowalbmin=', snowalbmin
3350  ! write (12, *) 'snowpacklimit=', snowpacklimit
3351  ! write (12, *) 'snowdens=', snowdens
3352  ! write (12, *) 'snowdensmax=', snowdensmax
3353  ! write (12, *) 'snowdensmin=', snowdensmin
3354  ! write (12, *) 'snowfallcum=', snowfallcum
3355  ! write (12, *) 'snowfrac=', snowfrac
3356  ! write (12, *) 'snowlimbldg=', snowlimbldg
3357  ! write (12, *) 'snowlimpaved=', snowlimpaved
3358  ! write (12, *) 'snowfrac_obs=', snowfrac_obs
3359  ! write (12, *) 'snowpack=', snowpack
3360  ! write (12, *) 'snowprof_24hr=', snowprof_24hr
3361  ! write (12, *) 'snowuse=', snowuse
3362  ! write (12, *) 'soildepth=', soildepth
3363  ! write (12, *) 'soilstore_id=', soilstore_id
3364  ! write (12, *) 'soilstorecap=', soilstorecap
3365  ! write (12, *) 'stabilitymethod=', stabilitymethod
3366  ! write (12, *) 'startdls=', startdls
3367  ! write (12, *) 'state_id=', state_id
3368  ! write (12, *) 'statelimit=', statelimit
3369  ! write (12, *) 'storageheatmethod=', storageheatmethod
3370  ! write (12, *) 'storedrainprm=', storedrainprm
3371  ! write (12, *) 'surfacearea=', surfacearea
3372  ! write (12, *) 'tair_av=', tair_av
3373  ! write (12, *) 'tau_a=', tau_a
3374  ! write (12, *) 'tau_f=', tau_f
3375  ! write (12, *) 'tau_r=', tau_r
3376  ! write (12, *) 'tmax_id=', tmax_id
3377  ! write (12, *) 'tmin_id=', tmin_id
3378  ! write (12, *) 't_critic_cooling=', t_critic_cooling
3379  ! write (12, *) 't_critic_heating=', t_critic_heating
3380  ! write (12, *) 'temp_c=', temp_c
3381  ! write (12, *) 'tempmeltfact=', tempmeltfact
3382  ! write (12, *) 'th=', th
3383  ! write (12, *) 'theta_bioco2=', theta_bioco2
3384  ! write (12, *) 'timezone=', timezone
3385  ! write (12, *) 'tl=', tl
3386  ! write (12, *) 'trafficrate=', trafficrate
3387  ! write (12, *) 'trafficunits=', trafficunits
3388  ! write (12, *) 'traffprof_24hr=', traffprof_24hr
3389  ! ! write (12, *) 'ts5mindata_ir=', ts5mindata_ir
3390  ! write (12, *) 'tstep=', tstep
3391  ! write (12, *) 'tstep_prev=', tstep_prev
3392  ! write (12, *) 'veg_type=', veg_type
3393  ! write (12, *) 'waterdist=', waterdist
3394  ! write (12, *) 'waterusemethod=', waterusemethod
3395  ! write (12, *) 'wetthresh=', wetthresh
3396  ! write (12, *) 'wu_m3=', wu_m3
3397  ! write (12, *) 'wuday_id=', wuday_id
3398  ! write (12, *) 'decidcap_id=', decidcap_id
3399  ! write (12, *) 'albdectr_id=', albdectr_id
3400  ! write (12, *) 'albevetr_id=', albevetr_id
3401  ! write (12, *) 'albgrass_id=', albgrass_id
3402  ! write (12, *) 'porosity_id=', porosity_id
3403  ! write (12, *) 'wuprofa_24hr=', wuprofa_24hr
3404  ! write (12, *) 'wuprofm_24hr=', wuprofm_24hr
3405  ! write (12, *) 'xsmd=', xsmd
3406  ! write (12, *) 'z=', z
3407  ! write (12, *) 'z0m_in=', z0m_in
3408  ! write (12, *) 'zdm_in=', zdm_in
3409  ! write (12, *) '/'
3410 
3411  ! WRITE (12, *) ''
3412 
3413  ! CLOSE (12)
3414  ! !================================================
3415 
3416  CALL suews_cal_main( &
3417  aerodynamicresistancemethod, ah_min, ahprof_24hr, ah_slope_cooling, & ! input&inout in alphabetical order
3418  ah_slope_heating, &
3419  alb, albmax_dectr, albmax_evetr, albmax_grass, &
3420  albmin_dectr, albmin_evetr, albmin_grass, &
3421  alpha_bioco2, alpha_enh_bioco2, alt, avkdn, avrh, avu1, baset, basete, &
3422  basethdd, beta_bioco2, beta_enh_bioco2, bldgh, capmax_dec, capmin_dec, &
3423  chanohm, co2pointsource, cpanohm, crwmax, crwmin, daywat, daywatper, &
3424  dectreeh, diagnose, diagqn, diagqs, drainrt, &
3425  dt_since_start, dqndt, qn1_av, dqnsdt, qn1_s_av, &
3426  ef_umolco2perj, emis, emissionsmethod, enef_v_jkm, enddls, evetreeh, faibldg, &
3427  faidectree, faievetree, faut, fcef_v_kgkm, fcld_obs, flowchange, &
3428  frfossilfuel_heat, frfossilfuel_nonheat, g1, g2, g3, g4, g5, g6, gdd_id, &
3429  gddfull, gridiv, gsmodel, hdd_id, humactivity_24hr, &
3430  icefrac, id, ie_a, ie_end, ie_m, ie_start, imin, &
3431  internalwateruse_h, irrfracconif, irrfracdecid, irrfracgrass, isec, it, evapmethod, &
3432  iy, kkanohm, kmax, lai_id, laicalcyes, laimax, laimin, lai_obs, &
3433  laipower, laitype, lat, lenday_id, ldown_obs, lng, maxconductance, maxfcmetab, maxqfmetab, &
3434  snowwater, metforcingdata_grid, minfcmetab, minqfmetab, min_res_bioco2, &
3435  narp_emis_snow, narp_trans_site, netradiationmethod, &
3436  ohm_coef, ohmincqf, ohm_threshsw, &
3437  ohm_threshwd, pipecapacity, popdensdaytime, &
3438  popdensnighttime, popprof_24hr, pormax_dec, pormin_dec, &
3439  precip, preciplimit, preciplimitalb, press_hpa, &
3440  qf0_beu, qf_a, qf_b, qf_c, &
3441  qn1_obs, qs_obs, qf_obs, &
3442  radmeltfact, raincover, rainmaxres, resp_a, resp_b, &
3443  roughlenheatmethod, roughlenmommethod, runofftowater, s1, s2, &
3444  sathydraulicconduct, sddfull, sdd_id, sfr, smdmethod, snowalb, snowalbmax, &
3445  snowalbmin, snowpacklimit, snowdens, snowdensmax, snowdensmin, snowfallcum, snowfrac, &
3446  snowlimbldg, snowlimpaved, snowfrac_obs, snowpack, snowprof_24hr, snowuse, soildepth, &
3447  soilstore_id, soilstorecap, stabilitymethod, startdls, state_id, statelimit, &
3448  storageheatmethod, storedrainprm, surfacearea, tair_av, tau_a, tau_f, tau_r, &
3449  tmax_id, tmin_id, &
3450  t_critic_cooling, t_critic_heating, temp_c, tempmeltfact, th, &
3451  theta_bioco2, timezone, tl, trafficrate, trafficunits, &
3452  traffprof_24hr, ts5mindata_ir, tstep, tstep_prev, veg_type, &
3453  waterdist, waterusemethod, wetthresh, wu_m3, &
3454  wuday_id, decidcap_id, albdectr_id, albevetr_id, albgrass_id, porosity_id, &
3455  wuprofa_24hr, wuprofm_24hr, xsmd, z, z0m_in, zdm_in, &
3456  datetimeline, dataoutlinesuews, dataoutlinesnow, dataoutlineestm, dataoutlinersl, &!output
3457  dailystateline)!output
3458 
3459  ! update dt_since_start_x for next iteration, dt_since_start_x is used for Qn averaging. TS 28 Nov 2018
3460  dt_since_start = dt_since_start + tstep
3461 
3462  !============ update DailyStateBlock ===============
3463  dailystateblock(ir, :) = [datetimeline, dailystateline]
3464 
3465  !============ write out results ===============
3466  ! works at each timestep
3467  CALL suews_update_output( &
3468  snowuse, storageheatmethod, &!input
3469  len_sim, 1, &
3470  ir, gridiv_x, datetimeline, dataoutlinesuews, dataoutlinesnow, dataoutlineestm, dataoutlinersl, &!input
3471  dataoutblocksuews_x, dataoutblocksnow_x, dataoutblockestm_x, dataoutblockrsl_x)!inout
3472 
3473  END DO
3474 
3475  dataoutblocksuews = dataoutblocksuews_x(:, :, 1)
3476  dataoutblocksnow = dataoutblocksnow_x(:, :, 1)
3477  dataoutblockestm = dataoutblockestm_x(:, :, 1)
3478  dataoutblockrsl = dataoutblockrsl_x(:, :, 1)
3479  ! DailyStateBlock=DailyStateBlock_X(:,:,1)
3480 
3481  END SUBROUTINE suews_cal_multitsteps
3482 
3483  ! a wrapper of NARP_cal_SunPosition used by supy
3484  SUBROUTINE suews_cal_sunposition( &
3485  year, idectime, UTC, locationlatitude, locationlongitude, locationaltitude, & !input
3486  sunazimuth, sunzenith) !output
3487  IMPLICIT NONE
3488 
3489  REAL(KIND(1D0)), INTENT(in) :: year, idectime, UTC, &
3490  locationlatitude, locationlongitude, locationaltitude
3491  REAL(KIND(1D0)), INTENT(out) ::sunazimuth, sunzenith
3492 
3493  CALL narp_cal_sunposition( &
3494  year, idectime, utc, locationlatitude, locationlongitude, locationaltitude, &
3495  sunazimuth, sunzenith)
3496 
3497  END SUBROUTINE suews_cal_sunposition
3498 
3499  ! function func(arg) result(retval)
3500  ! implicit none
3501  ! type :: arg
3502  ! type :: retval
3503 
3504  ! end function func
3505 
3506  function cal_tair_av(tair_av_prev, dt_since_start, tstep, temp_c) result(tair_av_next)
3507  ! calculate mean air temperature of past 24 hours
3508  ! TS, 17 Sep 2019
3509  implicit none
3510  real(KIND(1D0)), intent(in) :: tair_av_prev
3511  real(KIND(1D0)), intent(in) :: temp_c
3512  integer, intent(in) :: dt_since_start
3513  integer, intent(in) :: tstep
3514 
3515  real(KIND(1D0)) ::tair_av_next
3516 
3517  real(KIND(1D0)), parameter:: len_day_s = 24*3600 ! day length in seconds
3518  real(KIND(1D0)):: len_cal_s ! length of average period in seconds
3519  real(KIND(1D0)):: temp_k ! temp in K
3520 
3521  ! determine the average period
3522  if (dt_since_start > len_day_s) then
3523  ! if simulation has been running over one day
3524  len_cal_s = len_day_s
3525  else
3526  ! if simulation has been running less than one day
3527  len_cal_s = dt_since_start + tstep
3528  end if
3529  temp_k = temp_c + 273.15
3530  tair_av_next = tair_av_prev*(len_cal_s - tstep*1.)/len_cal_s + temp_k*tstep/len_cal_s
3531 
3532  end function cal_tair_av
3533 
3534  function cal_tsfc(qh, avdens, avcp, RA, temp_c) result(tsfc_C)
3535  ! calculate surface/skin temperature
3536  ! TS, 23 Oct 2019
3537  implicit none
3538  real(KIND(1D0)), intent(in) :: qh ! sensible heat flux [W m-2]
3539  real(KIND(1D0)), intent(in) :: avdens ! air density [kg m-3]
3540  real(KIND(1D0)), intent(in) :: avcp !air heat capacity [J m-3 K-1]
3541  real(KIND(1D0)), intent(in) :: RA !Aerodynamic resistance [s m^-1]
3542  real(KIND(1D0)), intent(in) :: temp_C ! air temperature [C]
3543 
3544  real(KIND(1D0)) ::tsfc_C ! surface temperature [C]
3545 
3546  tsfc_c = qh/(avdens*avcp)*ra + temp_c
3547  end function cal_tsfc
3548 
3549 END MODULE suews_driver
subroutine suews_update_soilmoist(NonWaterFraction, SoilStoreCap, sfr, soilstore_id, SoilMoistCap, SoilState, vsmd, smd)
subroutine anohm(tstep, dt_since_start, qn1, qn1_av_prev, dqndt_prev, qf, MetForcingData_grid, moist_surf, alb, emis, cpAnOHM, kkAnOHM, chAnOHM, sfr, nsurf, EmissionsMethod, id, Gridiv, qn1_av_next, dqndt_next, a1, a2, a3, qs, deltaQi)
High level wrapper for AnOHM calculation.
subroutine suews_cal_soilstate(SMDMethod, xsmd, NonWaterFraction, SoilMoistCap, SoilStoreCap, surf_chang_per_tstep, soilstore_id, soilstoreOld, sfr, smd, smd_nsurf, tot_chang_per_tstep, SoilState)
subroutine co2_biogen(alpha_bioCO2, alpha_enh_bioCO2, avkdn, beta_bioCO2, beta_enh_bioCO2, BSoilSurf, ConifSurf, DecidSurf, dectime, EmissionsMethod, gfunc, gfunc2, GrassSurf, gsmodel, id, it, ivConif, ivDecid, ivGrass, LAI_id, LAIMin, LAIMax, min_res_bioCO2, nsurf, NVegSurf, resp_a, resp_b, sfr, SnowFrac, t2, Temp_C, theta_bioCO2, Fc_biogen, Fc_photo, Fc_respi)
subroutine snowcalc(tstep, imin, it, dectime, is, EvapMethod, CRWmin, CRWmax, nsh_real, lvS_J_kg, lv_j_kg, avdens, avRh, Press_hPa, Temp_C, RAsnow, psyc_hPa, avcp, sIce_hPa, PervFraction, vegfraction, addimpervious, vpd_hPa, qn_e, s_hPa, ResistSurf, RA, rb, tlv, snowdensmin, SnowProf_24hr, precip, PipeCapacity, RunoffToWater, addVeg, SnowLimPaved, SnowLimBldg, FlowChange, drain, WetThresh, stateOld, mw_ind, SoilStoreCap, rainonsnow, freezmelt, freezstate, freezstatevol, Qm_Melt, Qm_rain, Tsurf_ind, sfr, dayofWeek_id, StoreDrainPrm, SnowPackLimit, AddWater, addwaterrunoff, soilstore_id, SnowPack, SurplusEvap, SnowFrac, SnowWater, iceFrac, SnowDens, runoffAGimpervious, runoffAGveg, surplusWaterBody, rss_nsurf, runoffSnow, runoff, runoffSoil, chang, changSnow, SnowToSurf, state_id, ev_snow, SnowDepth, SnowRemoval, swe, ev, chSnow_tot, ev_tot, qe_tot, runoff_tot, surf_chang_tot, runoffPipes, mwstore, runoffwaterbody)
subroutine redistributewater(snowUse, WaterDist, sfr, Drain, AddWaterRunoff, AddWater)
subroutine suews_cal_biogenco2(alpha_bioCO2, alpha_enh_bioCO2, avkdn, avRh, beta_bioCO2, beta_enh_bioCO2, BSoilSurf, ConifSurf, DecidSurf, dectime, Diagnose, EmissionsMethod, Fc_anthro, G1, G2, G3, G4, G5, G6, gfunc, GrassSurf, gsmodel, id, it, ivConif, ivDecid, ivGrass, Kmax, LAI_id, LAIMin, LAIMax, MaxConductance, min_res_bioCO2, nsurf, NVegSurf, Press_hPa, resp_a, resp_b, S1, S2, sfr, SMDMethod, SnowFrac, t2_C, Temp_C, theta_bioCO2, TH, TL, vsmd, xsmd, Fc, Fc_biogen, Fc_photo, Fc_respi)
integer, parameter ncolumnsdataoutestm
subroutine suews_init_qh(avdens, avcp, h_mod, qn1, dectime, H_init)
integer, parameter ivgrass
subroutine suews_cal_dailystate(iy, id, it, imin, isec, tstep, tstep_prev, dt_since_start, DayofWeek_id, Tmin_id_prev, Tmax_id_prev, lenDay_id_prev, WaterUseMethod, Ie_start, Ie_end, LAICalcYes, LAIType, nsh_real, avkdn, Temp_C, Precip, BaseTHDD, lat, Faut, LAI_obs, AlbMax_DecTr, AlbMax_EveTr, AlbMax_Grass, AlbMin_DecTr, AlbMin_EveTr, AlbMin_Grass, CapMax_dec, CapMin_dec, PorMax_dec, PorMin_dec, Ie_a, Ie_m, DayWatPer, DayWat, BaseT, BaseTe, GDDFull, SDDFull, LAIMin, LAIMax, LAIPower, DecidCap_id_prev, StoreDrainPrm_prev, LAI_id_prev, GDD_id_prev, SDD_id_prev, albDecTr_id_prev, albEveTr_id_prev, albGrass_id_prev, porosity_id_prev, HDD_id_prev, HDD_id_next, Tmin_id_next, Tmax_id_next, lenDay_id_next, albDecTr_id_next, albEveTr_id_next, albGrass_id_next, porosity_id_next, DecidCap_id_next, StoreDrainPrm_next, LAI_id_next, GDD_id_next, SDD_id_next, deltaLAI, WUDay_id)
subroutine narp_cal_sunposition(year, idectime, UTC, locationlatitude, locationlongitude, locationaltitude, sunazimuth, sunzenith)
integer, parameter ncolumnsdataoutdailystate
integer, parameter bsoilsurf
subroutine suews_cal_water(Diagnose, snowUse, NonWaterFraction, addPipes, addImpervious, addVeg, addWaterBody, state_id, soilstore_id, sfr, StoreDrainPrm, WaterDist, nsh_real, drain_per_tstep, drain, AddWaterRunoff, AdditionalWater, runoffPipes, runoff_per_interval, AddWater, stateOld, soilstoreOld)
subroutine estm(Gridiv, tstep, avkdn, avu1, temp_c, zenith_deg, avrh, press_hpa, ldown, bldgh, Ts5mindata_ir, Tair_av, dataOutLineESTM, QS)
real(kind(1d0)) function qa2rh(qa_gkg, pres_hPa, Ta_degC)
subroutine snow_cal_meltheat(snowUse, tstep, tau_r, SnowDensMax, lvS_J_kg, lv_J_kg, tstep_real, RadMeltFact, TempMeltFact, SnowAlbMax, SnowDensMin, Temp_C, Precip, PrecipLimit, PrecipLimitAlb, nsh_real, sfr, Tsurf_ind, Tsurf_ind_snow, state_id, qn1_ind_snow, kup_ind_snow, SnowWater, deltaQi, alb1, SnowPack_in, SnowFrac_in, SnowAlb_in, SnowDens_in, SnowfallCum_in, SnowPack_out, SnowFrac_out, SnowAlb_out, SnowDens_out, SnowfallCum_out, mwh, Qm, QmFreez, QmRain, veg_fr, snowCalcSwitch, Qm_melt, Qm_freezState, Qm_rain, FreezMelt, FreezState, FreezStateVol, rainOnSnow, SnowDepth, mw_ind, dataOutLineSnow)
subroutine radmethod(NetRadiationMethod, snowUse, NetRadiationMethod_use, AlbedoChoice, ldown_option)
subroutine suews_cal_sunposition(year, idectime, UTC, locationlatitude, locationlongitude, locationaltitude, sunazimuth, sunzenith)
elemental real(kind(1d0)) function set_nan(x)
subroutine boundarylayerresistance(zzd, z0m, avU1, UStar, rb)
subroutine suews_update_dailystate(id, datetimeline, Gridiv, NumberOfGrids, DailyStateLine, dataOutDailyState)
subroutine cal_atmmoist(Temp_C, Press_hPa, avRh, dectime, lv_J_kg, lvS_J_kg, es_hPa, Ea_hPa, VPd_hpa, VPD_Pa, dq, dens_dry, avcp, air_dens)
integer, parameter nsurf
subroutine lumps_cal_qhqe(veg_type, snowUse, qn1, qf, qs, Qm, Temp_C, Veg_Fr, avcp, Press_hPa, lv_J_kg, tstep_real, DRAINRT, nsh_real, Precip, RainMaxRes, RAINCOVER, sfr, LAI_id_prev, LAImax, LAImin, QH_LUMPS, QE_LUMPS, psyc_hPa, s_hPa, sIce_hpa, Veg_Fr_temp, VegPhenLumps)
subroutine suews_update_output(SnowUse, storageheatmethod, ReadLinesMetdata, NumberOfGrids, ir, gridiv, datetimeLine, dataOutLineSUEWS, dataOutLineSnow, dataOutLineESTM, dataoutLineURSL, dataOutSUEWS, dataOutSnow, dataOutESTM, dataOutRSL)
subroutine suews_cal_anthropogenicemission(AH_MIN, AHProf_24hr, AH_SLOPE_Cooling, AH_SLOPE_Heating, CO2PointSource, dayofWeek_id, Diagnose, DLS, EF_umolCO2perJ, EmissionsMethod, EnEF_v_Jkm, FcEF_v_kgkm, FrFossilFuel_Heat, FrFossilFuel_NonHeat, HDD_id, HumActivity_24hr, id, imin, it, MaxFCMetab, MaxQFMetab, MinFCMetab, MinQFMetab, nsh, PopDensDaytime, PopDensNighttime, PopProf_24hr, QF, QF0_BEU, Qf_A, Qf_B, Qf_C, QF_obs, QF_SAHP, sfr, SnowFrac, SurfaceArea, T_CRITIC_Cooling, T_CRITIC_Heating, Temp_C, TrafficRate, TrafficUnits, TraffProf_24hr, Fc_anthro, Fc_build, Fc_metab, Fc_point, Fc_traff)
real(kind(1d0)) lv_j_kg
subroutine snowupdate(tstep, Temp_C, tau_a, tau_f, tau_r, SnowDensMax, SnowDensMin, SnowAlbMax, SnowAlbMin, SnowPack_prev, SnowAlb_prev, SnowDens_prev, SnowAlb_next, SnowDens_next)
subroutine suews_cal_qn(NetRadiationMethod, snowUse, tstep, SnowPack_prev, tau_a, tau_f, SnowAlbMax, SnowAlbMin, Diagnose, snowFrac_obs, ldown_obs, fcld_obs, dectime, ZENITH_deg, Tsurf_0, avKdn, Temp_C, avRH, ea_hPa, qn1_obs, SnowAlb_prev, snowFrac_prev, DiagQN, NARP_TRANS_SITE, NARP_EMIS_SNOW, IceFrac, sfr, emis, alb_prev, albDecTr_id, albEveTr_id, albGrass_id, alb_next, ldown, fcld, qn1, qn1_snowfree, qn1_S, kclear, kup, lup, tsurf, qn1_ind_snow, kup_ind_snow, Tsurf_ind_snow, Tsurf_ind, alb1, snowFrac_next, SnowAlb_next)
subroutine ohm(qn1, qn1_av_prev, dqndt_prev, qn1_av_next, dqndt_next, qn1_S, qn1_s_av_prev, dqnsdt_prev, qn1_s_av_next, dqnsdt_next, tstep, dt_since_start, sfr, nsurf, Tair_mav_5d, OHM_coef, OHM_threshSW, OHM_threshWD, soilstore_id, SoilStoreCap, state_id, BldgSurf, WaterSurf, SnowUse, SnowFrac, DiagQS, a1, a2, a3, qs, deltaQi)
subroutine cal_stab(StabilityMethod, zzd, z0m, zdm, avU1, Temp_C, QH_init, avdens, avcp, L_MOD, TStar, UStar, zL)
integer, parameter conifsurf
subroutine suews_cal_tstep(tstep, nsh, nsh_real, tstep_real)
subroutine suews_cal_qh(QHMethod, qn1, qf, QmRain, qeOut, qs, QmFreez, qm, avdens, avcp, tsurf, Temp_C, RA, qh, qh_residual, qh_resist)
subroutine suews_cal_resistance(StabilityMethod, Diagnose, AerodynamicResistanceMethod, RoughLenHeatMethod, snowUse, id, it, gsModel, SMDMethod, avdens, avcp, QH_init, zzd, z0m, zdm, avU1, Temp_C, VegFraction, avkdn, Kmax, G1, G2, G3, G4, G5, G6, S1, S2, TH, TL, dq, xsmd, vsmd, MaxConductance, LAIMax, LAI_id, SnowFrac, sfr, UStar, TStar, L_mod, zL, gsc, ResistSurf, RA, RAsnow, rb)
subroutine surfaceresistance(id, it, SMDMethod, SnowFrac, sfr, avkdn, Temp_C, dq, xsmd, vsmd, MaxConductance, LAIMax, LAI_id, gsModel, Kmax, G1, G2, G3, G4, G5, G6, TH, TL, S1, S2, gfunc, gsc, ResistSurf)
real(kind(1d0)) avcp
real(kind(1d0)) function update_snow_albedo(tstep, SnowPack_prev, SnowAlb_prev, Temp_C, tau_a, tau_f, SnowAlbMax, SnowAlbMin)
subroutine suews_update_outputline(AdditionalWater, alb, avkdn, avU10_ms, azimuth, chSnow_per_interval, dectime, drain_per_tstep, E_mod, ev_per_tstep, ext_wu, Fc, Fc_build, fcld, Fc_metab, Fc_photo, Fc_respi, Fc_point, Fc_traff, FlowChange, h_mod, id, imin, int_wu, it, iy, kup, LAI_id, ldown, l_mod, lup, mwh, MwStore, nsh_real, NWstate_per_tstep, Precip, q2_gkg, qeOut, qf, qh, qh_resist, Qm, QmFreez, QmRain, qn1, qn1_S, qn1_snowfree, qs, RA, resistsurf, RH2, runoffAGimpervious, runoffAGveg, runoff_per_tstep, runoffPipes, runoffSoil_per_tstep, runoffWaterBody, sfr, smd, smd_nsurf, SnowAlb, SnowRemoval, state_id, state_per_tstep, surf_chang_per_tstep, swe, t2_C, tskin_C, tot_chang_per_tstep, tsurf, UStar, wu_DecTr, wu_EveTr, wu_Grass, z0m, zdm, zenith_deg, datetimeLine, dataOutLineSUEWS)
AnOHM: Analytical Objective Hysteresis Model.
real(kind(1d0)) function square(x)
subroutine update_dailystateline(it, imin, nsh_real, GDD_id, HDD_id, LAI_id, SDD_id, Tmin_id, Tmax_id, lenday_id, DecidCap_id, albDecTr_id, albEveTr_id, albGrass_id, porosity_id, WUDay_id, deltaLAI, VegPhenLumps, SnowAlb, SnowDens, a1, a2, a3, DailyStateLine)
subroutine soilstore(is, sfr, PipeCapacity, RunoffToWater, pin, wu_EveTr, wu_DecTr, wu_Grass, drain, AddWater, addImpervious, nsh_real, stateOld, AddWaterRunoff, PervFraction, addVeg, SoilStoreCap, addWaterBody, FlowChange, StateLimit, runoffAGimpervious, surplusWaterBody, runoffAGveg, runoffPipes, ev, soilstore_id, SurplusEvap, runoffWaterBody, p_mm, chang, runoff, state_id)
real(kind(1d0)) function cal_z0v(RoughLenHeatMethod, z0m, VegFraction, UStar)
real(kind(1d0)) function stab_psi_heat(StabilityMethod, ZL)
integer, parameter ncolumnsdataoutrsl
integer, parameter ncolumnsdataoutsnow
subroutine rslprofile(Zh, z0m, zdm, L_MOD, sfr, planF, StabilityMethod, avcp, lv_J_kg, avdens, avU1, Temp_C, avRH, Press_hPa, zMeas, qh, qe, T2_C, q2_gkg, U10_ms, RH2, dataoutLineRSL)
integer, parameter grasssurf
subroutine suews_cal_main(AerodynamicResistanceMethod, AH_MIN, AHProf_24hr, AH_SLOPE_Cooling, AH_SLOPE_Heating, alb, AlbMax_DecTr, AlbMax_EveTr, AlbMax_Grass, AlbMin_DecTr, AlbMin_EveTr, AlbMin_Grass, alpha_bioCO2, alpha_enh_bioCO2, alt, avkdn, avRh, avU1, BaseT, BaseTe, BaseTHDD, beta_bioCO2, beta_enh_bioCO2, bldgH, CapMax_dec, CapMin_dec, chAnOHM, CO2PointSource, cpAnOHM, CRWmax, CRWmin, DayWat, DayWatPer, DecTreeH, Diagnose, DiagQN, DiagQS, DRAINRT, dt_since_start, dqndt, qn1_av, dqnsdt, qn1_s_av, EF_umolCO2perJ, emis, EmissionsMethod, EnEF_v_Jkm, endDLS, EveTreeH, FAIBldg, FAIDecTree, FAIEveTree, Faut, FcEF_v_kgkm, fcld_obs, FlowChange, FrFossilFuel_Heat, FrFossilFuel_NonHeat, G1, G2, G3, G4, G5, G6, GDD_id, GDDFull, Gridiv, gsModel, HDD_id, HumActivity_24hr, IceFrac, id, Ie_a, Ie_end, Ie_m, Ie_start, imin, InternalWaterUse_h, IrrFracConif, IrrFracDecid, IrrFracGrass, isec, it, EvapMethod, iy, kkAnOHM, Kmax, LAI_id, LAICalcYes, LAIMax, LAIMin, LAI_obs, LAIPower, LAIType, lat, lenDay_id, ldown_obs, lng, MaxConductance, MaxFCMetab, MaxQFMetab, SnowWater, MetForcingData_grid, MinFCMetab, MinQFMetab, min_res_bioCO2, NARP_EMIS_SNOW, NARP_TRANS_SITE, NetRadiationMethod, OHM_coef, OHMIncQF, OHM_threshSW, OHM_threshWD, PipeCapacity, PopDensDaytime, PopDensNighttime, PopProf_24hr, PorMax_dec, PorMin_dec, Precip, PrecipLimit, PrecipLimitAlb, Press_hPa, QF0_BEU, Qf_A, Qf_B, Qf_C, qn1_obs, qs_obs, qf_obs, RadMeltFact, RAINCOVER, RainMaxRes, resp_a, resp_b, RoughLenHeatMethod, RoughLenMomMethod, RunoffToWater, S1, S2, SatHydraulicConduct, SDDFull, SDD_id, sfr, SMDMethod, SnowAlb, SnowAlbMax, SnowAlbMin, SnowPackLimit, SnowDens, SnowDensMax, SnowDensMin, SnowfallCum, SnowFrac, SnowLimBldg, SnowLimPaved, snowFrac_obs, SnowPack, SnowProf_24hr, snowUse, SoilDepth, soilstore_id, SoilStoreCap, StabilityMethod, startDLS, state_id, StateLimit, StorageHeatMethod, StoreDrainPrm, SurfaceArea, Tair_av, tau_a, tau_f, tau_r, Tmax_id, Tmin_id, T_CRITIC_Cooling, T_CRITIC_Heating, Temp_C, TempMeltFact, TH, theta_bioCO2, timezone, TL, TrafficRate, TrafficUnits, TraffProf_24hr, Ts5mindata_ir, tstep, tstep_prev, veg_type, WaterDist, WaterUseMethod, WetThresh, wu_m3, WUDay_id, DecidCap_id, albDecTr_id, albEveTr_id, albGrass_id, porosity_id, WUProfA_24hr, WUProfM_24hr, xsmd, Z, z0m_in, zdm_in, datetimeLine, dataOutLineSUEWS, dataOutLineSnow, dataOutLineESTM, dataoutLineRSL, DailyStateLine)
real(kind(1d0)) function, dimension(nsurf) update_snow_dens(tstep, SnowPack_prev, SnowDens_prev, tau_r, SnowDensMax, SnowDensMin)
integer, parameter ivconif
subroutine suews_cal_multitsteps(MetForcingBlock, len_sim, AerodynamicResistanceMethod, AH_MIN, AHProf_24hr, AH_SLOPE_Cooling, AH_SLOPE_Heating, alb, AlbMax_DecTr, AlbMax_EveTr, AlbMax_Grass, AlbMin_DecTr, AlbMin_EveTr, AlbMin_Grass, alpha_bioCO2, alpha_enh_bioCO2, alt, BaseT, BaseTe, BaseTHDD, beta_bioCO2, beta_enh_bioCO2, bldgH, CapMax_dec, CapMin_dec, chAnOHM, CO2PointSource, cpAnOHM, CRWmax, CRWmin, DayWat, DayWatPer, DecTreeH, Diagnose, DiagQN, DiagQS, DRAINRT, dt_since_start, dqndt, qn1_av, dqnsdt, qn1_s_av, EF_umolCO2perJ, emis, EmissionsMethod, EnEF_v_Jkm, endDLS, EveTreeH, FAIBldg, FAIDecTree, FAIEveTree, Faut, FcEF_v_kgkm, FlowChange, FrFossilFuel_Heat, FrFossilFuel_NonHeat, G1, G2, G3, G4, G5, G6, GDD_id, GDDFull, Gridiv, gsModel, HDD_id, HumActivity_24hr, IceFrac, Ie_a, Ie_end, Ie_m, Ie_start, InternalWaterUse_h, IrrFracConif, IrrFracDecid, IrrFracGrass, EvapMethod, kkAnOHM, Kmax, LAI_id, LAICalcYes, LAIMax, LAIMin, LAIPower, LAIType, lat, lng, MaxConductance, MaxFCMetab, MaxQFMetab, SnowWater, MinFCMetab, MinQFMetab, min_res_bioCO2, NARP_EMIS_SNOW, NARP_TRANS_SITE, NetRadiationMethod, OHM_coef, OHMIncQF, OHM_threshSW, OHM_threshWD, PipeCapacity, PopDensDaytime, PopDensNighttime, PopProf_24hr, PorMax_dec, PorMin_dec, PrecipLimit, PrecipLimitAlb, QF0_BEU, Qf_A, Qf_B, Qf_C, RadMeltFact, RAINCOVER, RainMaxRes, resp_a, resp_b, RoughLenHeatMethod, RoughLenMomMethod, RunoffToWater, S1, S2, SatHydraulicConduct, SDDFull, SDD_id, sfr, SMDMethod, SnowAlb, SnowAlbMax, SnowAlbMin, SnowPackLimit, SnowDens, SnowDensMax, SnowDensMin, SnowfallCum, SnowFrac, SnowLimBldg, SnowLimPaved, SnowPack, SnowProf_24hr, snowUse, SoilDepth, soilstore_id, SoilStoreCap, StabilityMethod, startDLS, state_id, StateLimit, StorageHeatMethod, StoreDrainPrm, SurfaceArea, Tair_av, tau_a, tau_f, tau_r, T_CRITIC_Cooling, T_CRITIC_Heating, TempMeltFact, TH, theta_bioCO2, timezone, TL, TrafficRate, TrafficUnits, Tmin_id, Tmax_id, lenday_id, TraffProf_24hr, Ts5mindata_ir, tstep, tstep_prev, veg_type, WaterDist, WaterUseMethod, WetThresh, WUDay_id, DecidCap_id, albDecTr_id, albEveTr_id, albGrass_id, porosity_id, WUProfA_24hr, WUProfM_24hr, Z, z0m_in, zdm_in, dataOutBlockSUEWS, dataOutBlockSnow, dataOutBlockESTM, dataOutBlockRSL, DailyStateBlock)
subroutine output_name_n(i, name, group, aggreg)
type(varattr), dimension(500) varlistall
real(kind(1d0)) function rh2qa(RH_dec, pres_hPa, Ta_degC)
subroutine suews_cal_horizontalsoilwater(sfr, SoilStoreCap, SoilDepth, SatHydraulicConduct, SurfaceArea, NonWaterFraction, tstep_real, soilstore_id, runoffSoil, runoffSoil_per_tstep)
subroutine suews_cal_dectime(id, it, imin, isec, dectime)
integer, parameter nvegsurf
subroutine output_size(nVar)
subroutine evap_suews(EvapMethod, state_is, WetThresh_is, capStore_is, vpd_hPa, avdens, avcp, qn_e, s_hPa, psyc_hPa, ResistSurf, RA, rb, tlv, rss, ev, qe)
subroutine suews_cal_weekday(iy, id, lat, dayofWeek_id)
subroutine narp(nsurf, sfr, SnowFrac, alb, emis, IceFrac, NARP_TRANS_SITE, NARP_EMIS_SNOW, DTIME, ZENITH_deg, tsurf_0, kdown, Temp_C, RH, Press_hPa, qn1_obs, SnowAlb, AlbedoChoice, ldown_option, NetRadiationMethod_use, DiagQN, QSTARall, QSTAR_SF, QSTAR_S, kclear, KUPall, LDOWN, LUPall, fcld, TSURFall, qn1_ind_snow, kup_ind_snow, Tsurf_ind_snow, Tsurf_ind, alb0, alb1)
subroutine suews_cal_wateruse(nsh_real, wu_m3, SurfaceArea, sfr, IrrFracConif, IrrFracDecid, IrrFracGrass, DayofWeek_id, WUProfA_24hr, WUProfM_24hr, InternalWaterUse_h, HDD_id, WUDay_id, WaterUseMethod, NSH, it, imin, DLS, wu_EveTr, wu_DecTr, wu_Grass, int_wu, ext_wu)
subroutine suews_cal_qe(Diagnose, snowuse, tstep, imin, it, EvapMethod, snowCalcSwitch, dayofWeek_id, CRWmin, CRWmax, dectime, lvS_J_kg, avRh, Press_hPa, Temp_C, RAsnow, psyc_hPa, sIce_hPa, PervFraction, vegfraction, addimpervious, qn1_snowfree, qf, qs, vpd_hPa, s_hPa, ResistSurf, RA, rb, snowdensmin, precip, PipeCapacity, RunoffToWater, NonWaterFraction, wu_EveTr, wu_DecTr, wu_Grass, addVeg, addWaterBody, SnowLimPaved, SnowLimBldg, SurfaceArea, FlowChange, drain, WetThresh, stateOld, mw_ind, SoilStoreCap, rainonsnow, freezmelt, freezstate, freezstatevol, Qm_Melt, Qm_rain, Tsurf_ind, sfr, StateLimit, AddWater, addwaterrunoff, StoreDrainPrm, SnowPackLimit, SnowProf_24hr, SnowPack_in, SnowFrac_in, SnowWater_in, iceFrac_in, SnowDens_in, runoff_per_interval_in, state_id_in, soilstore_id_in, state_id_out, soilstore_id_out, SnowPack_out, SnowFrac_out, SnowWater_out, iceFrac_out, SnowDens_out, runoffSoil, SnowRemoval, state_per_tstep, NWstate_per_tstep, qe, swe, chSnow_per_interval, ev_per_tstep, runoff_per_tstep, surf_chang_per_tstep, runoffPipes, mwstore, runoffwaterbody, runoffAGveg, runoffAGimpervious)
integer, parameter ncolumnsdataoutsuews
subroutine suews_cal_qs(StorageHeatMethod, qs_obs, OHMIncQF, Gridiv, id, tstep, dt_since_start, Diagnose, sfr, OHM_coef, OHM_threshSW, OHM_threshWD, soilstore_id, SoilStoreCap, state_id, SnowUse, SnowFrac, DiagQS, HDD_id, MetForcingData_grid, Ts5mindata_ir, qf, qn1, avkdn, avu1, temp_c, zenith_deg, avrh, press_hpa, ldown, bldgh, alb, emis, cpAnOHM, kkAnOHM, chAnOHM, EmissionsMethod, Tair_av, qn1_av_prev, dqndt_prev, qn1_s_av_prev, dqnsdt_prev, StoreDrainPrm, qn1_S, dataOutLineESTM, qs, qn1_av_next, dqndt_next, qn1_s_av_next, dqnsdt_next, deltaQi, a1, a2, a3)
real(kind(1d0)) function cal_tsfc(qh, avdens, avcp, RA, temp_c)
subroutine suews_cal_roughnessparameters(RoughLenMomMethod, sfr, bldgH, EveTreeH, DecTreeH, porosity_id, FAIBldg, FAIEveTree, FAIDecTree, z0m_in, zdm_in, Z, planF, Zh, z0m, zdm, ZZD)
real(kind(1d0)) function stab_psi_mom(StabilityMethod, ZL)
subroutine aerodynamicresistance(ZZD, z0m, AVU1, L_mod, UStar, VegFraction, AerodynamicResistanceMethod, StabilityMethod, RoughLenHeatMethod, RA)
integer, parameter decidsurf
integer, parameter pavsurf
real(kind(1d0)) avdens
real(kind(1d0)) function cal_tair_av(tair_av_prev, dt_since_start, tstep, temp_c)
subroutine errorhint(errh, ProblemFile, VALUE, value2, valueI)
subroutine suews_cal_surf(sfr, vegfraction, ImpervFraction, PervFraction, NonWaterFraction)
subroutine suews_cal_dls(id, startDLS, endDLS, DLS)
integer, parameter bldgsurf
integer, parameter watersurf
real(kind(1d0)) function qsatf(T, PMB)
subroutine drainage(is, state_is, StorCap, DrainEq, DrainCoef1, DrainCoef2, nsh_real, drain_is)
real function square_real(x)
integer, parameter ivdecid
subroutine anthropogenicemissions(CO2PointSource, EmissionsMethod, id, it, imin, DLS, nsh, DayofWeek_id, EF_umolCO2perJ, FcEF_v_kgkm, EnEF_v_Jkm, TrafficUnits, FrFossilFuel_Heat, FrFossilFuel_NonHeat, MinFCMetab, MaxFCMetab, MinQFMetab, MaxQFMetab, PopDensDaytime, PopDensNighttime, Temp_C, HDD_id, Qf_A, Qf_B, Qf_C, AH_MIN, AH_SLOPE_Heating, AH_SLOPE_Cooling, T_CRITIC_Heating, T_CRITIC_Cooling, TrafficRate, QF0_BEU, QF_SAHP, Fc_anthro, Fc_metab, Fc_traff, Fc_build, Fc_point, AHProf_24hr, HumActivity_24hr, TraffProf_24hr, PopProf_24hr, SurfaceArea)