Skip to content

Commit 2602590

Browse files
larsbuntemeyerlarsbuntemeyer
committed
added soil computation from ERA5
1 parent 77ab10d commit 2602590

File tree

4 files changed

+295
-17
lines changed

4 files changed

+295
-17
lines changed

‎pyremo/preproc/era5-params.yaml‎

Lines changed: 33 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -1,3 +1,34 @@
1+
# mo_gmgrib
2+
3+
# & 'FIB ', 129 , 1.0 , 0.0 , &
4+
# & 'BLA ', 172 , 1.0 , 0.0 , &
5+
# & 'PHI ', 253 , 57.296 , 0.0 , &
6+
# & 'LAM ', 254 , 57.296 , 0.0 , &
7+
# & 'PS ', 152 , 1.0 , 0.0 , &
8+
# & 'T ', 130 , 1.0 , 0.0 , &
9+
# & 'QD ', 133 , 1.0 , 0.0 , &
10+
# & 'QW ', 246 , 1.0 , 0.0 , &
11+
# & 'U ', 131 , 1.0 , 0.0 , &
12+
# & 'V ', 132 , 1.0 , 0.0 /
13+
14+
# DATA(GMNAME(I), GMGBNR(I), GMGBFK(I), GMGBBS(I) , I = 11,20) / &
15+
# & 'WB1 ', 39 , 1.0 , 0.0 , &
16+
# & 'WB2 ', 40 , 1.0 , 0.0 , &
17+
# & 'WB3 ', 41 , 1.0 , 0.0 , &
18+
# & 'SN ', 141 , 1.0 , 0.0 , &
19+
# & 'WL ', 198 , 1.0 , 0.0 , &
20+
# & 'TS ', 235 , 1.0 , 0.0 , &
21+
# & 'TD3 ', 139 , 1.0 , 0.0 , &
22+
# & 'TD4 ', 170 , 1.0 , 0.0 , &
23+
# & 'TD5 ', 183 , 1.0 , 0.0 , &
24+
# & 'TD ', 236 , 1.0 , 0.0 /
25+
26+
# DATA(GMNAME(I), GMGBNR(I), GMGBFK(I), GMGBBS(I) , I = 21,NFGM) / &
27+
# & 'TSN ', 238 , 1.0 , 0.0 , &
28+
# & 'SEAICE ', 31 , 1.0 , 0.0 , &
29+
# & 'SST ', 34 , 1.0 , 0.0 /
30+
31+
132
defaults:
233
era_id: E5
334
frequency: hourly
@@ -57,9 +88,9 @@ parameters:
5788
dataType: an
5889
era_id: E5
5990
frequency: hourly
60-
gc_name: siconc
91+
gc_name: sea_ice_cover
6192
level_type: surface
62-
short_name: ci
93+
short_name: siconc
6394
skt:
6495
code: 235
6596
dataType: an

‎pyremo/preproc/era5.py‎

Lines changed: 13 additions & 13 deletions
Original file line numberDiff line numberDiff line change
@@ -404,19 +404,19 @@ def era5_from_gcloud(time=None, chunks="auto", hfreq=6):
404404
The selected ERA5 dataset.
405405
"""
406406

407-
era5_to_cmip_cf = {
408-
"temperature": "ta",
409-
"surface_pressure": "ps",
410-
"u_component_of_wind": "ua",
411-
"v_component_of_wind": "va",
412-
"specific_humidity": "hus",
413-
"specific_cloud_liquid_water_content": "clw",
414-
"geopotential_at_surface": "orog",
415-
"land_sea_mask": "sftlf",
416-
"sea_ice_cover": "sic",
417-
"snow_depth": "snd",
418-
"sea_surface_temperature": "tos",
419-
}
407+
# era5_to_cmip_cf = {
408+
# "temperature": "ta",
409+
# "surface_pressure": "ps",
410+
# "u_component_of_wind": "ua",
411+
# "v_component_of_wind": "va",
412+
# "specific_humidity": "hus",
413+
# "specific_cloud_liquid_water_content": "clw",
414+
# "geopotential_at_surface": "orog",
415+
# "land_sea_mask": "sftlf",
416+
# "sea_ice_cover": "sic",
417+
# "snow_depth": "snd",
418+
# "sea_surface_temperature": "tos",
419+
# }
420420

421421
era5_to_cmip_cf = {v["gc_name"]: k for k, v in era5_params.items()}
422422

‎pyremo/preproc/physics.py‎

Lines changed: 184 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -7,12 +7,12 @@
77

88
from . import constants as const
99

10+
# Constants
1011
zds3 = (3.25 - 0.0) / (3.5 - 0.0)
1112
zds4 = (19.2 - 17.5) / (64.0 - 17.5)
1213
zds5 = (77.55 - 64.0) / (194.5 - 64.0)
1314
zdtfak = 0.00065
1415

15-
1616
# freezing and melting point for
1717
# seaice derivation
1818
frozen = 271.37
@@ -220,3 +220,186 @@ def derive_soil_temperatures(tds, ads):
220220
# tdhm(ij) = tdeh(ij) + zdts(ij)
221221
# tdclhm(ij) = tdcleh(ij) + zdts(ij)
222222
# ENDDO
223+
224+
225+
# def calculate_dp(ps, ak, bk):
226+
# return ps - (ak[-1] + bk[-1] * ps)
227+
228+
229+
def calculate_zdts(fibem, fibge):
230+
"""height correction"""
231+
return (fibem - fibge) * zdtfak
232+
233+
234+
def initialize_soil_temperatures(
235+
tdge, tswge, tslge, td3ge, td4ge, td5ge, fibem, fibge, iland
236+
):
237+
zdts = calculate_zdts(fibem, fibge)
238+
tslem = xr.where(iland == 0, tswge, tslge - zdts)
239+
tswem = tswge.copy()
240+
tsnem = xr.where(iland == 0, tswge, tslem)
241+
td3em = xr.where(iland == 0, tswge, tslge - ((tslge - td3ge) * zds3) - zdts)
242+
td4em = xr.where(iland == 0, tswge, td4ge - ((td4ge - td5ge) * zds4) - zdts)
243+
td5em = xr.where(iland == 0, tswge, td5ge - ((td5ge - tdge) * zds5) - zdts)
244+
tdem = td5em.copy()
245+
tdclem = td5em.copy()
246+
return tslem, tswem, tsnem, td3em, td4em, td5em, tdem, tdclem
247+
248+
249+
def calculate_seaice(tswem):
250+
seaem = xr.zeros_like(tswem)
251+
seaem = xr.where(tswem > melt, 0.0, seaem)
252+
seaem = xr.where(tswem < frozen, 1.0, seaem)
253+
seaem = xr.where(
254+
(tswem >= frozen) & (tswem <= melt), (melt - tswem) / (melt - frozen), seaem
255+
)
256+
return seaem
257+
258+
259+
def calculate_tsi(tswem):
260+
tsiem = tswem.copy()
261+
tswem = xr.where(tswem < frozen, frozen, tswem)
262+
tsiem = xr.where(tsiem > frozen, frozen, tsiem)
263+
return tswem, tsiem
264+
265+
266+
def diagnose_ice_mask(snem):
267+
return xr.where(snem > 9.50, 1.0, 0.0)
268+
269+
270+
def bodfld(
271+
psge,
272+
psem,
273+
akgm,
274+
bkgm,
275+
akem,
276+
bkem,
277+
snge,
278+
wsge,
279+
wsmxem,
280+
wlge,
281+
blaem,
282+
fibem,
283+
fibge,
284+
tswge,
285+
tslge,
286+
td3ge,
287+
td4ge,
288+
td5ge,
289+
tdge,
290+
prgpk=False,
291+
ipr=None,
292+
jpr=None,
293+
):
294+
# Calculate dpge and dpem
295+
# dpge = calculate_dp(psge, akgm, bkgm)
296+
# dpem = calculate_dp(psem, akem, bkem)
297+
298+
# fak1 = 0.07_DP
299+
# fak2 = 0.21_DP
300+
# fak3 = 0.72_DP
301+
302+
# WSGm(ij) = (fak1*zwb1(ij)) + (fak2*zwb2(ij)) + (fak3*zwb3(ij))
303+
304+
# Initialize arrays
305+
snem = snge.copy()
306+
wsem = xr.ufuncs.minimum(wsmxem, wsge * wsmxem)
307+
wlem = wlge.copy()
308+
iland = xr.ufuncs.round(blaem * 1.0).astype(
309+
int
310+
) # Assuming FAKINF is 1.0 for simplicity
311+
312+
# Initialize temperatures
313+
tslem, tswem, tsnem, td3em, td4em, td5em, tdem, tdclem = (
314+
initialize_soil_temperatures(
315+
tswge, tslge, td3ge, td4ge, td5ge, fibem, fibge, iland
316+
)
317+
)
318+
319+
# Calculate sea ice and tsi
320+
seaem = calculate_seaice(tswem)
321+
tswem, tsiem = calculate_tsi(tswem)
322+
323+
# Diagnose ice mask from snow height
324+
glacem = diagnose_ice_mask(snem)
325+
326+
return (
327+
tslem,
328+
tswem,
329+
tsnem,
330+
td3em,
331+
td4em,
332+
td5em,
333+
tdem,
334+
tdclem,
335+
seaem,
336+
tsiem,
337+
glacem,
338+
wsem,
339+
wlem,
340+
snem,
341+
)
342+
343+
# results = bodfld(psge, psem, akgm, bkgm, akem, bkem, snge, wsge, wsmxem, wlge, blaem, fibem, fibge, tswge, tslge, td3ge, td4ge, td5ge, tdge)
344+
# print(results)
345+
346+
347+
# & 'FIB ', 129 , 1.0 , 0.0 , &
348+
# & 'BLA ', 172 , 1.0 , 0.0 , &
349+
# & 'PHI ', 253 , 57.296 , 0.0 , &
350+
# & 'LAM ', 254 , 57.296 , 0.0 , &
351+
# & 'PS ', 152 , 1.0 , 0.0 , &
352+
# & 'T ', 130 , 1.0 , 0.0 , &
353+
# & 'QD ', 133 , 1.0 , 0.0 , &
354+
# & 'QW ', 246 , 1.0 , 0.0 , &
355+
# & 'U ', 131 , 1.0 , 0.0 , &
356+
# & 'V ', 132 , 1.0 , 0.0 /
357+
358+
# DATA(GMNAME(I), GMGBNR(I), GMGBFK(I), GMGBBS(I) , I = 11,20) / &
359+
# & 'WB1 ', 39 , 1.0 , 0.0 , &
360+
# & 'WB2 ', 40 , 1.0 , 0.0 , &
361+
# & 'WB3 ', 41 , 1.0 , 0.0 , &
362+
# & 'SN ', 141 , 1.0 , 0.0 , &
363+
# & 'WL ', 198 , 1.0 , 0.0 , &
364+
# & 'TS ', 235 , 1.0 , 0.0 , &
365+
# & 'TD3 ', 139 , 1.0 , 0.0 , &
366+
# & 'TD4 ', 170 , 1.0 , 0.0 , &
367+
# & 'TD5 ', 183 , 1.0 , 0.0 , &
368+
# & 'TD ', 236 , 1.0 , 0.0 /
369+
370+
# DATA(GMNAME(I), GMGBNR(I), GMGBFK(I), GMGBBS(I) , I = 21,NFGM) / &
371+
# & 'TSN ', 238 , 1.0 , 0.0 , &
372+
# & 'SEAICE ', 31 , 1.0 , 0.0 , &
373+
# & 'SST ', 34 , 1.0 , 0.0 /
374+
375+
376+
# CALL HIMBLA(TSGm, TSLge, 'TSL ')
377+
# CALL HIMBLA(TSGm, TSWge, 'TSW ')
378+
# CALL HIMBLA(TD3gm, TD3ge, 'TD3 ')
379+
# CALL HIMBLA(TD4gm, TD4ge, 'TD4 ')
380+
# CALL HIMBLA(TD5gm, TD5ge, 'TD5 ')
381+
# CALL HIMBLA(TDGm, TDGe, 'TD ')
382+
# CALL HIMBLA(TSNgm, TSNge, 'TSN ')
383+
# !
384+
# ! SKIN-RESERVOIR OF PLANTS HORIZONTAL INTERPOLIEREN
385+
# CALL HIMBLA(WLGm, WLGe, 'WL ')
386+
# !
387+
# ! ZONALEN WIND HORIZONTAL INTERPOLIEREN UND
388+
# !
389+
# ! SCHNEEHOEHE HORIZONTAL INTERPOLIEREN
390+
# CALL HIMBLA(SNGm, SNGe, 'SN ')
391+
# !
392+
# ! WASSERGEHALT (OBERE SCHICHT) HORIZONTAL INTERPOLIEREN
393+
# CALL HIMBLA(WSGm, WSGe, 'WS ')
394+
# !
395+
# ! SEE-EIS-BEDECKUNG HORIZONTAL INTERPOLIEREN
396+
# CALL HIMBLA(SEAgm, SEAem, 'SEAICE')
397+
# !
398+
# ! TEMPERATUR - DIFFERENZ TP - TB HORIZONTAL INTERPOLIEREN
399+
# CALL HIMBLA(DTPbgm, DTPbge, 'DTPB ')
400+
# !
401+
# ! GEOPOTENTIAL CHECK-FLAECHE HORIZONTAL INTERPOLIEREN
402+
# CALL HIOBLA(FICgm, FICge, 'FIC ', 1)
403+
# !
404+
# ! DRUCKTENDENZ AM BODEN HORIZONTAL INTERPOLIEREN
405+
# CALL HIOBLA(DPDtgm, DPDtge, 'DPDT ', 1)

‎pyremo/preproc/remap.py‎

Lines changed: 65 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -100,7 +100,7 @@ def broadcast_coords(ds, coords=("lon", "lat")):
100100
return lam, phi
101101

102102

103-
def remap(gds, domain_info, vc, surflib):
103+
def remap(gds, domain_info, vc, surflib, initial=False, lice=None):
104104
"""remapping workflow
105105
106106
This function should be similar to the ones in the
@@ -668,6 +668,70 @@ def update_soil_temperatures(ds):
668668
return ds
669669

670670

671+
def create_initial_soil(ds, domain_info, surflib):
672+
"""Create initial soil dataset from atmospheric dataset."""
673+
# rename vertical coordinate of input to avoid conflict with output lev
674+
ds = ds.copy()
675+
676+
fak1 = 0.07
677+
fak2 = 0.21
678+
fak3 = 0.72
679+
# WSGm(ij) = (fak1*zwb1(ij)) + (fak2*zwb2(ij)) + (fak3*zwb3(ij))
680+
wsgm = fak1 * ds.swvl1 + fak2 * ds.swvl2 + fak3 * ds.swvl3
681+
682+
# remove time dimension if there is one
683+
fibem = surflib.FIB.squeeze(drop=True) * const.grav_const
684+
blaem = surflib.BLA.squeeze(drop=True)
685+
686+
lamem, phiem = geo_coords(domain_info, fibem.rlon, fibem.rlat)
687+
688+
# broadcast 1d global coordinates
689+
lamgm, phigm = broadcast_coords(ds)
690+
691+
# compute remap matrix
692+
indii, indjj = intersect(lamgm, phigm, lamem, phiem) # .compute()
693+
694+
# horizontal interpolation
695+
wsge = interpolate_horizontal(
696+
wsgm, lamem, phiem, lamgm, phigm, "WS", indii=indii, indjj=indjj
697+
)
698+
td3ge = interpolate_horizontal(
699+
ds.tsl1, lamem, phiem, lamgm, phigm, "TD3", indii=indii, indjj=indjj
700+
)
701+
td4ge = interpolate_horizontal(
702+
ds.tsl2, lamem, phiem, lamgm, phigm, "TD4", indii=indii, indjj=indjj
703+
)
704+
td5ge = interpolate_horizontal(
705+
ds.tsl3, lamem, phiem, lamgm, phigm, "TD5", indii=indii, indjj=indjj
706+
)
707+
tdge = interpolate_horizontal(
708+
ds.tsl4, lamem, phiem, lamgm, phigm, "TD", indii=indii, indjj=indjj
709+
)
710+
snge = interpolate_horizontal(
711+
ds.snd, lamem, phiem, lamgm, phigm, "SN", indii=indii, indjj=indjj
712+
)
713+
wlge = interpolate_horizontal(
714+
ds.src, lamem, phiem, lamgm, phigm, "WL", indii=indii, indjj=indjj
715+
)
716+
tswge = remap_sst(
717+
ds.tos,
718+
lamem,
719+
phiem,
720+
lamgm,
721+
phigm,
722+
blagm=xr.where(ds.tos.isel(time=0).isnull(), 1.0, 0.0),
723+
blaem=blaem,
724+
)
725+
soil = xr.merge([wsge, td3ge, td4ge, td5ge, tdge, snge, wlge, tswge])
726+
soil.attrs["history"] = "preprocessing with pyremo = {}".format(pr.__version__)
727+
soil.attrs["domain_id"] = domain_info.get("domain_id", "no name")
728+
729+
soil = update_attrs(soil)
730+
731+
# transpose to remo convention
732+
return soil.transpose(..., "rlat", "rlon")
733+
734+
671735
# dpeh(ij) = pseh(ij) - GETP(akem(KEEM),bkem(KEEM),pseh(ij),akem(1))
672736
# dphm(ij) = pshm(ij) - GETP(akhm(KEHM),bkhm(KEHM),pshm(ij),akhm(1))
673737

0 commit comments

Comments
 (0)