Skip to content

Commit cbd8ad4

Browse files
gnrgomesgnrgomes
committed
Correct RG readmeteo and arcsin and unit conversion
1 parent eee6447 commit cbd8ad4

Some content is hidden

Large Commits have some content hidden by default. Use the searchbox below for content that may be hidden.

48 files changed

+336
-75
lines changed

‎settings_tpl.xml‎

Lines changed: 140 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -142,6 +142,140 @@ You can use $(SettingsDir) or $(SettingsPath) to refer directory containing the
142142
<setoption name="useRelHumidityMaps" choice="1"/>
143143
</lfoptions>
144144

145+
<lfconversions>
146+
<comment>
147+
**************************************************************
148+
UNITS CONVERSION FACTORS FOR INPUT VARIABLES
149+
**************************************************************
150+
</comment>
151+
152+
<setconversion name="TMaxMaps" value="">
153+
<comment>
154+
Maximum daily temperature
155+
Units: [K] if "TemperatureInKelvinFlag" is on, else [C].
156+
</comment>
157+
</setconversion>
158+
159+
<setconversion name="TMinMaps" value="">
160+
<comment>
161+
Minimum daily temperature
162+
Units: [K] if "TemperatureInKelvinFlag" is on, else [C].
163+
</comment>
164+
</setconversion>
165+
166+
<setconversion name="TDewMaps" value="">
167+
<comment>
168+
TDew - Dew point Temperature used for actual vap pressure
169+
Units: [K] if "TemperatureInKelvinFlag" is on, else [C].
170+
Used when the "GLOFAS" option is switched on.
171+
172+
</comment>
173+
</setconversion>
174+
175+
<setconversion name="RelHMaps" value="">
176+
<comment>
177+
RelHMaps - Relative Humidity used for actual vap pressure
178+
Units: [%]
179+
Used when the "EFAS" option is switched on.
180+
181+
</comment>
182+
</setconversion>
183+
184+
<setconversion name="EActMaps" value="">
185+
<comment>
186+
Actual vapour pressure [mbar]
187+
Used when the "EFAS" option is switched on.
188+
If the "CORDEX" option is on, "QairMaps" is used instead.
189+
</comment>
190+
</setconversion>
191+
192+
<setconversion name="WindMaps" value="">
193+
<comment>
194+
Wind speed at 10 m from surface [m/s]
195+
</comment>
196+
</setconversion>
197+
198+
<setconversion name="WindMapsU" value="">
199+
<comment>
200+
Wind speed U component at 10 m from surface [m/s]
201+
</comment>
202+
</setconversion>
203+
204+
<setconversion name="WindMapsV" value="">
205+
<comment>
206+
Wind speed V component at 10 m from surface [m/s]
207+
</comment>
208+
</setconversion>
209+
210+
<setconversion name="TAvgMaps" value="">
211+
<comment>
212+
average daily temperature
213+
Units: [K] if "TemperatureInKelvinFlag" is on, else [C].
214+
</comment>
215+
</setconversion>
216+
217+
<setconversion name="PSurfMaps" value="">
218+
<comment>
219+
Instantaneous sea level pressure [Pa]
220+
</comment>
221+
</setconversion>
222+
223+
<setconversion name="QAirMaps" value="">
224+
<comment>
225+
2 m instantaneous specific humidity [kg/kg]
226+
Used when the "CORDEX" option is switched on.
227+
If the "EFAS" option is on, "EActMaps" is used instead.
228+
</comment>
229+
</setconversion>
230+
231+
<setconversion name="RgdMaps" value="">
232+
<comment>
233+
rgd - calculated solar radiation [J/m2]
234+
Used when the "EFAS" option is switched on.
235+
If the "CORDEX" option is on, the following are used instead: "RdsMaps", "RdlMaps", "RusMaps", "RulMaps".
236+
</comment>
237+
</setconversion>
238+
239+
<setconversion name="RdsMaps" value="86400">
240+
<comment>
241+
rds - Downward short wave radiation [J/m2]
242+
Used when the "CORDEX" option is switched on.
243+
If the "EFAS" option is on, "RgdMaps" is used instead.
244+
</comment>
245+
</setconversion>
246+
247+
<setconversion name="RdlMaps" value="86400">
248+
<comment>
249+
rdl - Down long wave radiation [J/m2]
250+
Used when the "CORDEX" option is switched on.
251+
If the "EFAS" option is on, "RgdMaps" is used instead.
252+
</comment>
253+
</setconversion>
254+
255+
<setconversion name="RusMaps" value="86400">
256+
<comment>
257+
rus - up short wave radiation [J/m2]
258+
Used when the "CORDEX" option is switched on.
259+
If the "EFAS" option is on, "RgdMaps" is used instead.
260+
</comment>
261+
</setconversion>
262+
263+
<setconversion name="RulMaps" value="86400">
264+
<comment>
265+
rul - up long wave radiation [J/m2]
266+
Used when the "CORDEX" option is switched on.
267+
If the "EFAS" option is on, "RgdMaps" is used instead.
268+
</comment>
269+
</setconversion>
270+
271+
<setconversion name="RnlMaps" value="">
272+
<comment>
273+
rnl - net long wave radiation [J/m2]
274+
Used when the "GLOFAS" option is switched on.
275+
</comment>
276+
</setconversion>
277+
</lfconversions>
278+
145279

146280
<lfbinding>
147281
<!--
@@ -325,7 +459,7 @@ You can use $(SettingsDir) or $(SettingsPath) to refer directory containing the
325459

326460
<textvar name="RgdMaps" value="$(PathMeteoIn)/EMO-5-rg_1990_2017">
327461
<comment>
328-
rgd - calculated solar radiation [J/day/m2]
462+
rgd - calculated solar radiation [J/m2]
329463
Used when the "EFAS" option is switched on.
330464
If the "CORDEX" option is on, the following are used instead: "RdsMaps", "RdlMaps", "RusMaps", "RulMaps".
331465
</comment>
@@ -354,39 +488,39 @@ You can use $(SettingsDir) or $(SettingsPath) to refer directory containing the
354488

355489
<textvar name="RdsMaps" value="$(PathMeteoIn)/rsds_1_15">
356490
<comment>
357-
rds - Downward short wave radiation [W/m2]
491+
rds - Downward short wave radiation [J/m2]
358492
Used when the "CORDEX" option is switched on.
359493
If the "EFAS" option is on, "RgdMaps" is used instead.
360494
</comment>
361495
</textvar>
362496

363497
<textvar name="RdlMaps" value="$(PathMeteoIn)/rlds_1_15">
364498
<comment>
365-
rdl - Down long wave radiation [W/m2]
499+
rdl - Down long wave radiation [J/m2]
366500
Used when the "CORDEX" option is switched on.
367501
If the "EFAS" option is on, "RgdMaps" is used instead.
368502
</comment>
369503
</textvar>
370504

371505
<textvar name="RusMaps" value="$(PathMeteoIn)/rsus_1_15">
372506
<comment>
373-
rus - up short wave radiation [W/m2]
507+
rus - up short wave radiation [J/m2]
374508
Used when the "CORDEX" option is switched on.
375509
If the "EFAS" option is on, "RgdMaps" is used instead.
376510
</comment>
377511
</textvar>
378512

379513
<textvar name="RulMaps" value="$(PathMeteoIn)/rlus_1_15">
380514
<comment>
381-
rul - up long wave radiation [W/m2]
515+
rul - up long wave radiation [J/m2]
382516
Used when the "CORDEX" option is switched on.
383517
If the "EFAS" option is on, "RgdMaps" is used instead.
384518
</comment>
385519
</textvar>
386520

387521
<textvar name="RnlMaps" value="$(PathMeteoIn)/rlus_1_15">
388522
<comment>
389-
rnl - net long wave radiation [W/m2]
523+
rnl - net long wave radiation [J/m2]
390524
Used when the "GLOFAS" option is switched on.
391525
</comment>
392526
</textvar>

‎settings_tpl.yaml‎

Lines changed: 6 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -58,7 +58,7 @@ prefixes-input-variables:
5858
PrefixEAct: pd
5959
# prefix wind speed maps at 10m [m/s]
6060
PrefixWind: ws
61-
# prefix down solar radiation maps [W/m2]
61+
# prefix down solar radiation maps [J/m2]
6262
PrefixRgd: rg
6363
# prefix average temperature maps [k]
6464
PrefixTAvg: ta
@@ -139,23 +139,23 @@ input-maps:
139139
# Used when the "CORDEX" option is switched on.
140140
# If the "EFAS" option is on, "EActMaps" is used instead.
141141
QAirMaps: $(PathMeteoIn)/huss_1_15
142-
# rds - Downward short wave radiation [W/m2]
142+
# rds - Downward short wave radiation [J/m2]
143143
# Used when the "CORDEX" option is switched on.
144144
# If the "EFAS" option is on, "RgdMaps" is used instead.
145145
RdsMaps: $(PathMeteoIn)/rsds_1_15
146-
# rdl - Down long wave radiation [W/m2]
146+
# rdl - Down long wave radiation [J/m2]
147147
# Used when the "CORDEX" option is switched on.
148148
# If the "EFAS" option is on, "RgdMaps" is used instead.
149149
RdlMaps: $(PathMeteoIn)/rlds_1_15
150-
# rus - up short wave radiation [W/m2]
150+
# rus - up short wave radiation [J/m2]
151151
# Used when the "CORDEX" option is switched on.
152152
# If the "EFAS" option is on, "RgdMaps" is used instead.
153153
RusMaps: $(PathMeteoIn)/rsus_1_15
154-
# rul - up long wave radiation [W/m2]
154+
# rul - up long wave radiation [J/m2]
155155
# Used when the "CORDEX" option is switched on.
156156
# If the "EFAS" option is on, "RgdMaps" is used instead.
157157
RulMaps: $(PathMeteoIn)/rlus_1_15
158-
# rnl - net long wave radiation [W/m2]
158+
# rnl - net long wave radiation [J/m2]
159159
# Used when the "GLOFAS" option is switched on.
160160
RnlMaps: $(PathMeteoIn)/rlus_1_15
161161

‎src/lisvap/hydrological/miscinitial.py‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -87,6 +87,7 @@ def initial(self):
8787
self.var.Qair = None
8888
self.var.Wind = None
8989

90+
self.var.RelH = None
9091
self.var.ESat = None
9192

9293
self.var.Rds = None

‎src/lisvap/hydrological/readmeteo.py‎

Lines changed: 8 additions & 9 deletions
Original file line numberDiff line numberDiff line change
@@ -113,13 +113,13 @@ def dynamic(self):
113113
if self.settings.get_option('CORDEX'):
114114
self.var.Psurf = readnetcdf(self.settings.binding['PSurfMaps'], self.var.currentTimeStep(), variable_binding='PSurfMaps', splitIO=self.splitIO)
115115
self.var.Qair = readnetcdf(self.settings.binding['QAirMaps'], self.var.currentTimeStep(), variable_binding='QAirMaps', splitIO=self.splitIO)
116-
# Downward short wave radiation [W/m2]
116+
# Downward short wave radiation [J/m2]
117117
self.var.Rds = readnetcdf(self.settings.binding['RdsMaps'], self.var.currentTimeStep(), variable_binding='RdsMaps', splitIO=self.splitIO)
118-
# Down long wave radiation [W/m2]
118+
# Down long wave radiation [J/m2]
119119
self.var.Rdl = readnetcdf(self.settings.binding['RdlMaps'], self.var.currentTimeStep(), variable_binding='RdlMaps', splitIO=self.splitIO)
120-
# upward short wave radiation [W/m2]
120+
# upward short wave radiation [J/m2]
121121
self.var.Rus = readnetcdf(self.settings.binding['RusMaps'], self.var.currentTimeStep(), variable_binding='RusMaps', splitIO=self.splitIO)
122-
# upward long wave radiation [W/m2]
122+
# upward long wave radiation [J/m2]
123123
self.var.Rul = readnetcdf(self.settings.binding['RulMaps'], self.var.currentTimeStep(), variable_binding='RulMaps', splitIO=self.splitIO)
124124
else: # EFAS or GLOFAS
125125
self.read_vapor_pressure()
@@ -132,13 +132,12 @@ def dynamic(self):
132132
if self.settings.get_option('GLOFAS'):
133133
# set of forcings (rg, rn, ta, td, wu, wv)
134134
# Net long wave radiation [J/m2/day]
135-
self.var.Rnl = readnetcdf(self.settings.binding['RNMaps'], self.var.currentTimeStep(), variable_binding='RNMaps', splitIO=self.splitIO) * -1
135+
self.var.Rnl = readnetcdf(self.settings.binding['RnlMaps'], self.var.currentTimeStep(), variable_binding='RnlMaps', splitIO=self.splitIO)
136+
# For GLOFAS we need to invert the signal of the files that come from ERA5
137+
# so it fits the formulas we are using in Lisvap
138+
self.var.Rnl = self.var.Rnl * -1
136139

