Merge pull request #430 from wwieder/clm-diags

add CERES obs, SOILWATER, more regions

Author: wwieder

config_clm_unstructured_plots.yaml

@@ -140,29 +140,30 @@ diag_cam_climo:
     cam_overwrite_climo: false
 
     #Name of CAM case (or CAM run name):
-    cam_case_name: ctsm5.4.CMIP7_ciso_ctsm5.3.075_ne30_123_HIST_popDens
+    cam_case_name: b.e30_alpha07c_cesm.B1850C_LTso.ne30_t232_wgx3.247
 
     #Case nickname
     #NOTE: if nickname starts with '0' - nickname must be in quotes!
     # ie '026a' as opposed to 026a
     #If missing or left blank, will default to cam_case_name
-    case_nickname: '123_popDens'
+    case_nickname: 'B_247'
 
     #Location of CAM history (h0) files:
-    cam_hist_loc: /glade/derecho/scratch/wwieder/archive/${diag_cam_climo.cam_case_name}/lnd/hist/
+    cam_hist_loc: /glade/campaign/cesm/development/cross-wg/diagnostic_framework/CESM_output_for_testing/${diag_cam_climo.cam_case_name}/lnd/hist/
+    # /glade/derecho/scratch/wwieder/archive/${diag_cam_climo.cam_case_name}/lnd/hist/
 
     # If unstructured_plotting, a mesh file is required!
     mesh_file: /glade/campaign/cesm/cesmdata/inputdata/share/meshes/ne30pg3_ESMFmesh_cdf5_c20211018.nc
 
     #model year when time series files should start:
     #Note:  Leaving this entry blank will make time series
     #       start at earliest available year.
-    start_year: 1850
-    climo_start_year: 2004
+    start_year: 1
+    climo_start_year: 85
     #model year when time series files should end:
     #Note:  Leaving this entry blank will make time series
     #       end at latest available year.
-    end_year: 2023
+    end_year: 104
 
     #Do time series files exist?
     #If True, then diagnostics assumes that model files are already time series.
@@ -202,17 +203,18 @@ diag_cam_baseline_climo:
     cam_overwrite_climo: false
 
     #Name of CAM baseline case:
-    cam_case_name: ctsm5.4_5.3.068_PPEcal115_116_HIST
+    cam_case_name: b.e30_alpha07c_cesm.B1850C_LTso.ne30_t232_wgx3.234
 
     #Baseline case nickname
     #NOTE: if nickname starts with '0' - nickname must be in quotes!
     # ie '026a' as opposed to 026a
     #If missing or left blank, will default to cam_case_name
-    case_nickname: '116'
+    case_nickname: 'B_234'
 
     #Location of CAM baseline history (h0) files:
     #Example test files
-    cam_hist_loc: /glade/derecho/scratch/wwieder/archive/${diag_cam_baseline_climo.cam_case_name}/lnd/hist/
+    cam_hist_loc: /glade/campaign/cesm/development/cross-wg/diagnostic_framework/CESM_output_for_testing/${diag_cam_baseline_climo.cam_case_name}/lnd/hist/
+    # /glade/derecho/scratch/wwieder/archive/${diag_cam_baseline_climo.cam_case_name}/lnd/hist/
 
     # If unstructured_plotting, a mesh file is required!
     mesh_file: /glade/campaign/cesm/cesmdata/inputdata/share/meshes/ne30pg3_ESMFmesh_cdf5_c20211018.nc
@@ -223,13 +225,13 @@ diag_cam_baseline_climo:
     #model year when time series files should start:
     #Note:  Leaving this entry blank will make time series
     #       start at earliest available year.
-    start_year: 1850
-    climo_start_year: 2004
+    start_year: 1
+    climo_start_year: 61
 
     #model year when time series files should end:
     #Note:  Leaving this entry blank will make time series
     #       end at latest available year.
-    end_year: 2023
+    end_year: 80
 
     #Do time series files need to be generated?
     #If True, then diagnostics assumes that model files are already time series.
@@ -313,9 +315,10 @@ diag_var_list:
     - TOTECOSYSC
     - TOTSOMC_1m
     - ALTMAX
-    - FAREA_BURNED
     - TWS
+    - SOILWATER_10CM
     - GRAINC_TO_FOOD  
+    - FAREA_BURNED
     #- C13_GPP_pm
     # - C13_TOTVEGC_pm
     #- C14_GPP_pm
@@ -341,12 +344,14 @@ region_list:
     - Western US
     - Central US
     - Eastern US
+    - CONUS
     - Europe
     - Mediterranean
     - Central America
     - Amazonia
     - Central Africa
     - Indonesia
