Added directory anomalies (readme.md modified accordingly) and file nemo_bdy/prepare_ecearth_hindcasts_paraso.sh

anomalies/calculate_anomalies.sh

@@ -0,0 +1,40 @@
+#!/bin/bash
+set -u
+
+# This script takes ORAS5 and EC-Earth historical ocean temperature and salinity conditions and calculates anomalies from them.
+# Author: D. Verfaillie
+
+# Inputs
+in_folder="/storepelican/dverfail/bsc/anom"
+oras5_file="ORAS5_preinterp_opa1_grid_T_1985-2014"
+histo_1m_thetao="thetao_Omon_EC-Earth3_historical_ensmean_gn_198501-201412"
+histo_1m_so="so_Omon_EC-Earth3_historical_ensmean_gn_198501-201412"
+
+# Prepare ORAS5 file to have the same grid and format as the EC-Earth files
+module --force purge
+module load releases/2020b
+module load noarch
+module load NCO
+module load CDO
+
+cd ${in_folder}/"oras5_clim"
+
+#cdo remapbil,${in_folder}/oras5_clim/ECEarth_grid_T_new.grid ${oras5_file}.nc ${oras5_file}_remapped.nc
+#ncrename -O -d time_counter,time -d x,i -d y,j -d nvertex,vertices -d deptht,lev -v time_counter,time -v nav_lon,longitude -v nav_lat,latitude -v nav_lon_bnds,longitude_bnds -v nav_lat_bnds,latitude_bnds -v deptht,lev -v conservative_temperature,thetao -v absolute_salinity,so ${oras5_file}_remapped.nc
+#ncap2 -O -s 'defdim("bnds",2);make_bounds(lev,$bnds,"lev_bnds");' ${oras5_file}_remapped.nc ${oras5_file}_remapped_bounds.nc
+#ncks -A -v i,j,lev_bnds ../histo_clim/${histo_1m_thetao}.nc ${oras5_file}_remapped_bounds.nc
+ncap2 -s 'lev=double(lev); lev_bnds=double(lev_bnds); longitude=float(longitude); longitude_bnds=float(longitude_bnds); latitude=float(latitude); latitude_bnds=float(latitude_bnds);' ${oras5_file}_remapped_bounds.nc ${oras5_file}_remapped_bounds_reformatted.nc
+ncks -O -x -v so,sossheig ${oras5_file}_remapped_bounds_reformatted.nc thetao_ORAS5_preinterp_opa1_grid_T_1985-2014.nc
+ncks -O -x -v thetao,sossheig ${oras5_file}_remapped_bounds_reformatted.nc so_ORAS5_preinterp_opa1_grid_T_1985-2014.nc
+
+# Calculate monthly and climatological anomalies: histo - oras5 (these anomalies will then have to be SUBSTRACTED from the hindcasts)
+ncdiff ../histo_clim/${histo_1m_so}.nc thetao_ORAS5_preinterp_opa1_grid_T_1985-2014.nc ../anom_histo-oras5/thetao_anomaly_1985-2014_1m.nc
+ncdiff ../histo_clim/${histo_1m_thetao}.nc so_ORAS5_preinterp_opa1_grid_T_1985-2014.nc ../anom_histo-oras5/so_anomaly_1985-2014_1m.nc
+
+cd ${in_folder}/anom_histo-oras5
+cdo ymonmean thetao_anomaly_1985-2014_1m.nc thetao_anomaly_1985-2014_Mm.nc
+cdo ymonmean so_anomaly_1985-2014_1m.nc so_anomaly_1985-2014_Mm.nc
+
+# Reformat the anomaly files to have same format as EC-Earth hindcast ocean temperature and salinity files
+ncrename -O -d time,time_counter -d i,x -d j,y -d vertices,nvertex -d lev,olevel -d bnds,axis_nbounds -v time,time_counter -v time_bnds,time_counter_bounds -v latitude,nav_lat -v latitude_bnds,bounds_nav_lat -v longitude,nav_lon -v longitude_bnds,bounds_nav_lon -v lev,olevel -v lev_bnds,olevel_bounds thetao_anomaly_1985-2014_Mm.nc thetao_anomaly_1985-2014_Mm_reformatted.nc
+ncrename -O -d time,time_counter -d i,x -d j,y -d vertices,nvertex -d lev,olevel -d bnds,axis_nbounds -v time,time_counter -v time_bnds,time_counter_bounds -v latitude,nav_lat -v latitude_bnds,bounds_nav_lat -v longitude,nav_lon -v longitude_bnds,bounds_nav_lon -v lev,olevel -v lev_bnds,olevel_bounds so_anomaly_1985-2014_Mm.nc so_anomaly_1985-2014_Mm_reformatted.nc

anomalies/readme.md

@@ -0,0 +1,13 @@
+# Calculating ocean temperature and salinity anomalies of EC-Earth vs. ORAS5.
+
+Script for calculating ocean temperature and salinity anomalies of EC-Earth historical runs vs. ORAS5.
+
+Anomalies of EC-Earth compared to observations are calculated using NCO and CDO operators in the script `calculate_anomalies.sh`. Please read the NCO and CDO documentation for information on the commands used. The BSC provided us with EC-Earth historical files of ocean temperature and salinity on 75 depth levels, which are ensemble means of their 10 historical members over the period 1985-2014, at monthly resolution. The observations we use consist of the corresponding file (ocean temperature and salinity inside the same netcdf file) from the ORAS5 ocean reanalysis.
+
+The first step in `calculate_anomalies.sh` is to prepare the ORAS5 file to have the same grid (remapping) and format (renaming and reformatting of dimensions and variables, separation into 2 different files for temperature and salinity) as the EC-Earth files. 
+
+Then, the anomalies are calculated at a monthly timestep and on each depth level as the difference between the variable (temperature/salinity) in the EC-Earth historical file and the variable in ORAS5. These anomalies will thus have to be SUBSTRACTED from the hindcasts later on. Note that this step is extremely slow on pelican (it takes several hours).
+
+Climatological anomalies (for mean January, February, etc. over the 1985-2014 period) are then calculated.
+
+A final reformatting step then takes place to have the same format (renaming of dimensions and variables) as the EC-Earth hindcast ocean temperature and salinity files used for the boundary conditions.

nemo_bdy/prepare_ecearth_hindcasts_paraso.sh

