Skip to content
Snippets Groups Projects
Commit a63f7cc6 authored by Aiko Voigt's avatar Aiko Voigt
Browse files

Commits git repo for long-term archiving

parents
Branches
No related tags found
No related merge requests found
Showing
with 428 additions and 0 deletions
The repo is licensed under the Creative Commons 4.0 BY licence, https://creativecommons.org/licenses/by/4.0/legalcode.
# Code repository for publication "Tug-of-war on idealized midlatitude cyclones between radiative heating from low-level and high-level clouds
The paper is published in the AGU journal Geophysical Research Letters.
The paper is based on the MSc thesis of Klara Butz, 2022, The radiative impact of clouds on idealized extratropical cyclones, https://ubdata.univie.ac.at/AC16593523.
**Author:** Aiko Voigt, IMG, University of Vienna, aiko.voigt@univie.ac.at.
**The repository contains:**
* sims: ICON simulation logfiles
* literature-search: excel sheet for Web of Science search of previous studies of baroclinic life cyles
* python3: python3 scripts and notebooks, with the following subdirectories:
- plots4paper: jupyter notebooks for figures used in GRL paper, also figure pdfs postprocessed with inkscape and pdfcrop; mostly one jupyter notebook per figure, but some of the supplmementary figures are done in a joint notebook
- environment: yml file that documents the python environment used at DKRZ jupyterhub server
- helperfuncs: python functions to faciliate loading of data and analysis, e.g., domain means, cyclone metrics, averaging of the 6 cyclones into a mean cyclone etc.
To load the simulations, use the dictionary dict_expid defined in myfunctions.py, which links the simulation id to the simulation setup. For an illustration of how to use the dictionary, see for example figure_5.ipynb.
**Simulation summary:** (also found in python3/helperfuncs/myfunctions.py)
*ICON2.1:*
* i2100-0001: ICON2.1 ALLSKYRAD (all-sky radiative heating)
* i2100-0002: ICON2.1 NORAD (no radiative heating)
* i2100-0003: ICON2.1 CLRSKYRAD (clear-sky radiative heating)
* i2100-0004: ICON2.1 CRHONLY (cloud-radiative heating following the Keshtgar et al., 2023 method)
* i2100-0005: ICON2.1 PBLCRHONLY (cloud-radiative heating below 2 km)
* i2100-0007: ICON2.1 CRHONLY (same as i2100-0004 except for the output)
* i2100-0008: ICON2.1 FTCRHONLY (cloud-radiative heating above 2km)
*ICON2.6:*
* i2622-v2-0001: ICON2.6 ALLSKYRAD
* i2622-v2-0003: ICON2.6 NORAD
* i2622-v2-0004: ICON2.6 CLRSKYRAD
* i2622-v2-0005: ICON2.6 CRHONLY
* i2622-v2-0006: ICON2.6 PBLCRHONLY
* i2622-v2-0007: ICON2.6 FTCRHONLY
File added
This diff is collapsed.
%% Cell type:markdown id:6d1830bb-ca0c-48cc-a6f1-5c603e074bf8 tags:
# Save environment to yml file
The kernel used for the analysis is
![grafik.png](attachment:21ab7ce2-87e7-4fed-af18-f4f6224e2125.png)
%% Cell type:markdown id:38517e21-eb5c-4f6a-bfe3-3d61b2779081 tags:
Check python version.
%% Cell type:code id:4962fae4-cb52-4074-9cd1-a535876f6670 tags:
``` python
!which python3
```
%% Output
/sw/spack-levante/mambaforge-4.11.0-0-Linux-x86_64-sobz6z/bin/python3
%% Cell type:markdown id:e2bd0b66-769f-48cb-bcf4-2f869f6f15cb tags:
Available kernel(s).
%% Cell type:code id:a1a4afa6-970f-44f6-a9ad-bee47d682227 tags:
``` python
!jupyter kernelspec list
```
%% Output
Available kernels:
python3 /sw/spack-levante/mambaforge-4.11.0-0-Linux-x86_64-sobz6z/share/jupyter/kernels/python3
%% Cell type:markdown id:2e00bf50-b4ec-4723-acef-906ff91f7db9 tags:
Print available conda environments, in our case the only environment is "base".
%% Cell type:code id:b3e48b68-2b78-401b-840f-df447ff11ed3 tags:
``` python
!conda env list
```
%% Output
# conda environments:
#
base * /sw/spack-levante/mambaforge-4.11.0-0-Linux-x86_64-sobz6z
%% Cell type:markdown id:7e45809d-bd8f-4da0-9343-f94af4041b3f tags:
Save environment to yaml file and check it is there.
%% Cell type:code id:357a97ec-215f-4ef9-bbd6-f60d7c42f18e tags:
``` python
!conda env export --name base > ./environment.yml
!ls
```
%% Output
environment2yml.ipynb environment.yml
# a set of functions to analyse cyclone characteristics and for plotting
# author: Aiko Voigt, UNIVIE
import xarray as xr
import numpy as np
import matplotlib.pyplot as plt
# dictionary of expids and their meaning
dict_expid = dict([
("i2100-0001", "ICON2.1 ALLSKYRAD"),
("i2100-0002", "ICON2.1 NORAD"),
("i2100-0003", "ICON2.1 CLRSKYRAD"),
("i2100-0004", "ICON2.1 CRHONLY"),
("i2100-0005", "ICON2.1 PBLCRHONLY"),
("i2100-0007", "ICON2.1 CRHONLY"), # same as 0004 except for the output
("i2100-0008", "ICON2.1 FTCRHONLY"),
("i2622-v2-0001", "ICON2.6 ALLSKYRAD"),
("i2622-v2-0003", "ICON2.6 NORAD"),
("i2622-v2-0004", "ICON2.6 CLRSKYRAD"),
("i2622-v2-0005", "ICON2.6 CRHONLY"),
("i2622-v2-0006", "ICON2.6 PBLCRHONLY"),
("i2622-v2-0007", "ICON2.6 FTCRHONLY")
])
def domainmean(data):
"""Computes domain mean for a single level, not taking into account the density of air at that level.
Therefore, the function is appropriate for fields at the Earth's surface or at the top of the atmosphere.
The computation is done over entire domain. For a global dataset, this means that
a global mean will be computed. For a regional dataset, the regional domain will
be computed."
"""
if "lon" in data.dims: data = data.mean("lon")
weights = np.cos(np.deg2rad(data.lat))
weights.name = "weights"
data_weighted = data.weighted(weights)
return data_weighted.mean("lat")
def domainmean_per_level(data, rho):
"""Computes domain mean for a single level, taking into account the density of air at that level.
Therefore, the function is appropriate for fields at model levels, e.g., atmospheric temperature.
The function is the counterpart to the function "domainmean".
The computation is done over entire domain. For a global dataset, this means that
a global mean will be computed. For a regional dataset, the regional domain will
be computed."
"""
data.load(); rho.load()
weights = np.cos(np.deg2rad(data.lat)) * rho
weights.name = "weights"
data_weighted = data.weighted(weights)
if "lon" in data.dims:
mean = data_weighted.mean(["lat", "lon"])
else:
mean = data_weighted.mean("lat")
return mean
def make_meancyclone(da):
"""Constructs a mean cyclone from the 6 simulated cyclones for data array.
This is done by averaging over the 6 longitudes that are separated by dlon=60 deg."""
# make lon the last dimension
_da = da.transpose(...,"lon")
_shape = _da.shape # this is (something, lon)
_nlon = _shape[-1]
_mean = _da.values.reshape(_shape[:-1]+(6,int(_nlon/6))).mean(axis=len(_shape)-1)
# return as data array
_lon0 = _da["lon"][0] # to make lon start at 0
return xr.DataArray( _mean, dims = _da.dims,
coords= {**_da.drop("lon").coords, **{"lon":-_lon0+_da.coords["lon"][0:int(_nlon/6)]}} )
def eke(data, meancyclone=False):
"""Computes eddy kinetic energy,
Eddies are defined as the deviation from the zonal mean at a given timestep."""
u = data["u"]
v = data["v"]
if meancyclone:
u = make_meancyclone(u)
v = make_meancyclone(v)
aux = 0.5*( np.power(u-u.mean("lon"),2) + np.power(v-v.mean("lon"),2))
return aux.mean("lon")
def corepressure(data, meancyclone=False):
"""Computes cyclone core pressure as the global minimum of variable pres_sfc.
Default is to consider the surface pressure of the mean cyclone, if this is not desired the
optional argument meancyclone needs to be set to False."""
pres_sfc = data["pres_sfc"]
if meancyclone: pres_sfc = make_meancyclone(data["pres_sfc"])
return pres_sfc.min(("lat","lon"))
def zlevels():
"""Return height of z-levels in km above ground
As z-levels are the same for all simulations, we can just get them from one simulation."""
path = "/work/bb1135/b380459/Butz-MSc2022/consolidated"
z_ifc = xr.open_dataset(path+"/icon-lc1-i2100-0002-atmfields_ll_DOM01_ML.nc")["z_ifc"].squeeze().isel(time=0, lat=0, lon=0).values
z = 0.5*(z_ifc[1:91]+z_ifc[0:90]) * 1e-3
return z
def calculate_crh(ds, expid):
"""Calculates cloud-radiative heating within the atmosphere.
The calculation of CRH depends on the simulation, hence the expid is necessary. crhlw is the longwave part of the heating."""
if expid=="i2100-0001":
ds["crhlw"] = 86400*(ds["ddt_temp_radlw"]-ds["ddt_temp_radlwcs"])
ds["crh"] = 86400*(ds["ddt_temp_radlw"]-ds["ddt_temp_radlwcs"]+ds["ddt_temp_radsw"]-ds["ddt_temp_radswcs"])
elif expid=="i2100-0002":
ds["crhlw"] = 86400*(ds["ddt_temp_radlwnw"]-ds["ddt_temp_radlwcs"])
ds["crh"] = 86400*(ds["ddt_temp_radlwnw"]-ds["ddt_temp_radlwcs"]+ds["ddt_temp_radswnw"]-ds["ddt_temp_radswcs"])
elif expid=="i2100-0003":
ds["crhlw"] = 0*ds["ddt_temp_radlw"]
ds["crh"] = 0*ds["ddt_temp_radlw"]
elif expid=="i2100-0004":
ds["crhlw"] = 86400*(ds["ddt_temp_radlwnw"]-ds["ddt_temp_radlwcs"])
ds["crh"] = 86400*(ds["ddt_temp_radlwnw"]-ds["ddt_temp_radlwcs"]+ds["ddt_temp_radswnw"]-ds["ddt_temp_radswcs"])
elif expid=="i2100-0005":
ds["crhlw"] = 86400*ds["ddt_temp_radlw"]
ds["crh"] = 86400*(ds["ddt_temp_radlw"]+ds["ddt_temp_radsw"])
elif expid=="i2100-0007":
ds["crhlw"] = 86400*ds["ddt_temp_radlw"]
ds["crh"] = 86400*(ds["ddt_temp_radlw"]+ds["ddt_temp_radsw"])
elif expid=="i2100-0008":
ds["crhlw"] = 86400*ds["ddt_temp_radlw"]
ds["crh"] = 86400*(ds["ddt_temp_radlw"]+ds["ddt_temp_radsw"])
elif expid=="i2622-v2-0001":
ds["crhlw"] = 86400*(ds["ddt_temp_radlw"]-ds["ddt_temp_radlwcs"])
ds["crh"] = 86400*(ds["ddt_temp_radlw"]-ds["ddt_temp_radlwcs"]+ds["ddt_temp_radsw"]-ds["ddt_temp_radswcs"])
elif expid=="i2622-v2-0003":
ds["crhlw"] = 86400*(ds["ddt_temp_radlwnw"]-ds["ddt_temp_radlwcs"])
ds["crh"] = 86400*(ds["ddt_temp_radlwnw"]-ds["ddt_temp_radlwcs"]+ds["ddt_temp_radswnw"]-ds["ddt_temp_radswcs"])
elif expid=="i2622-v2-0004":
ds["crhlw"] = 0*ds["ddt_temp_radlw"]
ds["crh"] = 0*ds["ddt_temp_radlw"]
elif expid=="i2622-v2-0005":
ds["crhlw"] = 86400*ds["ddt_temp_radlw"]
ds["crh"] = 86400*(ds["ddt_temp_radlw"]+ds["ddt_temp_radsw"])
elif expid=="i2622-v2-0006":
ds["crhlw"] = 86400*ds["ddt_temp_radlw"]
ds["crh"] = 86400*(ds["ddt_temp_radlw"]+ds["ddt_temp_radsw"])
elif expid=="i2622-v2-0007":
ds["crhlw"] = 86400*ds["ddt_temp_radlw"]
ds["crh"] = 86400*(ds["ddt_temp_radlw"]+ds["ddt_temp_radsw"])
else:
print("ERROR: unknown expid in recalculate_crh")
return ds
def calculate_latentheating(ds):
"""Calculates latent heating in K/day"""
ds["latheat"] = 86400*ds["ddt_temp_mphy"]
return ds
def calculate_moistheating(ds):
"""Calculates heating related to moist processes of microphysics, saturation adjustment and moist convection in K/day"""
ds["moistheat"] = 86400*(ds["ddt_temp_mphy"]+ds["ddt_temp_pconv"])
return ds
def load_data_cyclonemetrics(expid: str, chunks=None):
"""Loads data for simulation "expid" to compute cyclone metrics of EKE at 300hPa and 925hPa and core pressure.
Returns an xarray dataset.
"chunks" can be set to use lazy loading with dask. If None, dask is nit used.
"expid" is added as an attribute to the dataset.
"setup" is an added as an attribute based on expid and dict_expid."""
path = "/work/bb1135/b380459/Butz-MSc2022/consolidated"
sfc = xr.open_dataset(path+"/icon-lc1-"+expid+"-sfcfields_ll_DOM01_ML.nc", chunks=chunks).squeeze()
atm = xr.open_dataset(path+"/icon-lc1-"+expid+"-atmfields_ll_DOM01_PL.nc", chunks=chunks).rename({"plev_3": "lev"}).squeeze()
ds = xr.merge([sfc,atm]).squeeze()
del sfc, atm
# reorder latitude if not from NP to SP
if ds.lat[0]<ds.lat[1]:
ds = ds.reindex(lat=list(reversed(ds.lat)))
# add expid and setup as global attributes
ds.attrs["expid"]=expid
ds.attrs["setup"]=dict_expid[expid]
return ds
def load_data(expid: str, chunks=None):
"""Loads icon model output for the simulation "expid" as an xarray dataset.
"chunks" can be set to use lazy loading with dask. If None, dask is not used.
"expid" is added as an attribute to the dataset.
"setup" is an added as an attribute based on expid and dict_expid."""
path = "/work/bb1135/b380459/Butz-MSc2022/consolidated"
sfc = xr.open_dataset(path+"/icon-lc1-"+expid+"-sfcfields_ll_DOM01_ML.nc", chunks=chunks).squeeze()
atm = xr.open_dataset(path+"/icon-lc1-"+expid+"-atmfields_ll_DOM01_ML.nc", chunks=chunks).squeeze()
dTdt = xr.open_dataset(path+"/icon-lc1-"+expid+"-phys_tendencies_ll_DOM01_ML.nc", chunks=chunks).squeeze()
ds = xr.merge([sfc,atm,dTdt]).squeeze()
del sfc, atm, dTdt
# reorder latitude if not from NP to SP
if ds.lat[0]<ds.lat[1]:
ds = ds.reindex(lat=list(reversed(ds.lat)))
# add cloud-radiative heating, latent heating and heating from moist processes
ds = calculate_crh(ds, expid=expid)
ds = calculate_latentheating(ds)
ds = calculate_moistheating(ds)
# vertical velocity at 1, 2, 5, 8 and 10 km height
z = zlevels()
def add_omega_at_km(height):
iheight = np.argmin(np.abs(z-height))
ds["w"+str(height)+"km"] = ds["omega"].isel(height=iheight)
ds["w"+str(height)+"km"] = ds["omega"].isel(height=iheight)
for height in [1,2,5,8,10]:
add_omega_at_km(height)
# add expid and setup as global attributes
ds.attrs["expid"]=expid
ds.attrs["setup"]=dict_expid[expid]
return ds
def beautify_timeseries(ax):
""" Makes plots of time series nicer in terms of axes and labeling"""
# adjust spines
ax = plt.gca()
ax.spines['top'].set_color('none')
ax.spines['right'].set_color('none')
ax.xaxis.set_ticks_position('bottom')
ax.spines['bottom'].set_position(('data',0))
plt.xlim(0,11)
plt.xlabel("day", loc="right", size=12)
plt.xticks([0,1,2,3,4,5,6,7,8,9,10,11],["", "1","","3","","5","","7","","9","","11"], size=10);
\ No newline at end of file
File added
This diff is collapsed.
File added
File added
This diff is collapsed.
File added
File added
File added
File added
This diff is collapsed.
File added
File added
This diff is collapsed.
File added
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment