SUEWS API Site
Documentation of SUEWS source code
suews_phys_ohm.f95
Go to the documentation of this file.
1 ! MODULE OHM_module
2 ! ! USE allocateArray
3 ! ! USE data_in
4 ! ! USE defaultNotUsed
5 ! ! USE gis_data
6 ! ! USE sues_data
7 ! ! USE time
8 !
9 ! IMPLICIT NONE
10 ! CONTAINS
11 !========================================================================================
12 SUBROUTINE ohm(qn1, qn1_av_prev, dqndt_prev, qn1_av_next, dqndt_next, &
13  qn1_S, qn1_s_av_prev, dqnsdt_prev, qn1_s_av_next, dqnsdt_next, &
14  tstep, dt_since_start, &
15  sfr, nsurf, &
16  Tair_mav_5d, &
17  OHM_coef, &
18  OHM_threshSW, OHM_threshWD, &
19  soilstore_id, SoilStoreCap, state_id, &
20  BldgSurf, WaterSurf, &
21  SnowUse, SnowFrac, &
22  DiagQS, &
23  a1, a2, a3, qs, deltaQi)
24  ! Made by HCW Jan 2015 to replace OHMnew (no longer needed).
25  ! Calculates net storage heat flux (QS) from Eq 4, Grimmond et al. 1991, Atm Env.
26  ! Accounts for variable timesteps in dQ*/dt term.
27  ! BareSoilSurfFraction removed so bare soil is now handled like other surfaces.
28  ! Snow part changed from summer wet to winter wet coefficients.
29  ! Changed -333 checks to -999 checks and added error handling
30  ! Gradient now calculated for t-1 (was previously calculated for t-2).
31  ! TS 28 Jun 2018:
32  ! improved and tested the phase-in method for calculating dqndt
33  ! TS & SG 30 Apr 2018:
34  ! a new calculation scheme of dqndt by using a phase-in approach that releases
35  ! the requirement for storeing multiple qn values for adapting SUEWS into WRF
36  ! TS 07 Aug 2017:
37  ! 1. interface changed to account for explict passing
38  ! 2. calculation refactorization.
39  ! HCW 25 Feb 2015:
40  ! Adapted q1,q2,q3 & r1,r2,r3 for multiple grids
41  ! HCW 14 Dec 2016:
42  ! Thresholds for Summer/Winter and Wet/Dry now provided in input files
43  ! Calculation of dqndt now uses hourly running mean rather than instantaneous values
44  ! To Do:
45  ! - No canyons implemented at the moment [OHM_coef(nsurf+1,,)]
46  !========================================================================================
47 
48  IMPLICIT NONE
49  INTEGER, INTENT(in) :: tstep ! time step [s]
50  INTEGER, INTENT(in) :: dt_since_start ! time since simulation starts [s]
51 
52  REAL(KIND(1d0)), INTENT(in)::qn1 ! net all-wave radiation
53  REAL(KIND(1d0)), INTENT(in)::qn1_S ! net all-wave radiation over snow
54  REAL(KIND(1d0)), INTENT(in)::sfr(nsurf) ! surface fractions
55  REAL(KIND(1d0)), INTENT(in)::SnowFrac(nsurf) ! snow fractions of each surface
56  REAL(KIND(1d0)), INTENT(in)::Tair_mav_5d ! Tair_mav_5d=HDD(id-1,4) HDD at the begining of today (id-1)
57  REAL(KIND(1d0)), INTENT(in)::OHM_coef(nsurf + 1, 4, 3) ! OHM coefficients
58  REAL(KIND(1d0)), INTENT(in)::OHM_threshSW(nsurf + 1), OHM_threshWD(nsurf + 1) ! OHM thresholds
59  REAL(KIND(1d0)), INTENT(in)::soilstore_id(nsurf) ! soil moisture
60  REAL(KIND(1d0)), INTENT(in)::SoilStoreCap(nsurf) ! capacity of soil store
61  REAL(KIND(1d0)), INTENT(in)::state_id(nsurf) ! wetness status
62 
63  INTEGER, INTENT(in)::nsurf ! number of surfaces
64  ! INTEGER,INTENT(in)::nsh ! number of timesteps in one hour
65  ! integer,intent(in) :: dt ! current timestep [second]
66  ! integer,intent(INOUT) :: dt0 ! period length for qn1 memory
67  INTEGER, INTENT(in)::BldgSurf ! code for specific surfaces
68  INTEGER, INTENT(in)::WaterSurf ! code for specific surfaces
69  INTEGER, INTENT(in)::SnowUse ! option for snow related calculations
70  INTEGER, INTENT(in)::DiagQS ! diagnostic option
71 
72  REAL(KIND(1d0)), INTENT(in)::qn1_av_prev
73  REAL(KIND(1d0)), INTENT(out)::qn1_av_next
74  REAL(KIND(1d0)), INTENT(in)::dqndt_prev ! Rate of change of net radiation [W m-2 h-1] at t-1
75  REAL(KIND(1d0)), INTENT(out)::dqndt_next ! Rate of change of net radiation [W m-2 h-1] at t-1
76  REAL(KIND(1d0)), INTENT(in)::qn1_s_av_prev
77  REAL(KIND(1d0)), INTENT(out)::qn1_s_av_next
78  REAL(KIND(1d0)), INTENT(in)::dqnsdt_prev ! Rate of change of net radiation [W m-2 h-1] at t-1
79  REAL(KIND(1d0)), INTENT(out)::dqnsdt_next ! Rate of change of net radiation [W m-2 h-1] at t-1
80 
81  ! REAL(KIND(1d0)), INTENT(inout)::qn1_av
82  ! REAL(KIND(1d0)), INTENT(inout)::dqndt !Rate of change of net radiation [W m-2 h-1] at t-1
83  ! REAL(KIND(1d0)), INTENT(inout)::qn1_s_av
84  ! REAL(KIND(1d0)), INTENT(inout)::dqnsdt !Rate of change of net radiation [W m-2 h-1] at t-1
85  ! REAL(KIND(1d0)),INTENT(inout)::qn1_store_grid(nsh)
86  ! REAL(KIND(1d0)),INTENT(inout)::qn1_av_store_grid(2*nsh+1)
87  ! REAL(KIND(1d0)),INTENT(inout)::qn1_S_store_grid(nsh)
88  ! REAL(KIND(1d0)),INTENT(inout)::qn1_S_av_store_grid(2*nsh+1)
89 
90  REAL(KIND(1d0)), INTENT(out):: qs ! storage heat flux
91  ! REAL(KIND(1d0)),INTENT(out)::deltaQi(nsurf+1) ! storage heat flux of snow surfaces
92  REAL(KIND(1d0)), INTENT(out)::deltaQi(nsurf) ! storage heat flux of snow surfaces
93 
94  REAL(KIND(1d0)), INTENT(out):: a1, a2, a3 ! OHM coefficients of grid
95 
96  ! REAL(KIND(1d0)):: nsh_nna ! number of timesteps per hour with non -999 values (used for spinup)
97 
98  ! REAL(KIND(1d0)):: dqndt !Rate of change of net radiation [W m-2 h-1] at t-1
99  ! REAL(KIND(1d0)):: surfrac !Surface fraction accounting for SnowFrac if appropriate
100 
101  ! REAL(KIND(1d0)):: qn1_av, qn1_S_av !Average net radiation over previous hour [W m-2]
102  REAL(KIND(1d0)):: deltaQi0 ! temporarily store
103 
104  ! REAL(KIND(1d0)):: qn1_store_grid0(nsh), qn1_av_store_grid0(2*nsh+1) ! temporarily store
105 
106  !These are now provided in SiteInfo (OHMthresh for Summer/Winter and Wet/Dry)
107  !!real(kind(1d0)):: OHM_TForSummer = 5 !Use summer coefficients if 5-day Tair >= 5 degC
108  !real(kind(1d0)):: OHM_TForSummer = 10 !Use summer coefficients if 5-day Tair >= 10 degC - modified for UK HCW 14 Dec 2015
109  !real(kind(1d0)):: OHM_SMForWet = 0.9 !Use wet coefficients if SM close to soil capacity
110 
111  CALL ohm_coef_cal(sfr, nsurf, &
112  tair_mav_5d, ohm_coef, ohm_threshsw, ohm_threshwd, &
113  soilstore_id, soilstorecap, state_id, &
114  bldgsurf, watersurf, &
115  snowuse, snowfrac, &
116  a1, a2, a3)
117  ! WRITE(*,*) '----- OHM coeffs new-----'
118  ! WRITE(*,*) a1,a2,a3
119 
120  ! Old OHM calculations (up to v2016a)
121  !! Calculate radiation part ------------------------------------------------------------
122  !qs=NAN !qs = Net storage heat flux [W m-2]
123  !if(qn1>-999) then !qn1 = Net all-wave radiation [W m-2]
124  ! !if(q1>-999.and.q3>-999) then
125  ! !dqndt = 0.5*(q3-q1)*nsh_real !gradient at t-2
126  ! dqndt = 0.5*(qn1-q2_grids(Gridiv))*nsh_real !gradient at t-1
127  !
128  ! !Calculate net storage heat flux
129  ! qs = qn1*a1 + dqndt*a2 + a3 !Eq 4, Grimmond et al. 1991
130  ! !endif
131  ! !q1=q2 !q1 = net radiation at t-2 (at t-3 when q1 used in next timestep)
132  ! !q2=q3 !q2 = net radiation at t-1
133  ! !q3=qn1 !q3 = net radiation at t (at t-1 when q3 used in next timestep)
134  ! q1_grids(Gridiv) = q2_grids(Gridiv) !q1 = net radiation at t-2 (at t-3 when q1 used in next timestep)
135  ! q2_grids(Gridiv) = q3_grids(Gridiv) !q2 = net radiation at t-1
136  ! q3_grids(Gridiv) = qn1 !q3 = net radiation at t (at t-1 when q3 used in next timestep)
137  !else
138  ! call ErrorHint(21,'Bad value for qn1 found during OHM calculation',qn1,NotUsed,notUsedI)
139  !endif
140 
141  ! New OHM calculations (v2017a onwards) using running mean (HCW Dec 2016)
142  ! Calculate radiation part ------------------------------------------------------------
143  qs = -999 !qs = Net storage heat flux [W m-2]
144  IF (qn1 > -999) THEN !qn1 = Net all-wave radiation [W m-2]
145  ! Store instantaneous qn1 values for previous hour (qn1_store_grid) and average (qn1_av)
146  ! print*,''
147  ! CALL OHM_dqndt_cal(nsh,qn1,qn1_store_grid,qn1_av_store_grid,dqndt)
148  ! print*, 'old dqndt',dqndt
149  CALL ohm_dqndt_cal_x(tstep, dt_since_start, qn1_av_prev, qn1, dqndt_prev, &
150  qn1_av_next, dqndt_next)
151  ! print*, 'new dqndt',dqndt
152 
153  ! Calculate net storage heat flux
154  CALL ohm_qs_cal(qn1, dqndt_next, a1, a2, a3, qs)
155  IF (diagqs == 1) WRITE (*, *) 'qs: ', qs, 'qn1:', qn1, 'dqndt: ', dqndt_next
156 
157  ELSE
158  CALL errorhint(21, 'In SUEWS_OHM.f95: bad value for qn1 found during qs calculation.', qn1, -55.55d0, -55)
159  ENDIF
160 
161  !write(*,*) qs
162  !write(*,*) '--------------------'
163 
164  ! Do snow calculations separately -----
165  ! Added by LJ in August 2013
166  IF (snowuse == 1) THEN
167  deltaqi = -999
168  IF (qn1_s > -999) THEN
169  ! Old OHM calculations (commented out HCW Dec 2016)
170  !!if(r1>-999.and.r3>-999) then
171  ! !dqndt = 0.5*(r3-r1)*nsh_real !gradient at t-2
172  ! dqndt = 0.5*(qn1_S-r2_grids(Gridiv))*nsh_real !gradient at t-1
173  ! ! Calculate net storage heat flux for snow surface (winter wet conditions HCW 15/01/2015)
174  ! deltaQi = qn1_S*OHM_coef(nsurf+1,3,1) + dqndt*OHM_coef(nsurf+1,3,2) + OHM_coef(nsurf+1,3,3)
175  !!endif
176  !r1_grids(Gridiv)=r2_grids(Gridiv)
177  !r2_grids(Gridiv)=r3_grids(Gridiv)
178  !r3_grids(Gridiv)=qn1_S
179  ! New OHM calculations
180  ! Store instantaneous qn1 values for previous hour (qn1_store_grid) and average (qn1_av)
181  ! CALL OHM_dqndt_cal(nsh,qn1_S,qn1_S_store_grid,qn1_S_av_store_grid,dqndt)
182 
183  CALL ohm_dqndt_cal_x(tstep, dt_since_start, qn1_s_av_prev, qn1_s, dqnsdt_prev, &
184  qn1_s_av_next, dqnsdt_next)
185 
186  ! Calculate net storage heat flux for snow surface (winter wet conditions)
187  CALL ohm_qs_cal(qn1_s, dqnsdt_next, &
188  ohm_coef(nsurf + 1, 3, 1), ohm_coef(nsurf + 1, 3, 2), ohm_coef(nsurf + 1, 3, 3), &
189  deltaqi0)
190  deltaqi = deltaqi0
191 
192  ELSE
193  CALL errorhint(21, 'In SUEWS_OHM.f95: bad value for qn1(snow) found during qs calculation.', qn1_s, -55.55d0, -55)
194  ENDIF
195 
196  ENDIF
197 
198  RETURN
199 END SUBROUTINE ohm
200 !========================================================================================
201 
202 SUBROUTINE ohm_coef_cal(sfr, nsurf, &
203  Tair_mav_5d, OHM_coef, OHM_threshSW, OHM_threshWD, &
204  soilstore_id, SoilStoreCap, state_id, &
205  BldgSurf, WaterSurf, &
206  SnowUse, SnowFrac, &
207  a1, a2, a3)
208  IMPLICIT NONE
209  INTEGER, INTENT(in) :: &
210  nsurf, & ! number of surfaces
211  SnowUse, & ! option for snow related calculations
212  BldgSurf, WaterSurf ! code for specific surfaces
213  REAL(KIND(1d0)), INTENT(in) :: &
214  sfr(nsurf), & ! surface cover fractions
215  SnowFrac(nsurf), & ! snow fractions of each surface
216  Tair_mav_5d, & ! Tair_mav_5d=HDD(id-1,4) HDD at the begining of today (id-1)
217  OHM_coef(nsurf + 1, 4, 3), &
218  OHM_threshSW(nsurf + 1), OHM_threshWD(nsurf + 1), & ! OHM thresholds
219  soilstore_id(nsurf), & ! soil moisture
220  SoilStoreCap(nsurf), &! capacity of soil store
221  state_id(nsurf) ! wetness status
222  REAL(KIND(1d0)), INTENT(out):: a1, a2, a3
223 
224  REAL(KIND(1d0)) :: surfrac
225  INTEGER :: i, ii, is
226 
227  ! OHM coefficients --------
228  ! Set to zero initially
229  a1 = 0 ![-]
230  a2 = 0 ![h]
231  a3 = 0 ![W m-2]
232  ! -------------------------
233 
234  ! Loop through surface types ----------------------------------------------------------
235  DO is = 1, nsurf
236  surfrac = sfr(is)
237 
238  ! Use 5-day running mean Tair to decide whether it is summer or winter ----------------
239  IF (tair_mav_5d >= ohm_threshsw(is)) THEN !Summer
240  ii = 0
241  ELSE !Winter
242  ii = 2
243  ENDIF
244 
245  IF (state_id(is) > 0) THEN !Wet surface
246  i = ii + 1
247  ELSE !Dry surface
248  i = ii + 2
249  ! If the surface is dry but SM is close to capacity, use coefficients for wet surfaces
250  IF (is > bldgsurf .AND. is /= watersurf) THEN !Wet soil (i.e. EveTr, DecTr, Grass, BSoil surfaces)
251  IF (soilstore_id(is)/soilstorecap(is) > ohm_threshwd(is)) THEN
252  i = ii + 1
253  ENDIF
254  ENDIF
255  ENDIF
256 
257  ! If snow, adjust surface fractions accordingly
258  IF (snowuse == 1 .AND. is /= bldgsurf .AND. is /= watersurf) THEN ! QUESTION: Why is BldgSurf excluded here?
259  surfrac = surfrac*(1 - snowfrac(is))
260  ENDIF
261 
262  ! Calculate the areally-weighted OHM coefficients
263  a1 = a1 + surfrac*ohm_coef(is, i, 1)
264  a2 = a2 + surfrac*ohm_coef(is, i, 2)
265  a3 = a3 + surfrac*ohm_coef(is, i, 3)
266 
267  ENDDO !end of loop over surface types ------------------------------------------------
268 END SUBROUTINE ohm_coef_cal
269 
270 ! Updated OHM calculations for WRF-SUEWS coupling (v2018b onwards) weighted mean (TS Apr 2018)
271 SUBROUTINE ohm_dqndt_cal_x(dt, dt_since_start, qn1_av_prev, qn1, dqndt_prev, qn1_av_next, dqndt_next)
272  IMPLICIT NONE
273  INTEGER, INTENT(in) :: dt ! time step [s]
274  INTEGER, INTENT(in) :: dt_since_start ! time since simulation starts [s]
275  REAL(KIND(1d0)), INTENT(in) :: qn1 ! new qn1 value [W m-2]
276  REAL(KIND(1d0)), INTENT(in) :: qn1_av_prev ! weighted average of qn1 [W m-2]
277  REAL(KIND(1d0)), INTENT(in) :: dqndt_prev ! dQ* per dt for 60 min [W m-2 h-1]
278  REAL(KIND(1d0)), INTENT(out) :: qn1_av_next ! weighted average of qn1 [W m-2]
279  REAL(KIND(1d0)), INTENT(out) :: dqndt_next ! dQ* per dt for 60 min [W m-2 h-1]
280  REAL(KIND(1d0)), PARAMETER :: dt0_thresh = 3600 ! threshold for period of dqndt0 [s]
281  REAL(KIND(1d0)), PARAMETER :: window_hr = 2 ! window size for Difference calculation [hr]
282 
283  INTEGER :: dt0 ! period of dqndt0 [s]
284 
285  REAL(KIND(1d0)) :: qn1_av_0 !, qn1_av_start,qn1_av_end
286 
287  ! if previous period shorter than dt0_thresh, expand the storage/memory period
288  IF (dt_since_start < dt0_thresh) THEN ! spinup period
289  dt0 = dt_since_start + dt
290 
291  ELSE ! effective period
292  dt0 = dt0_thresh
293  ENDIF
294 
295  ! get weighted average at a previous time specified by `window_hr`
296  qn1_av_0 = qn1_av_prev - dqndt_prev*(window_hr - dt/3600.)
297 
298  ! averaged qn1 for previous period = dt0_thresh
299  qn1_av_next = (qn1_av_prev*(dt0 - dt) + qn1*dt)/(dt0)
300 
301  ! do weighted average to calculate the difference by using the memory value and new forcing value
302  ! NB: keep the output dqndt in [W m-2 h-1]
303  dqndt_next = (qn1_av_next - qn1_av_0)/window_hr
304 
305 END SUBROUTINE ohm_dqndt_cal_x
306 
307 ! New OHM calculations (v2017a-v2018a) using running mean (HCW Dec 2016)
308 SUBROUTINE ohm_dqndt_cal(nsh, qn1, qn1_store_grid, qn1_av_store_grid, dqndt)
309  IMPLICIT NONE
310  INTEGER, INTENT(in) :: nsh ! number of timesteps in one hour
311  REAL(KIND(1d0)), INTENT(in) :: qn1
312  REAL(KIND(1d0)), INTENT(inout) :: qn1_store_grid(nsh) ! instantaneous qn1 values for previous hour
313  REAL(KIND(1d0)), INTENT(inout) :: qn1_av_store_grid(2*nsh + 1) ! average qn1 values for previous hour
314  REAL(KIND(1d0)), INTENT(out) :: dqndt !dQ* per dt for 60 min
315 
316  REAL(KIND(1d0)) :: qn1_av
317  INTEGER :: nsh_nna
318 
319  dqndt = -999 ! initialise as -999
320 
321  ! Store instantaneous qn1 values for previous hour (qn1_store_grid) and average (qn1_av)
322  IF (nsh > 1) THEN
323  qn1_store_grid = cshift(qn1_store_grid, 1) ! shift to left with one place
324  qn1_store_grid(nsh) = qn1
325  nsh_nna = count(qn1_store_grid /= -999, dim=1) !Find how many are not -999s !bug fixed HCW 08 Feb 2017
326  qn1_av = sum(qn1_store_grid, mask=qn1_store_grid /= -999)/nsh_nna
327  ELSEIF (nsh == 1) THEN
328  qn1_store_grid(:) = qn1
329  qn1_av = qn1
330  ENDIF
331  ! Store hourly average values (calculated every timestep) for previous 2 hours
332  IF (nsh > 1) THEN
333  qn1_av_store_grid = cshift(qn1_av_store_grid, 1)
334  qn1_av_store_grid(2*nsh + 1) = qn1_av
335  ELSEIF (nsh == 1) THEN
336  qn1_av_store_grid(:) = qn1_av
337  ENDIF
338  ! Calculate dQ* per dt for 60 min (using running mean Q* at t hours and (t-2) hours)
339  IF (any(qn1_av_store_grid == -999)) THEN
340  dqndt = 0 ! Set dqndt term to zero for spinup
341  ELSE
342  dqndt = 0.5*(qn1_av_store_grid((2*nsh + 1)) - qn1_av_store_grid(1))
343  ENDIF
344 
345 END SUBROUTINE ohm_dqndt_cal
346 
347 SUBROUTINE ohm_qs_cal(qn1, dqndt, a1, a2, a3, qs)
348  IMPLICIT NONE
349  REAL(KIND(1d0)), INTENT(in) :: qn1, dqndt, a1, a2, a3
350  REAL(KIND(1d0)), INTENT(out):: qs
351  qs = -999 ! initialise as -999
352  qs = qn1*a1 + dqndt*a2 + a3 !Eq 4, Grimmond et al. 1991
353 
354 END SUBROUTINE ohm_qs_cal
355 
356 ! END MODULE OHM_module
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 ohm_dqndt_cal(nsh, qn1, qn1_store_grid, qn1_av_store_grid, dqndt)
subroutine ohm_dqndt_cal_x(dt, dt_since_start, qn1_av_prev, qn1, dqndt_prev, qn1_av_next, dqndt_next)
subroutine errorhint(errh, ProblemFile, VALUE, value2, valueI)
subroutine ohm_coef_cal(sfr, nsurf, Tair_mav_5d, OHM_coef, OHM_threshSW, OHM_threshWD, soilstore_id, SoilStoreCap, state_id, BldgSurf, WaterSurf, SnowUse, SnowFrac, a1, a2, a3)
subroutine ohm_qs_cal(qn1, dqndt, a1, a2, a3, qs)