Additions for geochemical tracers

cmpitool/cmpisetup.py

@@ -29,6 +29,7 @@ def cmpisetup(reanalysis='ERA5'):
     mlotst                      Variable object for mixed layer depth
     thetao                      Variable object for 3D ocean temperature
     so                          Variable object for 3D ocean salinity
+    o2                          Variable object for 3D oxygen in ocean 
     '''
 
 
@@ -71,6 +72,13 @@ def __init__ (self, name, variables):
     mlotst = variable(name='mlotst', obs='C-GLORSv7', depths=['surface'], domain='oce')
     thetao = variable(name='thetao', obs='EN4', depths=['10m','100m','1000m'], domain='oce')
     so = variable(name='so', obs='EN4', depths=['10m','100m','1000m'], domain='oce')
-
-    return variable, region, climate_model, siconc, tas, clt, pr, rlut, uas, vas, ua, zg, zos, tos, mlotst, thetao, so
+    alk = variable(name='alk', obs='GLODAP', depths=['50m','400m','800m','2000m'], domain='oce')
+    co2 = variable(name='co2', obs='GLODAP', depths=['50m','400m','800m','2000m'], domain='oce')
+    n = variable(name='n', obs='WOA23', depths=['50m','400m','800m'], domain='oce')
+    o2 = variable(name='o2', obs='WOA23', depths=['50m','400m','1000m'], domain='oce')
+    p = variable(name='p', obs='WOA23', depths=['50m','400m','800m'], domain='oce')
+    zoo = variable(name='zoo', obs='MOPS', depths=['50m'], domain='oce')
+    phy = variable(name='phy', obs='MOPS', depths=['50m'], domain='oce')
+
+    return variable, region, climate_model, siconc, tas, clt, pr, rlut, uas, vas, ua, zg, zos, tos, mlotst, thetao, so, alk, co2, n, o2, p, zoo, phy
 

cmpitool/cmpitool.py

@@ -1,5 +1,5 @@
 def cmpitool(model_path, models, eval_models=None, out_path='output/', obs_path='obs/' , reanalysis='ERA5', 
-             eval_path=None, time='198912-201411', seasons=['MAM', 'JJA', 'SON', 'DJF'], 
+             eval_path=None, time='198912-201411', seasons=['MAM', 'JJA', 'SON', 'DJF', 'year'], # time changes from time='198912-201411'
              maskfixes=True, use_for_eval=False, complexity='boxes', verbose=False, biasmaps=False):
     '''
     AUTHORS:
@@ -12,7 +12,7 @@ def cmpitool(model_path, models, eval_models=None, out_path='output/', obs_path=
     INPUT:
     model_path                  Path pointing towards the output of your model,
                                 preprocessed to be read in by cmiptool
-    models		        List of climate model objects to be evaluated via cmiptool
+    models		                List of climate model objects to be evaluated via cmiptool
     eval_models                 List of climate model objects used as reference for evaluation
                                 By default this is set to None, which results in a set of 30 CMIP6
                                 being used
@@ -35,7 +35,6 @@ def cmpitool(model_path, models, eval_models=None, out_path='output/', obs_path=
                                 for simple lat/lon boxes (boxes) or continents & ocean
                                 basins (regions)
     verbose                     Boolean to activate verbose output
-    biasmaps                    Boolean to activate bias map plots
 
 
     RETURN:
@@ -53,9 +52,9 @@ def cmpitool(model_path, models, eval_models=None, out_path='output/', obs_path=
     else:
         eval_path=eval_path+'/'
 
-    variable, region, climate_model, siconc, tas, clt, pr, rlut, uas, vas, ua, zg, zos, tos, mlotst, thetao, so = cmpisetup()
+    variable, region, climate_model, siconc, tas, clt, pr, rlut, uas, vas, ua, zg, zos, tos, mlotst, thetao, so, alk, co2, n, o2, p, zoo, phy = cmpisetup()
 
-    obs = [siconc, tas, clt, pr, rlut, uas, vas, ua, zg, zos, tos, mlotst, thetao, so]
+    obs = [siconc, tas, clt, pr, rlut, uas, vas, ua, zg, zos, tos, mlotst, thetao, so, alk, co2, n, o2, p, zoo, phy]
 
     '''
     If you don't add all variables to obs for your analysis, the missing ones will be skipped.
@@ -65,7 +64,7 @@ def cmpitool(model_path, models, eval_models=None, out_path='output/', obs_path=
     - If you add more variables and generate new .csv files, increase the number 14 accordingly!
     - If you just skip a variable for your analysis, don't change number_of_implemented_variables!
     '''
-    n_implemented_var = 14 
+    n_implemented_var = 21
 
     #The CMIP6 models are set up by default in their own function
     cmip6_models = config_cmip6(climate_model, obs)
@@ -150,6 +149,6 @@ def cmpitool(model_path, models, eval_models=None, out_path='output/', obs_path=
     cmpi =  write_fractions(error_fraction, models, regions, seasons, out_path, verbose)
 
     plotting_heatmaps(models, regions, seasons, obs, error_fraction, cmpi, out_path, verbose)
-    
+
     if biasmaps == True:
         plotting_biasmaps(ds_model, ds_obs , models, seasons, obs, out_path, verbose)

cmpitool/loading_models.py

@@ -28,7 +28,12 @@ def loading_models(models, model_path, seasons, time, verbose):
     for model in tqdm(models):
         for var in model.variables:
             for depth in var.depths:
-                for seas in seasons:
+                if var.name in ('alk', 'co2', 'zoo', 'phy'):
+                    season = ['year']
+                else:
+                    season = seasons
+
+                for seas in season:
                     if verbose:
                         print('loading '+model_path+var.name+'_'+model.name+'_'+time+'_'+depth+'_'+seas+'.nc')
                     intermediate = xr.open_mfdataset(model_path+var.name+'_'+model.name+'_'+time+'_'+depth+'_'+seas+'.nc')

cmpitool/loading_obs.py

@@ -32,7 +32,12 @@ def loading_obs(obs, obs_path, seasons, verbose):
 
     for var in tqdm(obs):
         for depth in var.depths:
-            for seas in seasons:
+            if var.name in ('alk', 'co2', 'zoo', 'phy'):
+                season = ['year']
+            else:
+                season = seasons
+
+            for seas in season:
                 if verbose:
                     print('loading '+obs_path+var.name+'_'+var.obs+'_'+depth+'_'+seas+'.nc')
 

preprocessing_examples/Preprocessing_geochemical.sh

@@ -0,0 +1,160 @@
+#!/bin/bash
+
+Help()
+{
+   # Display Help
+	echo "##############################################################################"
+	echo "# This is an example script preparing climate output without CMOR complient  #"
+	echo "# output from the cmip-tool                                                  #"
+	echo "# Author: Jan Streffing/Lina Lucas 2024-03-10					   #"
+	echo "##############################################################################"
+	echo "Positional arguments:"
+	echo "#1    The original path to folder of ocean and atmosphere model data"  
+	echo "#2    Path where the Pre-preprocessed data should be stored"
+	echo "#3    Path where preprocessed data should be stored"
+	echo "#4    Name of the climate model (your free choice!)"
+	echo "#5    Start year of analysed period"
+	echo "#6    End year of analysed period"
+	echo "Positional optional argument:"
+	echo "#7    Set to true when temporary data (tmpstr) should be deleted"
+	echo "#################################################"
+	echo "# example: ./Preprocessing_geochemical.sh /home/llucas /home/llucas/CMPI_data/BM /home/llucas/cmpitool/preprocessed_data FOCI1.20.0-CC104_RCP_ESM_spinup2089 2009 2010 true" 
+}
+
+# Check if the script was called with the help option
+if [[ $1 == "--help" ]]; then
+  Help
+  exit 0  # Exit successfully after displaying help
+fi
+
+printf "##########################\n"
+printf "# read command line args #\n"
+printf "##########################\n"
+origdir=$1    # The original path to folder of ocean and atmosphere model data
+predata=$2    # Path where the Pre-preprocessed data should be stored
+outdir=$3     # Path where preprocessed data should be stored
+model_name=$4 # Name of the climate model
+starty=$5     # Start year of analysed period
+endy=$6       # End year of analysed period
+deltmp=${7}   # Set to true when temporary data (tmpstr) should be deleted
+
+
+# Create folder of predata and go into the path of the Pre-preprocessed data
+mkdir -p $predata
+mkdir -p $outdir
+cd $predata
+tmpstr="analysis_cmpi_period"
+
+
+# Delete files in the predata and out directory so the code can work properly
+printf "##############################################\n"
+printf "# clean up so we can create new files #\n"
+printf "##############################################\n"
+
+rm -f ${predata}/* 
+rm -f ${predata}/*
+rm -f ${outdir}/*${model_name}_${starty}-${endy}*
+wait
+
+
+printf "#################\n"
+printf "# Loading files and selecting variables and timeperiod #\n"
+printf "#################\n"
+
+# Atmosphere 1
+for var in ALK DIC DIN O2 PO4 ZOO PHY;
+do
+    cdo -select,startdate=${starty}-01-01T00:00:00,enddate=${endy}-01-01T00:00:00,name=${var} ${origdir}/${model_name}_1m_*_ptrc_T.nc ${var}_${starty}-${endy}.nc 
+done
+wait
+
+
+cd $outdir
+
+
+printf "#################\n"
+printf "# ALK and DIC #\n"
+printf "#################\n"
+
+# Interpolating levels
+for var in ALK DIC;
+do
+    cdo -L -intlevel,50,400,800,2000 -setctomiss,0 ${predata}/${var}_${starty}-${endy}.nc ${var}_${tmpstr}_int.nc 
+done
+wait 
+
+# Split levels, change names, remap
+cdo -L -splitlevel -chname,ALK,alk -remapbil,r180x91 ALK_${tmpstr}_int.nc alk_${tmpstr}_remap_
+cdo -L -splitlevel -chname,DIC,co2 -remapbil,r180x91 DIC_${tmpstr}_int.nc co2_${tmpstr}_remap_
+wait
+
+# Adjust names
+for lvl in 50 400 800 2000;
+do
+    mv -f alk_${tmpstr}_remap_$(printf "%06d" $lvl).nc alk_${model_name}_${starty}-${endy}_${lvl}m_year.nc
+    mv -f co2_${tmpstr}_remap_$(printf "%06d" $lvl).nc co2_${model_name}_${starty}-${endy}_${lvl}m_year.nc
+done
+
+
+
+printf "#################\n"
+printf "# DIN, O2, PO4 for seasonal means #\n"
+printf "#################\n"
+
+# Interpolating levels
+cdo -L -intlevel,50,400,800 -setctomiss,0 ${predata}/DIN_${starty}-${endy}.nc DIN_${tmpstr}_int.nc 
+cdo -L -intlevel,50,400,800 -setctomiss,0 ${predata}/PO4_${starty}-${endy}.nc PO4_${tmpstr}_int.nc 
+cdo -L -intlevel,50,400,1000 -setctomiss,0 ${predata}/O2_${starty}-${endy}.nc O2_${tmpstr}_int.nc 
+wait
+
+# Split levels, change names, remap
+cdo -L -splitlevel -chname,DIN,n -remapbil,r180x91 DIN_${tmpstr}_int.nc n_${tmpstr}_remap_
+cdo -L -splitlevel -chname,PO4,p -remapbil,r180x91 PO4_${tmpstr}_int.nc p_${tmpstr}_remap_
+cdo -L -splitlevel -chname,O2,o2 -remapbil,r180x91 O2_${tmpstr}_int.nc o2_${tmpstr}_remap_
+wait
+
+# Split seasons
+for lvl in 50 400 800;
+do
+    cdo -L splitseas -yseasmean n_${tmpstr}_remap_$(printf "%06d" $lvl).nc n_${model_name}_${starty}-${endy}_${lvl}m_ &
+    cdo -L splitseas -yseasmean p_${tmpstr}_remap_$(printf "%06d" $lvl).nc p_${model_name}_${starty}-${endy}_${lvl}m_ &
+done
+
+for lvl in 50 400 1000;
+do
+    cdo -L splitseas -yseasmean o2_${tmpstr}_remap_$(printf "%06d" $lvl).nc o2_${model_name}_${starty}-${endy}_${lvl}m_ &
+done
+wait
+
+
+printf "#################\n"
+printf "# ZOO and PHY #\n"
+printf "#################\n"
+
+# Interpolating levels
+for var in ZOO PHY;
+do
+    cdo -L -intlevel,50 -setctomiss,0 ${predata}/${var}_${starty}-${endy}.nc ${var}_${tmpstr}_int.nc 
+done
+wait 
+
+# Split levels, change names, remap
+cdo -L -splitlevel -chname,ZOO,zoo -remapbil,r180x91 ZOO_${tmpstr}_int.nc zoo_${tmpstr}_remap_
+cdo -L -splitlevel -chname,PHY,phy -remapbil,r180x91 PHY_${tmpstr}_int.nc phy_${tmpstr}_remap_
+wait
+
+# Adjust names
+for lvl in 50;
+do
+    mv -f zoo_${tmpstr}_remap_$(printf "%06d" $lvl).nc zoo_${model_name}_${starty}-${endy}_${lvl}m_year.nc
+    mv -f phy_${tmpstr}_remap_$(printf "%06d" $lvl).nc phy_${model_name}_${starty}-${endy}_${lvl}m_year.nc
+done
+
+
+
+# If deltmp set to true, deleting all the temporary files
+if $deltmp; then
+	printf "Deleting tmp data\n"
+	rm -rf *${tmpstr}*
+	rm -rf *_int*
+fi