137140
if self.settings.get_option('CORDEX'):
138-
self.var.Rds = self.var.Rds * 86400
139-
self.var.Rdl = self.var.Rdl * 86400
140-
self.var.Rus = self.var.Rus * 86400
141-
self.var.Rul = self.var.Rul * 86400
142141
self.var.EAct = (self.var.Psurf * self.var.Qair) / 62.2
143142
# [KPA] * [kg/kg] = KPa
144143

‎src/lisvap/lisvapdynamic.py‎

Lines changed: 36 additions & 11 deletions
Original file line numberDiff line numberDiff line change
@@ -19,8 +19,10 @@
1919

2020
import sys
2121
import datetime
22+
import cftime
2223

2324
import numpy as np
25+
from netCDF4 import num2date, date2num
2426

2527
from .utils.operators import exp, maximum, cos, sin, ifthenelse, asin, scalar, cover, tan, sqr, sqrt, abs
2628
from .utils import LisSettings, TimeProfiler, DynamicModel
@@ -34,6 +36,8 @@ def __init__(self):
3436
initialization of the hydrological modules
3537
"""
3638
# super(LisvapModelDyn, self).__init__()
39+
self.gregorian_calendars_id = [u'standard', u'gregorian', u'proleptic_gregorian']
40+
self.calendar_day_start_gregorian = None
3741
self.calendar_date = None
3842
self.calendar_day = None
3943
self.time_since_start = None
@@ -80,24 +84,23 @@ def angot_radiation(self):
8084
# solar constant at top of the atmosphere [J/m2/s]
8185
solar_constant = self.AvSolarConst * (1 + (0.033 * cos(360. * self.calendar_day / 365.)))
8286
bld = ((-sin(self.PD / 180.)) + sin(declin) * sin(self.Lat)) / (cos(declin) * cos(self.Lat))
83-
84-
tmp2 = ifthenelse(bld < 0, scalar(asin(bld))-360., scalar(asin(bld)))
85-
87+
asin_bld = asin(bld)
88+
tmp2 = ifthenelse(bld < 0, asin_bld-360., asin_bld)
8689
abs_bld = abs(bld)
8790
# daylength [hour]
8891
# abs(bld) > 1. corrects the day length at higher altitudes to 24h
8992
day_length = ifthenelse(abs_bld > 1., scalar(24.), 12. + (24. / 180.) * tmp2)
9093
# Daylength equation can produce MV at high latitudes, this statements sets day length to 0 in that case
91-
day_length = cover(day_length, 0.0)
94+
day_length = cover(day_length, mv=0.0)
9295
# integral of solar height [s] over the day
9396
# abs(bld) > 1. allows correcting the integral of solar height at higher altitudes (north pole)
9497
int_solar_height = ifthenelse(abs_bld > 1., self.int_solar_height_north_pole(day_length, declin), self.int_solar_height_main(day_length, declin))
9598
# Integral of solar height cannot be negative, so truncate at 0
9699
int_solar_height = maximum(int_solar_height, scalar(0.0))
97-
int_solar_height = cover(int_solar_height, 0.0)
100+
int_solar_height = cover(int_solar_height, mv=0.0)
98101
# daily extra-terrestrial radiation (Angot radiation) [J/m2/d]
99-
RadiationAngot = int_solar_height * solar_constant
100-
return RadiationAngot
102+
radiation_angot = int_solar_height * solar_constant
103+
return radiation_angot
101104

102105
def net_absorbed_radiation(self, solar_radiation, radiation_angot):
103106
"""
@@ -115,7 +118,7 @@ def net_absorbed_radiation(self, solar_radiation, radiation_angot):
115118
EmNet = (0.56 - 0.079 * sqrt(self.EAct))
116119
Rso = radiation_angot * (0.75 + (2 * 10 ** -5 * self.Dem))
117120
TransAtm_Allen = solar_radiation / Rso
118-
TransAtm_Allen = cover(TransAtm_Allen, 0)
121+
TransAtm_Allen = cover(TransAtm_Allen, mv=0)
119122
AdjCC = 1.8 * TransAtm_Allen - 0.35
120123
AdjCC = ifthenelse(AdjCC < 0, 0.05, AdjCC)
121124
AdjCC = ifthenelse(AdjCC > 1, 1, AdjCC)
@@ -126,6 +129,26 @@ def net_absorbed_radiation(self, solar_radiation, radiation_angot):
126129

