Skip to content

GitLab

Commit 94b6072d authored by Paul Gierz's avatar Paul Gierz
Browse files
Options

merge fixup

parents f6309102 1da98a63
Expand all Hide whitespace changes
Inline Side-by-side
Showing with 540 additions and 144 deletions
functions/iterative_coupling/coupling_fesom2ice.functions
View file @ 94b6072d
......@@ -301,15 +301,6 @@ function fesom_ice_fesom_build_grid_with_lonlat_depth {
fi
done
# Several things might be needed for pyfesom. If it isn't installed;
# get it
which pip && \
for required_python_external in joblib cmocean; do
python -c "import $required_python_external" \
|| pip install --user ${required_python_external} \
&& echo -e "\t\t - Using existing ${required_python_external} for python"
done
processed_fesom1x_file3D_T_list_for_ice=""
processed_fesom1x_file3D_S_list_for_ice=""
processed_fesom1x_file2D_list_for_ice=""
......@@ -372,6 +363,15 @@ function fesom_ice_fesom_build_grid_with_lonlat_depth {
echo " ## Do not activate HOT fix"
echo "=====## ---------------------------------------------------------------"
fi
# Several things might be needed for pyfesom. If it isn't installed;
# get it
which pip && \
for required_python_external in joblib cmocean seawater; do
python -c "import $required_python_external" \
|| pip install --user ${required_python_external} \
&& echo -e "\t\t - Using existing ${required_python_external} for python"
done
# --------------------------------------------------------
# ------------------------------------------------
......
functions/iterative_coupling/coupling_ice2echam.functions
View file @ 94b6072d
......@@ -344,7 +344,7 @@ set_glacial_mask_jsbach_update_init_file() {
res_oce=GENERIC
ntiles=11 # number of jsbach tiles
dynveg=true # setup for dynamic vegetation
dynveg=${ldynveg_esm} # setup for dynamic vegetation
c3c4crop=true # differentiate between C3 and C4 crops
lpasture=true # distinguish pastures from grasses
......@@ -481,7 +481,7 @@ EOF
echo "error in jsbach_init_file"
exit 1
fi
mv jsbach_${res_atm}_11tiles_5layers_1850.nc ${COUPLE_DIR}/jsbach_init_file_${CHUNK_START_DATE_echam}-${CHUNK_END_DATE_echam}.nc
mv jsbach_${res_atm}${res_oce}_11tiles_5layers_1850${jsbach_nc}.nc ${COUPLE_DIR}/jsbach_init_file_${CHUNK_START_DATE_echam}-${CHUNK_END_DATE_echam}.nc
echo :> ${COUPLE_DIR}/jsbach_init_override_${CHUNK_START_DATE_echam}-${CHUNK_END_DATE_echam}.dat
echo "LAND_BOUNDARY_CONDITIONS_jsbach=${COUPLE_DIR}/jsbach_init_file_${CHUNK_START_DATE_echam}-${CHUNK_END_DATE_echam}.nc" > ${COUPLE_DIR}/jsbach_init_override_${CHUNK_START_DATE_echam}-${CHUNK_END_DATE_echam}.dat
add_to ${COUPLE_DIR}/jsbach_init_override_${CHUNK_START_DATE_echam}-${CHUNK_END_DATE_echam}.dat jsbach_init_override.dat config jsbach
......@@ -685,7 +685,18 @@ set_in_jsbach_restart_file() {
echo " in: $SOURCE_FILENAME"
if [ ! -z $unpack ]; then
echo " - with unpacking from landpoint grid to latlon grid"
python ${FUNCTION_PATH}/unpack_jsbach.py $TARGET_VARNAME $RESTART_FILE
###################################################
# LA: There is a syntax error with the asteriks in the 'pack_jsbach.py' for newer
# python 2.7 versions
##################################################
module unload python3
module load python3
echo "########################"
echo "python version:"
pyv="$(python3 -V 2>&1)"
echo "$pyv"
python3 ${FUNCTION_PATH}/unpack_jsbach.py $TARGET_VARNAME $RESTART_FILE
################################################################################
# $ cdo -h replace
# NAME
......@@ -705,7 +716,7 @@ set_in_jsbach_restart_file() {
${TARGET_VARNAME}_lonlat_grid_replaced.nc
echo " - and repacking from latlon grid to landpoint grid"
# FIXME: The next line is probably wrong, or incomplete.
python ${FUNCTION_PATH}/pack_jsbach.py ${TARGET_VARNAME} ${TARGET_VARNAME}_lonlat_grid_replaced.nc ${RESTART_FILE}
python3 ${FUNCTION_PATH}/pack_jsbach.py ${TARGET_VARNAME} ${TARGET_VARNAME}_lonlat_grid_replaced.nc ${RESTART_FILE}
else
cdo -s replace \
${RESTART_FILE} \
......
functions/iterative_coupling/fesom_scalar_array_to_LonLat.py
View file @ 94b6072d
......@@ -134,30 +134,41 @@ time_out = FID.variables['time'][0:no_timesteps]
#
# Full hydrographic fields
#
for itime in np.arange(0, no_timesteps, 1, dtype=np.int32):
##time_out[itime] = FID.variables['time'][itime]
time_read = time_out[itime]
# Some information
print('* TIME('+str(itime)+') = '+str(time_read))
ilevel = -1
for depth in mesh.zlevs:
#
# Check shape of input file to process old (2D array) and new (3D array) FESOM files
#
print('FID.variables[args.FESOM_VARIABLE[0]][itime, ilevel].shape',
FID.variables[args.FESOM_VARIABLE[0]].shape)
if len(FID.variables[args.FESOM_VARIABLE[0]].shape) == 2:
for itime in np.arange(0, no_timesteps, 1, dtype=np.int32):
##time_out[itime] = FID.variables['time'][itime]
time_read = time_out[itime]
# Some information
idepth = np.int(depth)
print('* depth='+str(idepth)+' ('+str(depth)+\
') ++ time('+str(itime)+') = '+str(time_read))
#
# Prepare data for final netcdf output
#
ilevel = ilevel + 1
flag_verbose=False
level_data, elem_no_nan = \
pf.get_data(FID.variables[args.FESOM_VARIABLE[0]][itime, :], mesh, idepth, flag_verbose)
level_data[np.where(np.isnan(level_data))] = NAN_REPLACE
TempFields_out[itime, ilevel, :] = level_data
print('* TIME('+str(itime)+') = '+str(time_read))
ilevel = -1
for depth in mesh.zlevs:
# Some information
idepth = np.int(depth)
print('* depth='+str(idepth)+' ('+str(depth)+\
') ++ time('+str(itime)+') = '+str(time_read))
#
# Prepare data for final netcdf output
#
ilevel = ilevel + 1
flag_verbose=False
level_data, elem_no_nan = \
pf.get_data(FID.variables[args.FESOM_VARIABLE[0]][itime, :], mesh, idepth, flag_verbose)
level_data[np.where(np.isnan(level_data))] = NAN_REPLACE
TempFields_out[itime, ilevel, :] = level_data
elif len(FID.variables[args.FESOM_VARIABLE[0]].shape) == 3:
TempFields_out = FID.variables[args.FESOM_VARIABLE[0]][:]
else:
print('Wrong shape of input file: ', FID)
# ----------------------------------------------------------------
#
......
functions/jsbach/jsbach.functions
View file @ 94b6072d
......@@ -17,6 +17,7 @@ jsbach_set_defaults()
if [[ "x${DYNVEG_jsbach}" = "xdynveg" ]]; then
jsbach_nc="_dynveg"
ldynveg_esm=true
jsbach_message="\t\tRunning WITH dynamic vegetation"
# CD: workaround for wrong filenames in pool
jsbach_nc=""
......@@ -29,7 +30,12 @@ jsbach_set_defaults()
jsbach_nc=""
fi
fi
jsbach_message="\t\tRunning WITHOUT dynamic vegetation"
jsbach_message="\t\tRunning WITHOUT dynamic vegetation"
ldynveg_esm=false
jsbach_nc="_no-dynveg"
# PG: The next line below is conflicting, but I don't know how
# to set it correctly right now....
jsbach_message="\t\tRunning WITHOUT dynamic vegetation"
fi
case $jsbach_DATASET in
......
functions/pism/index.dat 0 → 100644
View file @ 94b6072d
-120 1
-119 0.99
-118 0.98
-117 0.97
-116 0.96
-115 0.95
-114 0.94
-113 0.93
-112 0.92
-111 0.91
-110 0.9
-109 0.89
-108 0.88
-107 0.87
-106 0.86
-105 0.85
-104 0.84
-103 0.83
-102 0.82
-101 0.81
-100 0.8
-99 0.79
-98 0.78
-97 0.77
-96 0.76
-95 0.75
-94 0.74
-93 0.73
-92 0.72
-91 0.71
-90 0.7
-89 0.69
-88 0.68
-87 0.67
-86 0.66
-85 0.65
-84 0.64
-83 0.63
-82 0.62
-81 0.61
-80 0.6
-79 0.59
-78 0.58
-77 0.57
-76 0.56
-75 0.55
-74 0.54
-73 0.53
-72 0.52
-71 0.51
-70 0.5
-69 0.49
-68 0.48
-67 0.47
-66 0.46
-65 0.45
-64 0.44
-63 0.43
-62 0.42
-61 0.41
-60 0.4
-59 0.39
-58 0.38
-57 0.37
-56 0.36
-55 0.35
-54 0.34
-53 0.33
-52 0.32
-51 0.31
-50 0.3
-49 0.29
-48 0.28
-47 0.27
-46 0.26
-45 0.25
-44 0.24
-43 0.23
-42 0.22
-41 0.21
-40 0.2
-39 0.19
-38 0.18
-37 0.17
-36 0.16
-35 0.15
-34 0.14
-33 0.13
-32 0.12
-31 0.11
-30 0.1
-29 0.09
-28 0.08
-27 0.07
-26 0.06
-25 0.05
-24 0.04
-23 0.03
-22 0.02
-21 0.01
-20 0
functions/pism/pism.functions
View file @ 94b6072d
......@@ -216,6 +216,22 @@ pism_resolution() {
RESOLUTION_OPT_pism="-Mx ${XRES_pism:?Miss variable, Please provide XRES_pism} -My ${YRES_pism:?Miss variable, Please provide YRES_pism} -Mz ${ZRES_pism:?Miss variable, Please provide ZRES_pism} -Mbz ${BZ_pism:?Miss variable, Please provide BZ_pism} -Lz $LZ_pism -Lbz $LBZ_pism $PERIODICITY_OPT_pism"
}
pism_protocol() {
case ${PROTOCOL_pism} in
"ISMIP6")
TS_VARS_pism=${TS_VARS_pism:-"ice_mass,limnsw,tendency_of_ice_mass_due_to_basal_mass_flux,tendency_of_ice_mass_due_to_discharge,tendency_of_ice_mass_due_to_surface_mass_flux"}
EX_INTERVAL_pism=${EX_INTERVAL_pism:-monthly}
EX_VARS_pism=${EX_VARS_pism:-"climatic_mass_balance,ice_surface_temp,sftflf,sftgrf,shelfbtemp,surface_accumulation_flux,surface_melt_fux,surface_runoff_flux,tempbase,thk,topg,usurf,velbar,velbase,velsurf,wvelbase,wvelsurf"}
;;
*)
echo "#=============================================="
echo "#"
echo "# UNKNOWN protocol : >>${PROTOCOL_pism}<<"
echo "#"
echo "#=============================================="
esac
}
pism_prepare_exe() {
add_to ${BIN_DIR_pism}/${EXE_pism} ${EXE_pism}
echo -e "\t\tTaking pism executable BIN_DIR_pism/EXE_pism=${BIN_DIR_pism}/${EXE_pism}"
......@@ -294,9 +310,9 @@ pism_prepare_config() {
if [[ "x$iterative_coupling" == "x1" ]] || [[ "x${offline_index}" == "x1" ]]; then
for couple_type in coupler config_value forcing; do
mecho "Looking for ${COUPLE_DIR}/pism_${couple_type}_${CHUNK_START_DATE_pism}-${CHUNK_END_DATE_pism}.dat"
echo "Looking for ${COUPLE_DIR}/pism_${couple_type}_${CHUNK_START_DATE_pism}-${CHUNK_END_DATE_pism}.dat"
if [ -f ${COUPLE_DIR}/pism_${couple_type}_${CHUNK_START_DATE_pism}-${CHUNK_END_DATE_pism}.dat ]; then
mecho "Found ${COUPLE_DIR}/pism_${couple_type}_${CHUNK_START_DATE_pism}-${CHUNK_END_DATE_pism}.dat"
echo "Found ${COUPLE_DIR}/pism_${couple_type}_${CHUNK_START_DATE_pism}-${CHUNK_END_DATE_pism}.dat"
source ${COUPLE_DIR}/pism_${couple_type}_${CHUNK_START_DATE_pism}-${CHUNK_END_DATE_pism}.dat
else
echo "Could not find ${COUPLE_DIR}/pism_${couple_type}_${CHUNK_START_DATE_pism}-${CHUNK_END_DATE_pism}.dat"
......@@ -401,11 +417,12 @@ pism_prepare_forcing() {
COUPLE_DIR=$(pwd)/fake_couple
cd $COUPLE_DIR
index_file=${index_file:-/home/ollie/ukrebska/python_scripts/index.dat}
index_file=${index_file:-${FUNCTION_PATH}/index.dat}
echo "Index File: $index_file"
echo "Index values at $RUN_NUMBER_pism"
index_values=$(sed "${RUN_NUMBER_pism}q;d" $index_file)
echo $index_values | read kyr index
index_script=${index_script:-${HOME}/first.py}
index_script=${index_script:-${FUNCTION_PATH}/pism_index_script.py}
echo "Using PI_file=$PI_file"
echo "Using LGM_file=$LGM_file"
......
functions/pism/pism_index_script.py 0 → 100755
View file @ 94b6072d
#!/usr/bin/env python
# Import from Standard Library
import argparse
import logging
import sys
# Import Third-Party Packages
from climlab.solar.insolation import daily_insolation
from climlab.solar.orbital import OrbitalTable
import climlab
import numpy as np
import xarray as xr
def parse_arguments():
parser = argparse.ArgumentParser(description='Generates a pseudo ECHAM6 file given an index and two climate endmembers')
parser.add_argument("warm_file")
parser.add_argument("cold_file")
parser.add_argument("index", type=float)
parser.add_argument("kyr", type=int)
parser.add_argument("outfile_name")
parser.add_argument("-v", action="store_true", dest="debug")
return parser.parse_args()
def index_forcing(kyr, index, infile1, infile2, outfile):
"""
# synthesizes atm. forcing files for PISM, based on one index, two climatologies
# (a "warm" and a "cold" climate in infile1 and infile2 respectively) and a specific time in kyrs
# (LGM is -21, does not seem to work for future!!!)
#
# input file is expected to have climatologic monthly temp2, precip and geosp on a lonxlatx12 grid
#
# output forcing goes to netcdf outfile and includes
# temperature T = index*warm.temp2+(1-index)*cold.temp2
# precipitation PP analogue to temp2
# surface elevation z analogue to temp2
# and solar radiation is calculated from climlab for each month on the same grid
"""
logging.debug("kyr=%s" % str(kyr))
logging.debug("type(kyr)=%s" % str(type(kyr)))
logging.debug("index=%s" % str(index))
logging.debug("infile1=%s" % str(infile1))
logging.debug("infile2=%s" % str(infile2))
logging.debug("outfile=%s" % str(outfile))
warm = xr.open_dataset(infile1)
lat = warm.lat
lon = warm.lon
nx = lon.shape[0]
ny = lat.shape[0]
dtype = 1
if dtype == 1:
SOLm = np.zeros((365, ny))
daya = np.arange(1.78, 366.8) # (2.3,367.3)
else:
SOLm = np.zeros((360, ny))
daya = np.arange(0, 360)
srad0dI = np.zeros((12, ny, nx))
sradtest = np.zeros((12, ny, nx))
w_temp = warm.temp2.values
w_aprc = warm.aprc.values
w_aprl = warm.aprl.values
w_z = warm.geosp.values
w_srad0d = warm.srad0d.values
w_sradsu = warm.sradsu.values
w_srads = warm.srads.values
w_alb = warm.albedo.values # albedo
w_sradsd = w_srads - w_sradsu # downward solar radiation at the surface
w_tao = w_sradsd / w_srad0d # else surface/TOA downward solar radiation
cold = xr.open_dataset(infile2)
c_temp = cold.temp2.values
c_aprc = cold.aprc.values
c_aprl = cold.aprl.values
c_z = cold.geosp.values
c_srad0d = cold.srad0d.values
c_sradsu = cold.sradsu.values
c_srads = cold.srads.values
c_alb = cold.albedo.values # albedo
c_sradsd = c_srads - c_sradsu # downward solar radiation at the surface
c_tao = c_sradsd / c_srad0d # else surface/TOA downward solar radiation
# subset of orbital parameters for a specified year in kyr (LGM is -21, careful: not yet defined for kyear>0!)
orb = OrbitalTable.interp(kyear=kyr)
orb2 = OrbitalTable.interp(kyear=0)
orb_PI = {"ecc": 0.016724, "long_peri": 282.157, "obliquity": 23.4468}
orb_LGM = {"ecc": 0.018994, "long_peri": 294.42, "obliquity": 22.949}
# insolation values for the specific year for all days of a year and all latitudes of the atmospheric grid
SOL = daily_insolation(lat=lat, day=daya, orb=orb, S0=1365.0, day_type=dtype)
# compute monthly mean insolation from daily values
if dtype == 1:
days_per_month = [31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31]
firstday = 0
else:
days_per_month = 30
firstday = 360 - 80
# firstday=0
for i in range(0, 12):
if dtype == 1:
days = np.arange(firstday, firstday + days_per_month[i]) # dtype =1
firstday = firstday + days_per_month[i]
else:
if (firstday + days_per_month) < 360:
days = np.arange(firstday, firstday + days_per_month)
firstday = firstday + days_per_month
else:
first_day_of_next_month = firstday + days_per_month - 360
days = np.concatenate((range(firstday, 360), range(0, first_day_of_next_month)))
firstday = first_day_of_next_month
SOLm = SOL[days, :].mean(axis=0)
for j in range(0, ny):
srad0dI[i, j, :] = SOLm[j]
# combine "warm" and "cold" fields according to index
T2I = (index * w_temp) + ((1 - index) * c_temp)
PPCI = index * w_aprc + (1 - index) * c_aprc
PPLI = index * w_aprl + (1 - index) * c_aprl
ZI = index * w_z + (1 - index) * c_z
taoI = index * w_tao + (1 - index) * c_tao # scaled transmissivity taoI
albI = index * w_alb + (1 - index) * c_alb # scaled albedo
sradsdI = (
taoI * srad0dI
) # downward surf. radiation from taoI and TOA insolation from climlab
sradsuI = -albI * sradsdI # but we actually only need UPWARD surf. radiation
# and net surface radiation... I hope we can use albedo here...
sradsI = ( 1 - albI) * sradsdI
out = xr.open_dataset(infile1)
out.temp2.values = T2I
out.aprc.values = PPCI
out.aprl.values = PPLI
out.geosp.values = ZI
out.srads.values = sradsI
out.sradsu.values = sradsuI
out.srad0d.values = srad0dI
out.albedo_vis.values = albI
out.to_netcdf(path=outfile, mode="w", compute=True)
if __name__ == "__main__":
args = parse_arguments()
if args.debug:
logging.basicConfig(level=logging.DEBUG)
index_forcing(args.kyr, args.index, args.warm_file, args.cold_file, args.outfile_name)
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment