diff --git a/.gitignore b/.gitignore index 72364f99f..42b70ac71 100644 --- a/.gitignore +++ b/.gitignore @@ -1,3 +1,7 @@ +# User defined config files +/config.analysis.* +/config.analysis_* + # Byte-compiled / optimized / DLL files __pycache__/ *.py[cod] diff --git a/config.analysis b/config.analysis index c39304248..e58bb0146 100644 --- a/config.analysis +++ b/config.analysis @@ -3,7 +3,7 @@ casename = 20160805v0atm.A_WCYCL1850_v0atm.ne30_oEC.edison.alpha7_00 native_res = ne30 short_term_archive = 0 -casename_model_tocompare = B1850C5_ne30_v0.4 +ref_casename_v0 = B1850C5_ne30_v0.4 [paths] # paths to simulation and observational datasets @@ -15,14 +15,14 @@ log_dir = /global/project/projectdirs/acme/xylar/coupled_diagnostics_20160805v0a obs_ocndir = /global/project/projectdirs/acme/observations/Ocean obs_sstdir = /global/project/projectdirs/acme/observations/Ocean/SST obs_seaicedir = /global/project/projectdirs/acme/observations/SeaIce -ocndir_model_tocompare = /global/project/projectdirs/acme/ACMEv0_lowres/B1850C5_ne30_v0.4/ocn/postprocessing -seaicedir_model_tocompare = /global/project/projectdirs/acme/ACMEv0_lowres/B1850C5_ne30_v0.4/ice/postprocessing +ref_archive_v0_ocndir = /global/project/projectdirs/acme/ACMEv0_lowres/B1850C5_ne30_v0.4/ocn/postprocessing +ref_archive_v0_seaicedir = /global/project/projectdirs/acme/ACMEv0_lowres/B1850C5_ne30_v0.4/ice/postprocessing [data] # paths to mesh and mapping files mpas_meshfile = /global/project/projectdirs/acme/milena/MPAS-grids/ocn/gridfile.oEC60to30.nc mpas_remapfile = /global/project/projectdirs/acme/mapping/maps/map_oEC60to30_TO_0.5x0.5degree_blin.160412.nc -model_tocompare_remapfile = /global/project/projectdirs/acme/mapping/maps/map_gx1v6_TO_0.5x0.5degree_blin.160413.nc +pop_remapfile = /global/project/projectdirs/acme/mapping/maps/map_gx1v6_TO_0.5x0.5degree_blin.160413.nc # path to climotology dataset mpas_climodir = /global/project/projectdirs/acme/xylar/20160805v0atm.A_WCYCL1850_v0atm.ne30_oEC.edison.alpha7_00.test.pp @@ -66,8 +66,8 @@ generate = 1 generate = 1 [ohc_timeseries] -# compare to output from another model run? -compare_with_model = True +## compare to output from another model run? +#compare_with_model = True # compare to observations? compare_with_obs = False # list of region indices to plot from the region list below @@ -77,8 +77,8 @@ regionIndicesToPlot = [6] N_movavg = 12 [sst_timeseries] -# compare to output from another model run? -compare_with_model = True +## compare to output from another model run? +#compare_with_model = True # compare to observations? compare_with_obs = True # list of region indices to plot from the region list below @@ -88,8 +88,8 @@ regionIndicesToPlot = [6] N_movavg = 12 [nino34_timeseries] -# compare to output from another model run? -compare_with_model = True +## compare to output from another model run? +#compare_with_model = True # compare to observations? compare_with_obs = True # Number of points over which to compute moving average (e.g., for monthly @@ -97,8 +97,8 @@ compare_with_obs = True N_movavg = 12 [mht_timeseries] -# compare to output from another model run? -compare_with_model = True +## compare to output from another model run? +#compare_with_model = True # compare to observations? compare_with_obs = True # Number of points over which to compute moving average (e.g., for monthly @@ -106,8 +106,8 @@ compare_with_obs = True N_movavg = 12 [moc_timeseries] -# compare to output from another model run? -compare_with_model = True +## compare to output from another model run? +#compare_with_model = True # compare to observations? compare_with_obs = True # Number of points over which to compute moving average (e.g., for monthly @@ -115,8 +115,8 @@ compare_with_obs = True N_movavg = 12 [seaice_timeseries] -# compare to output from another model run? -compare_with_model = True +## compare to output from another model run? +#compare_with_model = True # compare to observations? compare_with_obs = True # Number of points over which to compute moving average (e.g., for monthly diff --git a/mpas_analysis/ocean/ohc_timeseries.py b/mpas_analysis/ocean/ohc_timeseries.py index 15d479dd1..c476b54a6 100644 --- a/mpas_analysis/ocean/ohc_timeseries.py +++ b/mpas_analysis/ocean/ohc_timeseries.py @@ -14,10 +14,9 @@ def ohc_timeseries(config): indir = config.get('paths','archive_dir_ocn') meshfile = config.get('data','mpas_meshfile') casename = config.get('case','casename') - casename_model_tocompare = config.get('case','casename_model_tocompare') - indir_model_tocompare = config.get('paths','ocndir_model_tocompare') + ref_casename_v0 = config.get('case','ref_casename_v0') + indir_v0data = config.get('paths','ref_archive_v0_ocndir') - compare_with_model = config.getboolean('ohc_timeseries','compare_with_model') compare_with_obs = config.getboolean('ohc_timeseries','compare_with_obs') plots_dir = config.get('paths','plots_dir') @@ -77,19 +76,18 @@ def ohc_timeseries(config): time_start = datetime.datetime(year_start,1,1) time_end = datetime.datetime(year_end,12,31) - if compare_with_model: - print " Load OHC for model_to_compare and make plots..." - # load in other model run data - infiles_model_tocompare = "".join([indir_model_tocompare,'/OHC.',casename_model_tocompare,'.year*.nc']) - ds_model_tocompare = xr.open_mfdataset(infiles_model_tocompare,preprocess=lambda x: preprocess_mpas(x, yearoffset=yr_offset)) - ds_model_tocompare = remove_repeated_time_index(ds_model_tocompare) - ds_model_tocompare_tslice = ds_model_tocompare.sel(Time=slice(time_start,time_end)) + if ref_casename_v0 != "None": + print " Load in OHC for ACMEv0 case..." + infiles_v0data = "".join([indir_v0data,'/OHC.',ref_casename_v0,'.year*.nc']) + ds_v0 = xr.open_mfdataset(infiles_v0data,preprocess=lambda x: preprocess_mpas(x, yearoffset=yr_offset)) + ds_v0 = remove_repeated_time_index(ds_v0) + ds_v0_tslice = ds_v0.sel(Time=slice(time_start,time_end)) sumLayerMaskValue = ds.time_avg_avgValueWithinOceanLayerRegion_sumLayerMaskValue avgLayerArea = ds.time_avg_avgValueWithinOceanLayerRegion_avgLayerArea avgLayerThickness = ds.time_avg_avgValueWithinOceanLayerRegion_avgLayerThickness - print " Compute OHC..." + print " Compute OHC and make plots..." for index in range(len(iregions)): iregion = iregions[index] @@ -118,25 +116,22 @@ def ohc_timeseries(config): xlabel = "Time [years]" ylabel = "[x$10^{22}$ J]" - if compare_with_model: - figname = "%s/ohc_%s_%s_%s.png" % (plots_dir,regions[iregion],casename,casename_model_tocompare) - # load in other model run data - ohc_model_tocompare_tot = ds_model_tocompare_tslice.ohc_tot - ohc_model_tocompare_700m = ds_model_tocompare_tslice.ohc_700m - ohc_model_tocompare_2000m = ds_model_tocompare_tslice.ohc_2000m - ohc_model_tocompare_btm = ds_model_tocompare_tslice.ohc_btm - title = "".join([title," (r), ",casename_model_tocompare," (b)"]) + if ref_casename_v0 != "None": + figname = "%s/ohc_%s_%s_%s.png" % (plots_dir,regions[iregion],casename,ref_casename_v0) + ohc_v0_tot = ds_v0_tslice.ohc_tot + ohc_v0_700m = ds_v0_tslice.ohc_700m + ohc_v0_2000m = ds_v0_tslice.ohc_2000m + ohc_v0_btm = ds_v0_tslice.ohc_btm + title = "".join([title," (r), ",ref_casename_v0," (b)"]) timeseries_analysis_plot(config, [ohc_tot, ohc_700m, ohc_2000m, ohc_btm, - ohc_model_tocompare_tot, ohc_model_tocompare_700m, - ohc_model_tocompare_2000m, ohc_model_tocompare_btm], + ohc_v0_tot, ohc_v0_700m, ohc_v0_2000m, ohc_v0_btm], N_movavg, title, xlabel, ylabel, figname, lineStyles = ['r-', 'r-', 'r--', 'r-.', 'b-', 'b-', 'b--', 'b-.'], lineWidths = [2, 1, 1.5, 1.5, 2, 1, 1.5, 1.5]) - if not compare_with_obs and not compare_with_model: + if not compare_with_obs and ref_casename_v0 == "None": figname = "%s/ohc_%s_%s.png" % (plots_dir,regions[iregion],casename) timeseries_analysis_plot(config, [ohc_tot, ohc_700m, ohc_2000m, ohc_btm], N_movavg, title, xlabel, ylabel, figname, lineStyles = ['r-', 'r-', 'r--', 'r-.'], lineWidths = [2, 1, 1.5, 1.5]) - diff --git a/mpas_analysis/ocean/sst_timeseries.py b/mpas_analysis/ocean/sst_timeseries.py index a632dc2cc..abd1da750 100644 --- a/mpas_analysis/ocean/sst_timeseries.py +++ b/mpas_analysis/ocean/sst_timeseries.py @@ -14,10 +14,9 @@ def sst_timeseries(config): infiles = '%s/am.mpas-o.timeSeriesStats.????-??*nc'%indir casename = config.get('case','casename') - casename_model_tocompare = config.get('case','casename_model_tocompare') - indir_model_tocompare = config.get('paths','ocndir_model_tocompare') + ref_casename_v0 = config.get('case','ref_casename_v0') + indir_v0data = config.get('paths','ref_archive_v0_ocndir') - compare_with_model = config.getboolean('sst_timeseries','compare_with_model') compare_with_obs = config.getboolean('sst_timeseries','compare_with_obs') plots_dir = config.get('paths','plots_dir') @@ -34,10 +33,8 @@ def sst_timeseries(config): ds = xr.open_mfdataset(infiles, preprocess=lambda x: preprocess_mpas(x, yearoffset=yr_offset, timeSeriesStats=True, timestr='time_avg_daysSinceStartOfSim', onlyvars=['time_avg_avgValueWithinOceanRegion_avgSurfaceTemperature'])) - ds = remove_repeated_time_index(ds) - SSTregions = ds.time_avg_avgValueWithinOceanRegion_avgSurfaceTemperature year_start = (pd.to_datetime(ds.Time.min().values)).year @@ -46,17 +43,14 @@ def sst_timeseries(config): time_end = datetime.datetime(year_end,12,31) - if compare_with_model: - # load in other model run data - infiles_model_tocompare = "".join([indir_model_tocompare,'/SST.',casename_model_tocompare,'.year*.nc']) - ds_model_tocompare = xr.open_mfdataset(infiles_model_tocompare,preprocess=lambda x: preprocess_mpas(x, yearoffset=yr_offset)) - ds_model_tocompare = remove_repeated_time_index(ds_model_tocompare) - ds_model_tocompare_tslice = ds_model_tocompare.sel(Time=slice(time_start,time_end)) + if ref_casename_v0 != "None": + print " Load in SST for ACMEv0 case..." + infiles_v0data = "".join([indir_v0data,'/SST.',ref_casename_v0,'.year*.nc']) + ds_v0 = xr.open_mfdataset(infiles_v0data,preprocess=lambda x: preprocess_mpas(x, yearoffset=yr_offset)) + ds_v0 = remove_repeated_time_index(ds_v0) + ds_v0_tslice = ds_v0.sel(Time=slice(time_start,time_end)) - - # Load data and make plot for every region print " Make plots..." - for index in range(len(iregions)): iregion = iregions[index] @@ -67,13 +61,12 @@ def sst_timeseries(config): SST = SSTregions[:,iregion] - if compare_with_model: - # load in other model run data - figname = "%s/sst_%s_%s_%s.png" % (plots_dir,regions[iregion],casename,casename_model_tocompare) - SST_model_tocompare = ds_model_tocompare_tslice.SST + if ref_casename_v0 != "None": + figname = "%s/sst_%s_%s_%s.png" % (plots_dir,regions[iregion],casename,ref_casename_v0) + SST_v0 = ds_v0_tslice.SST - title = "%s\n %s (b-)" % (title, casename_model_tocompare) - timeseries_analysis_plot(config, [SST,SST_model_tocompare], N_movavg, + title = "%s\n %s (b-)" % (title, ref_casename_v0) + timeseries_analysis_plot(config, [SST,SST_v0], N_movavg, title, xlabel, ylabel, figname, lineStyles = ['r-','b-'], lineWidths = [1.2,1.2]) @@ -81,6 +74,3 @@ def sst_timeseries(config): figname = "%s/sst_%s_%s.png" % (plots_dir,regions[iregion],casename) timeseries_analysis_plot(config, [SST], N_movavg, title, xlabel, ylabel, figname, lineStyles = ['r-'], lineWidths = [1.2]) - - - diff --git a/mpas_analysis/sea_ice/timeseries.py b/mpas_analysis/sea_ice/timeseries.py index 30ee9db3b..6209731cb 100644 --- a/mpas_analysis/sea_ice/timeseries.py +++ b/mpas_analysis/sea_ice/timeseries.py @@ -32,10 +32,9 @@ def seaice_timeseries(config): meshfile = config.get('data','mpas_meshfile') casename = config.get('case','casename') - casename_model_tocompare = config.get('case','casename_model_tocompare') - indir_model_tocompare = config.get('paths','seaicedir_model_tocompare') + ref_casename_v0 = config.get('case','ref_casename_v0') + indir_v0data = config.get('paths','ref_archive_v0_seaicedir') - compare_with_model = config.getboolean('seaice_timeseries','compare_with_model') compare_with_obs = config.getboolean('seaice_timeseries','compare_with_obs') plots_dir = config.get('paths','plots_dir') @@ -44,7 +43,6 @@ def seaice_timeseries(config): N_movavg = config.getint('seaice_timeseries','N_movavg') - print " Load sea-ice data..." # Load mesh dsmesh = xr.open_dataset(meshfile) @@ -65,12 +63,11 @@ def seaice_timeseries(config): time_start = datetime.datetime(year_start,1,1) time_end = datetime.datetime(year_end,12,31) - if compare_with_model: - infiles_model_tocompare = "".join([indir_model_tocompare,'/icevol.',casename_model_tocompare,'.year*.nc']) - ds_model_tocompare = xr.open_mfdataset(infiles_model_tocompare, + if ref_casename_v0 != "None": + infiles_v0data = "".join([indir_v0data,'/icevol.',ref_casename_v0,'.year*.nc']) + ds_v0 = xr.open_mfdataset(infiles_v0data, preprocess=lambda x: preprocess_mpas(x, yearoffset=yr_offset)) - ds_model_tocompare_tslice = ds_model_tocompare.sel(Time=slice(time_start,time_end)) - + ds_v0_tslice = ds_v0.sel(Time=slice(time_start,time_end)) # Make Northern and Southern Hemisphere partition: areaCell = ds.areaCell @@ -119,9 +116,9 @@ def seaice_timeseries(config): xlabel = "Time [years]" - if compare_with_model: - figname_nh = "%s/%sNH_%s_%s.png" % (plots_dir,varname,casename,casename_model_tocompare) - figname_sh = "%s/%sSH_%s_%s.png" % (plots_dir,varname,casename,casename_model_tocompare) + if ref_casename_v0 != "None": + figname_nh = "%s/%sNH_%s_%s.png" % (plots_dir,varname,casename,ref_casename_v0) + figname_sh = "%s/%sSH_%s_%s.png" % (plots_dir,varname,casename,ref_casename_v0) else: figname_nh = "%s/%sNH_%s.png" % (plots_dir,varname,casename) figname_sh = "%s/%sSH_%s.png" % (plots_dir,varname,casename) @@ -137,9 +134,9 @@ def seaice_timeseries(config): title_nh = "%s\nPIOMAS, annual cycle (k)" % title_nh title_sh = "%s\n" % title_sh - if compare_with_model: - title_nh = "%s\n %s (b)" % (title_nh,casename_model_tocompare) - title_sh = "%s\n %s (b)" % (title_sh,casename_model_tocompare) + if ref_casename_v0 != "None": + title_nh = "%s\n %s (b)" % (title_nh,ref_casename_v0) + title_sh = "%s\n %s (b)" % (title_sh,ref_casename_v0) if varname == "iceAreaCell": @@ -157,18 +154,17 @@ def seaice_timeseries(config): var_sh_obs = ds_obs.IceArea var_sh_obs = replicate_cycle(var_sh,var_sh_obs) - if compare_with_model: - infiles_model_tocompare = "".join([indir_model_tocompare,'/icearea.',casename_model_tocompare,'.year*.nc']) - ds_model_tocompare = xr.open_mfdataset(infiles_model_tocompare, + if ref_casename_v0 != "None": + infiles_v0data = "".join([indir_v0data,'/icearea.',ref_casename_v0,'.year*.nc']) + ds_v0 = xr.open_mfdataset(infiles_v0data, preprocess=lambda x: preprocess_mpas(x, yearoffset=yr_offset)) - ds_model_tocompare_tslice = ds_model_tocompare.sel(Time=slice(time_start,time_end)) - var_nh_model_tocompare = ds_model_tocompare_tslice.icearea_nh - var_sh_model_tocompare = ds_model_tocompare_tslice.icearea_sh + ds_v0_tslice = ds_v0.sel(Time=slice(time_start,time_end)) + var_nh_v0 = ds_v0_tslice.icearea_nh + var_sh_v0 = ds_v0_tslice.icearea_sh elif varname == "iceVolumeCell": if compare_with_obs: - ds_obs = xr.open_mfdataset(obs_filenameNH, preprocess=lambda x: preprocess_mpas(x, yearoffset=yr_offset)) ds_obs = remove_repeated_time_index(ds_obs) @@ -177,36 +173,35 @@ def seaice_timeseries(config): var_sh_obs = None - if compare_with_model: - - infiles_model_tocompare = "".join([indir_model_tocompare,'/icevol.',casename_model_tocompare,'.year*.nc']) - ds_model_tocompare = xr.open_mfdataset(infiles_model_tocompare, + if ref_casename_v0 != "None": + infiles_v0data = "".join([indir_v0data,'/icevol.',ref_casename_v0,'.year*.nc']) + ds_v0 = xr.open_mfdataset(infiles_v0data, preprocess=lambda x: preprocess_mpas(x, yearoffset=yr_offset)) - ds_model_tocompare_tslice = ds_model_tocompare.sel(Time=slice(time_start,time_end)) - var_nh_model_tocompare = ds_model_tocompare_tslice.icevolume_nh - var_sh_model_tocompare = ds_model_tocompare_tslice.icevolume_sh + ds_v0_tslice = ds_v0.sel(Time=slice(time_start,time_end)) + var_nh_v0 = ds_v0_tslice.icevolume_nh + var_sh_v0 = ds_v0_tslice.icevolume_sh if varname in ["iceAreaCell", "iceVolumeCell"]: if compare_with_obs: - if compare_with_model: - vars_nh = [var_nh, var_nh_obs, var_nh_model_tocompare] - vars_sh = [var_sh, var_sh_obs, var_sh_model_tocompare] + if ref_casename_v0 != "None": + vars_nh = [var_nh, var_nh_obs, var_nh_v0] + vars_sh = [var_sh, var_sh_obs, var_sh_v0] lineStyles = ['r-', 'k-', 'b-'] lineWidths = [1.2, 1.2, 1.2] else: - # just model and obs + # just v1 model and obs vars_nh = [var_nh, var_nh_obs] vars_sh = [var_sh, var_sh_obs] lineStyles = ['r-', 'k-'] lineWidths = [1.2, 1.2] - elif compare_with_model: - # just - vars_nh = [var_nh, var_nh_model_tocompare] - vars_sh = [var_sh, var_sh_model_tocompare] - lineStyles = ['r-', 'b-'] - lineWidths = [1.2, 1.2] - - if compare_with_obs or compare_with_model: + elif ref_casename_v0 != "None": + # just v1 and v0 models + vars_nh = [var_nh, var_nh_v0] + vars_sh = [var_sh, var_sh_v0] + lineStyles = ['r-', 'b-'] + lineWidths = [1.2, 1.2] + + if compare_with_obs or ref_casename_v0 != "None": # separate plots for nothern and southern hemispheres timeseries_analysis_plot(config, vars_nh, N_movavg, title_nh, xlabel, units, figname_nh, @@ -228,7 +223,6 @@ def seaice_timeseries(config): lineWidths=[1.2, 1.2], title_font_size=title_font_size) - elif varname == "iceThickCell": figname = "%s/%s.%s.png" % (plots_dir,casename,varname) @@ -243,7 +237,6 @@ def seaice_timeseries(config): raise SystemExit("varname variable %s not supported for plotting" % varname) - def replicate_cycle(ds,ds_toreplicate): dsshift = ds_toreplicate.copy() shiftT = ((dsshift.Time.max() - dsshift.Time.min()) @@ -258,5 +251,3 @@ def replicate_cycle(ds,ds_toreplicate): # constrict replicated ds_short to same time dimension as ds_long: dsshift = dsshift.sel(Time=ds.Time.values, method='nearest') return dsshift - - diff --git a/run_analysis.py b/run_analysis.py index 23dac43d3..898d16079 100644 --- a/run_analysis.py +++ b/run_analysis.py @@ -19,22 +19,15 @@ if not os.path.isdir(indir): raise SystemExit("Model directory %s not found. Exiting..." % indir) -any_compare_with_model = False -for section in ['ohc_timeseries', 'sst_timeseries', 'nino34_timeseries', - 'mht_timeseries', 'moc_timeseries', 'seaice_timeseries']: - if config.getboolean(section, 'generate') and config.getboolean(section,'compare_with_model'): - any_compare_with_model = True -if config.getboolean('seaice_modelvsobs', 'generate'): - any_compare_with_model = True - -if any_compare_with_model: +ref_casename_v0 = config.get('case','ref_casename_v0') +if ref_casename_v0 != "None": # we will need model data. Make sure it's there - indir_model_tocompare = config.get('paths','ocndir_model_tocompare') - if not os.path.isdir(indir_model_tocompare): - raise SystemExit("ocndir_model_tocompare directory %s not found. Exiting..." % indir_model_tocompare) - indir_model_tocompare = config.get('paths','seaicedir_model_tocompare') - if not os.path.isdir(indir_model_tocompare): - raise SystemExit("seaicedir_model_tocompare directory %s not found. Exiting..." % indir_model_tocompare) + indir_v0data = config.get('paths','ref_archive_v0_ocndir') + if not os.path.isdir(indir_v0data): + raise SystemExit("ref_archive_v0_ocndir directory %s not found. Exiting..." % indir_v0data) + indir_v0data = config.get('paths','ref_archive_v0_seaicedir') + if not os.path.isdir(indir_v0data): + raise SystemExit("ref_archive_v0_seaicedir directory %s not found. Exiting..." % indir_v0data) if ((config.getboolean('seaice_timeseries', 'generate') and config.getboolean('seaice_timeseries', 'compare_with_obs'))