127130
# =========== DYNAMIC ====================================================
128131

132+
def convert_calendar_date_to_gregorian(self, init_t_cal):
133+
"""
134+
Convert calendar_day_start to the default gregorian calendar.
135+
The calculation of the day of the year need to be done using
136+
the gregorian calendar for the angot equation.
137+
"""
138+
calendar_date_gregorian = self.calendar_date
139+
if (isinstance(self.calendar_date, cftime.datetime) and
140+
not init_t_cal in self.gregorian_calendars_id):
141+
if self.calendar_day_start_gregorian is None:
142+
self.calendar_day_start_gregorian = datetime.datetime(year=self.calendar_date.year,
143+
month=self.calendar_date.month,
144+
day=self.calendar_date.day,
145+
hour=self.calendar_date.hour,
146+
minute=self.calendar_date.minute,
147+
second=self.calendar_date.second)
148+
calendar_date_gregorian = (self.calendar_day_start_gregorian +
149+
datetime.timedelta(seconds=(self.currentTimeStep()-1) * self.DtSec))
150+
return calendar_date_gregorian
151+
129152
def dynamic(self):
130153
""" Dynamic part of LISVAP
131154
calls the dynamic part of the hydrological modules
@@ -136,12 +159,14 @@ def dynamic(self):
136159
tp.timemeasure('Start dynamic')
137160
# CM: date corresponding to the model time step (yyyy-mm-dd hh:mm:ss)
138161
self.calendar_date = self.calendar_day_start + datetime.timedelta(seconds=(self.currentTimeStep()) * self.DtSec)
162+
163+
# To get the day of the year we need to convert to gregorian whatever calendar type (360_day, no_leap) so the
164+
# radiation is calculated for the full range 1..365
165+
calendar_date_gregorian = self.convert_calendar_date_to_gregorian(settings.binding['internal.time.calendar'])
139166
# CM: day of the year corresponding to the model time step
140-
self.calendar_day = int(self.calendar_date.strftime("%j"))
167+
self.calendar_day = int(calendar_date_gregorian.strftime("%j"))
141168

142-
# correct method to calculate the day of the year
143169
# CM: model time step
144-
145170
i = self.currentTimeStep()
146171
self.time_since_start = self.currentTimeStep() - self.firstTimeStep() + 1
147172

0 commit comments

Comments
 (0)