Adding geochemical tracer evaluations.

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,6 +1,6 @@
 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'], 
-             maskfixes=True, use_for_eval=False, complexity='boxes', verbose=False, biasmaps=False):
+             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, obs=None):
     '''
     AUTHORS:
     Jan Streffing		2022-12-01	Split off from main tool
@@ -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()
-
-    obs = [siconc, tas, clt, pr, rlut, uas, vas, ua, zg, zos, tos, mlotst, thetao, so]
+    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()
+    if obs == None:
+        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')
 

cmpitool/read_errors.py

@@ -42,7 +42,7 @@ def read_errors(obs, eval_models, regions, seasons, out_path, eval_path, n_imple
     regions_names = []
     for region in regions:
         regions_names.append(region.name)
-        
+
     collect = np.empty([len(eval_models),n_implemented_var,len(regions),max_depth,len(seasons)])*np.nan
     i=0
     for eval_model in tqdm(eval_models):
@@ -68,18 +68,22 @@ def read_errors(obs, eval_models, regions, seasons, out_path, eval_path, n_imple
                     m=0
                     for seas in seasons:
                         if any(a): # Check if variable appears in csv. If not, skip it.
-                            if regions_csv[r] not in regions_names: # Check if region from csv part of the analysis. Else advance
+                            if season_csv[r] not in seasons: # Check if region from csv part of the analysis. Else advance
                                 while True:
+                                    if verbose:
+                                        print('skipping season',season_csv[r],r,m,l,k,j,i)
                                     r+=1
-                                    if regions_csv[r] in regions_names:
+                                    if season_csv[r] in seasons:
                                         break
-                            if season_csv[r] not in seasons: # Check if region from csv part of the analysis. Else advance
+                            if regions_csv[r] not in regions_names: # Check if region from csv part of the analysis. Else advance
                                 while True:
+                                    if verbose:
+                                        print('skipping region',regions_csv[r],r,m,l,k,j,i)
                                     r+=1
-                                    if season_csv[r] in seasons:
+                                    if regions_csv[r] in regions_names:
                                         break
                             if verbose:
-                                print(eval_model.name,var.name,region.name,depth,seas,values[r])
+                                print('reading',r,m,l,k,j,i,eval_model.name,var.name,region.name,depth,seas,values[r])
                             collect[i,j,k,l,m]=values[r]
                             r+=1
                         m+=1