+    - S America
     - Brazil
     - Sahel
     - Southern Africa

lib/ldf_variable_defaults.yaml

@@ -122,12 +122,18 @@ FLDS: # atmospheric longwave radiation
       label : "W m$^{-2}$"
   pct_diff_contour_levels: [-100,-75,-50,-40,-30,-20,-10,-8,-6,-4,-2,0,2,4,6,8,10,20,30,40,50,75,100]
   pct_diff_colormap: "PuOr_r"
+  obs_file: "/glade/campaign/cgd/tss/people/oleson/FROM_LMWG/diag/lnd_diag4.2/obs_data/CERESed4.2-rlds_ALLMONS_climo.nc"
+  obs_name: "CERESed4.2"
+  obs_var_name: "rlds" #W/m2
 
 FSDS: # atmospheric incident solar radiation
   category: "Atmosphere"
   pct_diff_contour_levels: [-100,-75,-50,-40,-30,-20,-10,-8,-6,-4,-2,0,2,4,6,8,10,20,30,40,50,75,100]
   pct_diff_colormap: "PuOr_r"
   new_unit: "W m$^{-2}$"
+  obs_file: "/glade/campaign/cgd/tss/people/oleson/FROM_LMWG/diag/lnd_diag4.2/obs_data/CERESed4.2-rsds_ALLMONS_climo.nc"
+  obs_name: "CERESed4.2"
+  obs_var_name: "rsds" #W/m2
 
 WIND: # atmospheric air temperature
   category: "Atmosphere"
@@ -296,6 +302,10 @@ TOTRUNOFF: # total liquid runoff
     colorbar:
       label : "mm d$^{-1}$"
 
+SOILWATER_10CM: # soil liquid water + ice in top 10cm of soil 
+  category: "Hydrology"
+  colormap: "Blues"
+  
 TWS: # Terrestrial water storage
   category: "Hydrology"
   colormap: "Blues"

lib/regions_lnd.yaml

@@ -83,6 +83,10 @@ Eastern US:
   lat_bounds: [30, 50]
   lon_bounds: [-90, -70]
   region_category: Temperate 
+CONUS: 
+  lat_bounds: [25, 55]
+  lon_bounds: [-125, -70]
+  region_category: Temperate
 Europe:
   lat_bounds: [45, 60]
   lon_bounds: [-10, 30]
@@ -107,6 +111,10 @@ Indonesia:
   lat_bounds: [-10, 10]
   lon_bounds: [90, 150]
   region_category: Tropical 
+S America:
+  lat_bounds:  [-30, 5]
+  lon_bounds: [-85, -35]
+  region_category: Tropical 
 Brazil:
   lat_bounds: [-23.5, -10]
   lon_bounds: [-65, -30]

scripts/plotting/regional_climatology.py

@@ -78,7 +78,7 @@ def regional_climatology(adfobj):
     regional_climo_var_list = ['TSA','PREC','ELAI',
                                'FSDS','FLDS','SNOWDP','ASA',
                                'FSH','QRUNOFF_TO_COUPLER','ET','FCTR',
-                               'GPP','TWS','FCEV','FAREA_BURNED',
+                               'GPP','TWS','SOILWATER_10CM','FAREA_BURNED',
                                ]
 
     # Extract variables:
@@ -332,6 +332,10 @@ def regional_climatology(adfobj):
                     map_ax.set_extent([-180, 179, -89, 3],crs=ccrs.PlateCarree())
                 elif region_list[iReg]=='Polar':
                     map_ax.set_extent([-180, 179, 45, 90],crs=ccrs.PlateCarree())
+                elif region_list[iReg]=='CONUS':
+                    map_ax.set_extent([-140, -55, 10, 70],crs=ccrs.PlateCarree())
+                elif region_list[iReg]=='S America':
+                    map_ax.set_extent([-100, -20, -45, 20],crs=ccrs.PlateCarree())
                 else: 
                     if ((box_south >= 30) & (box_east<=-5) ):
                         map_ax.set_extent([-180, 0, 30, 90],crs=ccrs.PlateCarree())
@@ -358,17 +362,18 @@ def regional_climatology(adfobj):
                 print('Missing file for ', field)
                 continue
             else:
-                axs[plt_counter].plot(np.arange(12)+1, case_var_wgtd,
-                                      label=case_nickname, linewidth=2)
                 axs[plt_counter].plot(np.arange(12)+1, base_var_wgtd,
                                       label=base_nickname, linewidth=2)
+                axs[plt_counter].plot(np.arange(12)+1, case_var_wgtd,
+                                      label=case_nickname, linewidth=2)
                 if plot_obs[field] == True:
                     axs[plt_counter].plot(np.arange(12)+1, obs_var_wgtd,
                                           label=obs_name[field], color='black', linewidth=2)
+                    axs[plt_counter].legend()
+
                 axs[plt_counter].set_title(field)
                 axs[plt_counter].set_ylabel(base_data[field].units)
                 axs[plt_counter].set_xticks(np.arange(1, 13, 2))
-                axs[plt_counter].legend()
 
 
             plt_counter = plt_counter+1
@@ -437,6 +442,12 @@ def getRegion_xarray(varDS, varName,
 
     if varName not in varDS:
         varName = obs_var_name
+    
+    if varDS.lon.values.min() < 0:
+        # Convert lon to [0,360] if necessary
+        longitude = varDS['lon']
+        varDS = varDS.assign_coords(lon= (longitude + 180) % 360)
+        print(f"Converted lon to [0,360] for variable {varName}")
 
     # TODO is there a less brittle way to do this?
     if (area is not None) and (landfrac is not None):
@@ -455,7 +466,7 @@ def getRegion_xarray(varDS, varName,
             weight = weight * varDS['datamask']  
     else:
         raise ValueError("No valid weight, area, or landmask found in {varName} dataset.")
-    
+
     # check we have a data array for the variable
     if isinstance(varDS, xr.Dataset):
         varDS = varDS[varName]
@@ -478,7 +489,6 @@ def getRegion_xarray(varDS, varName,
                                    lon=(west_of_0 | east_of_0))
 
     wgt_subset = weight_subset / weight_subset.sum() 
-    weight_subset = weight.sel
     return domain_subset,wgt_subset
 
 def get_region_boundaries(regions, region_name):

scripts/plotting/regional_timeseries.py

@@ -79,15 +79,15 @@ def regional_timeseries(adfobj):
     regional_ts_var_list = ['TSA','PREC','ELAI',
                             'FSDS','FLDS','SNOWDP','ASA',
                             'FSH','QRUNOFF_TO_COUPLER','ET','FCTR',
-                            'GPP','TWS','FCEV','FAREA_BURNED',
+                            'GPP','TWS','SOILWATER_10CM','FAREA_BURNED',
                             ]
 
     # Extract variables:
-    baseline_name        = adfobj.get_baseline_info("cam_case_name", required=True)
-    input_ts_baseline = Path(adfobj.get_baseline_info("cam_ts_loc", required=True))
+    #baseline_name        = adfobj.get_baseline_info("cam_case_name", required=True)
+    #input_ts_baseline = Path(adfobj.get_baseline_info("cam_ts_loc", required=True))
     # TODO hard wired for single case name:
     case_name        = adfobj.get_cam_info("cam_case_name", required=True)[0]
-    input_ts_case = Path(adfobj.get_cam_info("cam_ts_loc", required=True)[0]) 
+    #input_ts_case = Path(adfobj.get_cam_info("cam_ts_loc", required=True)[0]) 
 
     # Get grid file 
     mesh_file = adfobj.mesh_files["baseline_mesh_file"]
@@ -132,7 +132,7 @@ def regional_timeseries(adfobj):
     obs_var_name = {}
     plot_obs = {}
 
-    var_obs_dict = adfobj.var_obs_dict
+    #var_obs_dict = adfobj.var_obs_dict
 
     # First, load all variable data once (instead of inside nested loops)
     for field in regional_ts_var_list:
@@ -333,6 +333,10 @@ def regional_timeseries(adfobj):
                     map_ax.set_extent([-180, 179, -89, 3],crs=ccrs.PlateCarree())
                 elif region_list[iReg]=='Polar':
                     map_ax.set_extent([-180, 179, 45, 90],crs=ccrs.PlateCarree())
+                elif region_list[iReg]=='CONUS':
+                    map_ax.set_extent([-140, -55, 10, 70],crs=ccrs.PlateCarree())
+                elif region_list[iReg]=='S America':
+                    map_ax.set_extent([-100, -20, -45, 20],crs=ccrs.PlateCarree())
                 else: 
                     if ((box_south >= 30) & (box_east<=-5) ):
                         map_ax.set_extent([-180, 0, 30, 90],crs=ccrs.PlateCarree())
@@ -370,7 +374,8 @@ def regional_timeseries(adfobj):
                 #                          label=obs_name[field], color='black', linewidth=2)
                 axs[plt_counter].set_title(field)
                 axs[plt_counter].set_ylabel(base_data[field].units)
-                axs[plt_counter].legend()                
+                if plt_counter == 3 or plt_counter == 7 or plt_counter ==11 or plt_counter ==15:
+                    axs[plt_counter].legend()                
 
             plt_counter = plt_counter+1