@@ -0,0 +1,182 @@
+#!/bin/bash
+set -u
+
+# This script takes EC-Earth outputs, reinterpolates them to the NEMO grid, and then generates BDY dta from it.
+# Author: C. Pelletier
+
+# Inclusive year range
+year_start=1985
+year_stop=2013
+
+# Input prefix (with folder) for input data files (in .tar format).
+in_prefix="/storepelican/dverfail/bsc/anom/hindcasts"
+# Name of experience
+name_exp="a3du"
+fc_num="fc00"
+# Frequency
+freq="1m"
+
+# Working folder where pre-interpolated datafiles will be put
+dest_folder="/storepelican/dverfail/nemo_bdy/bdydta"
+dest_prefix="ECEarth_${name_exp}_${fc_num}_${freq}_eORCA025-SO_BDYDTA_nbr1_"
+
+pwd=`pwd`
+# Python script for reformatting BDYDTA (final step)
+reformat_script="${pwd}/generate_new_bdy_dta_onefile.py"
+
+# min_lon, max_lon, min_lat, max_lat
+# To precut the input data prior to interpolation (not required, but less expensive if your source code is much bigger than the destination one, which is the case from global to the Southern Ocean.
+# Take a generous margin (4-5 nodes bigger than the destination grid limit)
+latlon_cut="-180,180,-90,-28"
+
+# # Script for converting EOS80 to TEOS-10 (look at the script comments for why it exists and what it does).
+# not needed for EC-Earth, which is already in TEOS-10
+# curr_dir=`pwd`
+# eos10_script="${curr_dir}/convert_eos10_new.py"
+
+# NEMO destination CDO grid description (see `cdo_get_griddes.sh`).
+dest_gridfile_t="/storepelican/pelletie/nemo_bdy/grids/eORCA025-SO_grid_T.grid"
+dest_gridfile_u="/storepelican/pelletie/nemo_bdy/grids/eORCA025-SO_grid_U.grid"
+dest_gridfile_v="/storepelican/pelletie/nemo_bdy/grids/eORCA025-SO_grid_V.grid"
+
+pid="$$"
+
+mkdir -p ${dest_folder}
+cd ${dest_folder}
+
+for((year=year_start;year<=year_stop;year++)); do
+
+    # Loop over input years.
+    
+    module purge
+    module load 2018b CDO NCO
+
+    # 1) Extract T, S and ssh (T grid).
+
+    tar -xvf ${in_prefix}/${fc_num}/MMO_${name_exp}_19841101_${fc_num}_${year}0101-${year}1231.tar ${name_exp}_${freq}_${year}0101_${year}1231_opa_grid_T_3D.nc ${name_exp}_${freq}_${year}0101_${year}1231_opa_grid_T_2D.nc
+    
+    cdo selname,thetao,so ${name_exp}_${freq}_${year}0101_${year}1231_opa_grid_T_3D.nc tmp_pid${pid}.nc
+    cdo selname,zos ${name_exp}_${freq}_${year}0101_${year}1231_opa_grid_T_2D.nc tmp2_pid${pid}.nc
+    cdo merge tmp_pid${pid}.nc tmp2_pid${pid}.nc tmp3_pid${pid}.nc
+
+    rm -f ${name_exp}_${freq}_${year}0101_${year}1231_opa_grid_T_3D.nc ${name_exp}_${freq}_${year}0101_${year}1231_opa_grid_T_2D.nc
+    
+    # precut
+    cdo -O sellonlatbox,${latlon_cut} tmp3_pid${pid}.nc tmp_pid${pid}.nc
+
+    ncrename -v olevel,deptht tmp_pid${pid}.nc
+    ncrename -d olevel,deptht tmp_pid${pid}.nc
+    
+    # remap
+    cdo remapbil,${dest_gridfile_t} tmp_pid${pid}.nc tmp2_pid${pid}.nc
+    cdo setmisstonn tmp2_pid${pid}.nc tmp_pid${pid}.nc
+
+    module purge
+    module load ELIC/0.1-foss-2018b-Python-3.6.6
+
+    input_full=`readlink -f tmp_pid${pid}.nc`
+    output_full="${dest_folder}/${dest_prefix}grid_T_y${year}.nc"
+    
+    python ${reformat_script} ${input_full} t ${output_full}
+
+    module purge
+    module load CDO NCO    
+    
+    # 2) Extract u velocities.
+
+    if [ "${freq}" == "1d" ]; then
+	
+	tar -xvf ${in_prefix}/${name_exp}/NEMO/MMO_${name_exp}_18500101_fc0_${year}0101-${year}1231.tar ${name_exp}_${freq}_${year}0101_${year}1231_opa_grid_U_3D.nc ${name_exp}_${freq}_${year}0101_${year}1231_opa_grid_U_2D.nc
+	cdo selname,ubar ${name_exp}_${freq}_${year}0101_${year}1231_opa_grid_U_2D.nc tmp2_pid${pid}.nc
+
+    else
+	# ubar is not stored in monthly outputs, so take it from daily and monmean
+	tar -xvf ${in_prefix}/${name_exp}/NEMO/MMO_${name_exp}_18500101_fc0_${year}0101-${year}1231.tar ${name_exp}_${freq}_${year}0101_${year}1231_opa_grid_U_3D.nc ${name_exp}_1d_${year}0101_${year}1231_opa_grid_U_2D.nc
+	cdo selname,ubar ${name_exp}_1d_${year}0101_${year}1231_opa_grid_U_2D.nc tmp3_pid${pid}.nc
+	cdo monmean tmp3_pid${pid}.nc tmp2_pid${pid}.nc
+	
+    fi    
+    
+    cdo selname,uo ${name_exp}_${freq}_${year}0101_${year}1231_opa_grid_U_3D.nc tmp_pid${pid}.nc
+    cdo merge tmp_pid${pid}.nc tmp2_pid${pid}.nc tmp3_pid${pid}.nc
+
+    rm -f ${name_exp}_${freq}_${year}0101_${year}1231_opa_grid_U_3D.nc ${name_exp}_??_${year}0101_${year}1231_opa_grid_U_2D.nc
+    
+    # precut
+    cdo -O sellonlatbox,${latlon_cut} tmp3_pid${pid}.nc tmp_pid${pid}.nc
+
+    ncap2 -O -s 'uobaroclinic = uo - ubar;' tmp_pid${pid}.nc tmp2_pid${pid}.nc
+    ncks -x -v uo tmp2_pid${pid}.nc -O tmp_pid${pid}.nc
+    ncrename -v ubar,uobarotropic tmp_pid${pid}.nc
+    ncrename -v olevel,depthu -v olevel_bnds,depthu_bnds tmp_pid${pid}.nc
+    ncatted -a bounds,depthu,o,c,"depthu_bnds" -a name,depthu,o,c,"depthu" tmp_pid${pid}.nc
+    ncrename -d olevel,depthu tmp_pid${pid}.nc
+    ncpdq -O -a time_counter,depthu,y,x tmp_pid${pid}.nc tmp2_pid${pid}.nc
+    
+    # remap
+    cdo remapbil,${dest_gridfile_u} tmp2_pid${pid}.nc tmp_pid${pid}.nc
+    cdo setmisstonn tmp_pid${pid}.nc tmp2_pid${pid}.nc
+
+    module purge
+    module load ELIC/0.1-foss-2018b-Python-3.6.6
+
+    input_full=`readlink -f tmp2_pid${pid}.nc`
+    output_full="${dest_folder}/${dest_prefix}grid_U_y${year}.nc"
+    
+    python ${reformat_script} ${input_full} u ${output_full}
+
+    module purge
+    module load CDO NCO    
+
+    
+    # 3) Extract v velocities.
+
+    if [ "${freq}" == "1d" ]; then
+	
+	tar -xvf ${in_prefix}/${name_exp}/NEMO/MMO_${name_exp}_18500101_fc0_${year}0101-${year}1231.tar ${name_exp}_${freq}_${year}0101_${year}1231_opa_grid_V_3D.nc ${name_exp}_${freq}_${year}0101_${year}1231_opa_grid_V_2D.nc
+	cdo selname,vbar ${name_exp}_${freq}_${year}0101_${year}1231_opa_grid_V_2D.nc tmp2_pid${pid}.nc
+
+    else
+	# vbar is not stored in monthly outputs, so take it from daily and monmean
+	tar -xvf ${in_prefix}/${name_exp}/NEMO/MMO_${name_exp}_18500101_fc0_${year}0101-${year}1231.tar ${name_exp}_${freq}_${year}0101_${year}1231_opa_grid_V_3D.nc ${name_exp}_1d_${year}0101_${year}1231_opa_grid_V_2D.nc
+	cdo selname,vbar ${name_exp}_1d_${year}0101_${year}1231_opa_grid_V_2D.nc tmp3_pid${pid}.nc
+	cdo monmean tmp3_pid${pid}.nc tmp2_pid${pid}.nc
+	
+    fi    
+    
+    cdo selname,vo ${name_exp}_${freq}_${year}0101_${year}1231_opa_grid_V_3D.nc tmp_pid${pid}.nc
+    cdo merge tmp_pid${pid}.nc tmp2_pid${pid}.nc tmp3_pid${pid}.nc
+
+    rm -f ${name_exp}_${freq}_${year}0101_${year}1231_opa_grid_V_3D.nc ${name_exp}_??_${year}0101_${year}1231_opa_grid_V_2D.nc
+    
+    # precut
+    cdo -O sellonlatbox,${latlon_cut} tmp3_pid${pid}.nc tmp_pid${pid}.nc
+
+    ncap2 -O -s 'vobaroclinic = vo - vbar;' tmp_pid${pid}.nc tmp2_pid${pid}.nc
+    ncks -x -v vo tmp2_pid${pid}.nc -O tmp_pid${pid}.nc
+    ncrename -v vbar,vobarotropic tmp_pid${pid}.nc
+    ncrename -v olevel,depthv -v olevel_bnds,depthv_bnds tmp_pid${pid}.nc
+    ncatted -a bounds,depthv,o,c,"depthv_bnds" -a name,depthv,o,c,"depthv" tmp_pid${pid}.nc
+    ncrename -d olevel,depthv tmp_pid${pid}.nc
+    ncpdq -O -a time_counter,depthv,y,x tmp_pid${pid}.nc tmp2_pid${pid}.nc
+    
+    # remap
+    cdo remapbil,${dest_gridfile_v} tmp2_pid${pid}.nc tmp_pid${pid}.nc
+    cdo setmisstonn tmp_pid${pid}.nc tmp2_pid${pid}.nc 
+
+    module purge
+    module load ELIC/0.1-foss-2018b-Python-3.6.6
+
+    input_full=`readlink -f tmp2_pid${pid}.nc`
+    output_full="${dest_folder}/${dest_prefix}grid_V_y${year}.nc"
+    
+    python ${reformat_script} ${input_full} v ${output_full}
+
+    module purge
+    module load CDO NCO    
+
+    
+done
+
+
+rm -rf ./*pid${pid}*

readme.md

@@ -1,3 +1,5 @@
+  - `anomalies`: calculating ocean temperature and salinity anomalies of EC-Earth vs. ORAS5.
+
   - `interp`: general guidelines and toolbox for performing horizontal and vertical interpolations with CDO.
 
   - `nemo_bdy`: generating ocean lateral boundary conditions for NEMO from EC-Earth (and also ORAS5).