diff --git a/cases/ctdas-test/config.py b/cases/ctdas-test/config.py new file mode 100644 index 00000000..1504a547 --- /dev/null +++ b/cases/ctdas-test/config.py @@ -0,0 +1,206 @@ +import os +""" +Configuration file for the 'icon-art-oem-vprm-test' case with ICON-ART +""" + +# GENERAL SETTINGS =========================================================== +user = os.environ['USER'] +compute_account = 's1152' +compute_host = 'daint' +compute_queue = 'normal' # 'normal' / 'debug' +constraint = 'mc' # 'mc' / 'gpu' + +Init_from_ICON = False + +target = 'icon-art-oem' +restart_step = 2 * 8760 # hours +restart_cycle_window = 60 * 60 * 24 * 10 #10 days in sec + +if constraint == 'gpu': + ntasks_per_node = 12 +elif constraint == 'mc': + ntasks_per_node = 36 + +#Path to the CTDAS root directory +ctdas_root = '/scratch/snx3000/ekoene/ctdas-icon' +Init_from_ICON = False + +# case name = pathname in cases/ +path = os.path.realpath(__file__) +casename = os.path.basename(os.path.dirname(path)) + +# Root directory of the sourcecode of the chain (where run_chain.py is) +chain_src_dir = os.path.join('/scratch/snx3000/ekoene/processing-chain') + +# Root directory of the working space of the chain +work_root = os.path.join(chain_src_dir, 'work') + +# Directory where executables are stored +# exe_dir = "/scratch/snx3000/ekoene/processing-chain/work/backupVPRM_EU_ERA5_22/2018010100_0_72/icon/run/" # This one works with the current setup... +# exe_dir = "/scratch/snx3000/nponomar/processing_chain_python/processing-chain/work/CTDAS_ZH/2022071400_0_504/icon/run/" # Nikolai's latest file +exe_dir = "/scratch/snx3000/nponomar/icon-vprm-try2/icon-vprm/bin/" +# exe_dir = "/scratch/snx3000/msteiner/ICON_TRANSCOM_SIMULATIONS/executables" + +# Case directory +case_dir = os.path.join(chain_src_dir, 'cases', casename) + +# PREPARE_DATA --------------------------------------------------------------- +input_root_meteo = '/scratch/snx3000/ekoene/ERA5' +meteo_prefix = 'era5_' +meteo_nameformat = meteo_prefix + '%Y%m%d%H' +meteo_suffix = '.nc' +meteo_inc = 3 + +input_root_chem = '/store/empa/em05/dbrunner/icon-art/icbc/' +input_root_icbc = input_root_chem +chem_prefix = 'cams_gqpe_' +# input_root_chem = '/scratch/snx3000/ekoene/MERGE2/' +# input_root_icbc = input_root_chem +# chem_prefix = 'CAMS_' + +chem_nameformat = chem_prefix + '%Y%m%d%H' +chem_suffix = '.nc' +chem_inc = 3 + +icontools_runjobs = [ + # 'icontools_remap_ic_runjob.cfg', + # 'icontools_remap_00_lbc_runjob.cfg', + # 'icontools_remap_lbc_rest_runjob.cfg', + # 'icontools_remap_ic_chem_runjob.cfg', + # 'icontools_remap_lbc_chem_runjob.cfg', +] + +# Icontools executables +#icontools_dir = '/project/s903/mjaehn/spack-install/daint/icontools/master/cce/ldcbgsjjzq2p73xbei7ws4wce5ivzxer/bin/' +icontools_dir = '/scratch/snx3000/nponomar/spack-install/daint/icontools/c2sm-master/gcc/hv7e5pklc6hyntvowrgkywb6rrwzdevb/bin' +#icontools_dir = '/scratch/snx3000/nponomar/icon-vprm/bin' +iconremap_bin = os.path.join(icontools_dir, "iconremap") +iconsub_bin = os.path.join(icontools_dir, "iconsub") + +# Input data for runscript---------------------------------------------------- +# Grid +input_root = '/store/empa/em05/input_iconart_processing_chain_example/' +input_root_grid_icon = '/users/ekoene/CTDAS_inputs/' +radiation_grid_filename = os.path.join(input_root_grid_icon, + "icon_europe_DOM01.parent.nc") +dynamics_grid_filename = os.path.join(input_root_grid_icon, + "icon_europe_DOM01.nc") + +input_root_mapping = '/users/ekoene/CTDAS_inputs' +map_file_ana = os.path.join(input_root_mapping, "map_file.ana") +map_file_latbc = os.path.join(input_root_mapping, "map_file.latbc") +extpar_filename = os.path.join(input_root_grid_icon, + "icon_extpar_EriksGrid.nc") +input_root_rad = os.path.join(input_root, 'rad') +cldopt_filename = os.path.join(input_root_rad, 'rrtm_cldopt.nc') +lrtm_filename = os.path.join(input_root_rad, 'rrtmg_lw.nc') + +# File names ----------------------------------------------------------------- +latbc_filename = "era5__lbc.nc" +inidata_prefix = "era5_init_" +inidata_nameformat = inidata_prefix + '%Y%m%d%H' +inidata_filename_suffix = ".nc" + +output_filename = "icon-art-test" +filename_format = "_DOM_" + +# ART settings---------------------------------------------------------------- +input_root_tracers = '/users/ekoene/CTDAS_inputs/input/xml/' +chemtracer_xml_filename = os.path.join(input_root_tracers, + 'tracers_CoCO2ens_restart.xml') +pntSrc_xml_filename = os.path.join(input_root_tracers, 'pntSrc_example.xml') +# art_input_folder = os.path.join(input_root, 'ART') +# input_root_tracers = '/scratch/snx3000/nponomar/processing_chain_python/processing-chain/work/CTDAS/2022070200_0_672/icon/input/xml/' +# chemtracer_xml_filename = '/scratch/snx3000/nponomar/processing_chain_python/processing-chain/work/VPRM_ENSEMBLE_testing/2022070100_0_24/icon/input/xml/vprm_ensemble_co2_icos_cities_full.xml' + +# pntSrc_xml_filename = os.path.join(input_root_tracers, 'vprm_ensemble_co2_icos_cities_full_restart.xml') +art_input_folder = '/scratch/snx3000/nponomar/ICON_ctdas_msteiner/ART' +# THESE TWO MUST BE AVAILABLE ON SCRATCH! +vprm_regions_synth_nc = '/scratch/snx3000/ekoene/lambdaregions.nc' +vprm_lambdas_synth_nc = '/scratch/snx3000/ekoene/prior_all_ones.nc' + +# OAE ------------------------------------------------------------------------ +# Online anthropogenic emissions +oae_dir = '/users/ekoene/CTDAS_inputs/' +oae_gridded_emissions_nc = 'INV_20180101.nc' +oae_vertical_profiles_nc = 'vertical_profiles_t.nc' +oae_hourofday_nc = 'hourofday.nc' +oae_dayofweek_nc = 'dayofweek.nc' +oae_monthofyear_nc = 'monthofyear.nc' +oae_hourofyear_nc = 'hourofyear8784.nc' + +# VPRM ------------------------------------------------------------------------ +# ICON-ART VPRM coefficients calculated using MODIS data +online_vprm_dir = '/users/ekoene/CTDAS_inputs/' +vprm_coeffs_nc = 'VPRM_indices_ICON.nc' + +# CTDAS ------------------------------------------------------------------------ +ctdas_restart = False +ctdas_BG_run = False +# CTDAS cycle length in days +ctdas_cycle = int(restart_cycle_window / 60 / 60 / 24) +ctdas_nlag = 2 +ctdas_tracer = 'co2' +# CTDAS number of regions and cells->regions file +ctdas_nreg_params = 21184 +ctdas_regionsfile = vprm_regions_synth_nc +# Number of boundaries, and boundaries mask/regions file +ctdas_bg_params = 8 +ctdas_boundary_mask_file = '/scratch/snx3000/ekoene/boundary_mask_bg.nc' +# Number of ensemble members (make this consistent with your XML file!) +ctdas_optimizer_nmembers = 180 +# CTDAS path +ctdas_dir = '/scratch/snx3000/ekoene/ctdas-icon/exec' +# Distance file from region to region (shape: [N_reg x N_reg]), for statevector localization +ctdas_sv_distances = '/scratch/snx3000/ekoene/CTDAS_cells2cells.nc' +# Distance file from region to stations (shape: [N_reg x N_obs]), for observation localization +ctdas_op_loc_coeffs = '/scratch/snx3000/ekoene/cells2stations.nc' +# Directory containing a file with all the observations +ctdas_datadir = '/scratch/snx3000/ekoene/ICOS_extracted/2018/' +# CTDAS localization setting +ctdas_system_localization = 'spatial' +# CTDAS statevector length for one window +ctdas_nparameters = 2 * ctdas_nreg_params + 8 # 2 (A , VPRM) * ctdas_nreg_params + ctdas_bg_params +# Extraction template +ctdas_extract_template = '/scratch/snx3000/ekoene/processing-chain/cases/VPRM_EU_ERA5_22/extract_template_icos_EU' +# ICON runscript template +ctdas_ICON_template = '/scratch/snx3000/ekoene/processing-chain/cases/VPRM_EU_ERA5_22/runscript_template_restart_icos_EU' +# Full path to SBATCH template that can submit extraction scripts +ctdas_sbatch_extract_template = '/scratch/snx3000/ekoene/processing-chain/cases/VPRM_EU_ERA5_22/sbatch_extract_template' +# Full path to possibly time-varying emissionsgrid (if not time-varying, supply a filename without {}!) +ctdas_oae_grid = "/scratch/snx3000/ekoene/inventories/INV_{}.nc" +ctdas_oae_grid_fname = '%Y%m%d' # Specifies the naming scheme to use for the emission grids +# Spinup time length +ctdas_restart_init_time = 60 * 60 * 24 # 1 day in seconds +# Restart file for the first simulation +ctdas_first_restart_init = '/scratch/snx3000/ekoene/processing-chain/work/VPRM_EU_ERA5_22/2018010100_0_240/icon/output_INIT' +# Number of vertical levels +nvlev = 60 +# NOT NEEDED FOR ANYTHING, EXCEPT TO MAKE CTDAS RUN +ctdas_obsoperator_home = '/scratch/snx3000/msteiner/ctdas_test/exec/da/rc/stilt' +ctdas_obsoperator_rc = os.path.join(ctdas_obsoperator_home, 'stilt_0.rc') +ctdas_regtype = 'olson19_oif30' + +# SIMULATION ================================================================= +# ICON ----------------------------------------------------------------------- +# Executable +# icon_bin = os.path.join(exe_dir, "icon.exe") +icon_bin = os.path.join(exe_dir, "icon_vprm_ens_categories_hoy_bg") + +# icon_bin = os.path.join(exe_dir, "20220822_icon_art_offset_nuding_radioactm2") +#icon_bin = os.path.join(exe_dir, "icon_oem_emissions") +# Namelists and slurm runscript templates +icon_runjob = os.path.join(case_dir, 'icon_runjob.cfg') +icon_namelist_master = os.path.join(case_dir, 'icon_master.namelist.cfg') +icon_namelist_nwp = os.path.join(case_dir, 'icon_NAMELIST_NWP.cfg') + +# Walltimes and domain decomposition +if compute_queue == "normal": + icon_walltime = "24:00:00" + icon_np_tot = 32 +elif compute_queue == "debug": + icon_walltime = "00:30:00" + icon_np_tot = 2 +else: + logging.error("Unknown queue name: %s" % compute_queue) + sys.exit(1) diff --git a/cases/ctdas-test/extract_template_icos_EU b/cases/ctdas-test/extract_template_icos_EU new file mode 100644 index 00000000..4c0a69b2 --- /dev/null +++ b/cases/ctdas-test/extract_template_icos_EU @@ -0,0 +1,408 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- + +#%% +"""Import""" +from shapely.geometry import Point, Polygon +from netCDF4 import Dataset +import os +import numpy as np +import datetime +from multiprocessing import Pool +from math import sin, cos, sqrt, atan2, radians +import xarray as xr +from multiprocessing import Process, Manager +from unidecode import unidecode + +#%% + +mountain_stations = [ + 'Jungfraujoch_5', 'Monte Cimone_8', 'Puy de Dome_10', 'Pic du Midi_28', + 'Zugspitze_3', 'Hohenpeissenberg_50', 'Hohenpeissenberg_93', + 'Hohenpeissenberg_131', 'Schauinsland_12', 'Plateau Rosa_10', + 'Laegern-Hochwacht_32' +] +"""Interpolation function""" + + +# def intp_icon_data(args): +def intp_icon_data(iloc, gridinfo, datainfo, latitudes, longitudes, asl, elev, + station_name): + + nn_sel = np.zeros(gridinfo.nn) + u = np.zeros(gridinfo.nn) + + R = 6373.0 # approximate radius of earth in km + + if (radians(longitudes[iloc]) < np.nanmin(gridinfo.clon)) or (radians( + longitudes[iloc]) > np.nanmax(gridinfo.clon)): + u[:] = np.nan + return np.zeros((gridinfo.nn)), np.zeros( + (gridinfo.nn)).astype(int), np.zeros( + (gridinfo.nn)).astype(int), nn_sel[:], u[:] + + if (radians(latitudes[iloc]) < np.nanmin(gridinfo.clat)) or (radians( + latitudes[iloc]) > np.nanmax(gridinfo.clat)): + u[:] = np.nan + return np.zeros((gridinfo.nn)), np.zeros( + (gridinfo.nn)).astype(int), np.zeros( + (gridinfo.nn)).astype(int), nn_sel[:], u[:] + + #% + lat1 = radians(latitudes[iloc]) + lon1 = radians(longitudes[iloc]) + + #% + """FIND 4 CLOSEST CENTERS""" + distances = np.zeros((len(gridinfo.clon))) + for icell in np.arange(len(gridinfo.clon)): + lat2 = gridinfo.clat[icell] + lon2 = gridinfo.clon[icell] + dlon = lon2 - lon1 + dlat = lat2 - lat1 + a = sin(dlat / 2)**2 + cos(lat1) * cos(lat2) * sin(dlon / 2)**2 + c = 2 * atan2(sqrt(a), sqrt(1 - a)) + distances[icell] = R * c + nn_sel[:] = [ + x for _, x in sorted(zip(distances, np.arange(len(gridinfo.clon)))) + ][0:gridinfo.nn] + nn_sel = nn_sel.astype(int) + #print('---nn_sel:',nn_sel) + #print('---distances[0:gridinfo.nn]:',distances[0:gridinfo.nn]) + u[:] = [1. / distances[y] for y in nn_sel] + print('---distances:', [distances[y] for y in nn_sel]) + print('---weights:', u) + + #% + """Calculate vertical interpolation factor""" + idx_above = -1 * np.ones((len(nn_sel))).astype(int) + idx_below = -1 * np.ones((len(nn_sel))).astype(int) + target_asl = np.zeros((len(nn_sel))) + for nnidx in np.arange(len(nn_sel)): + model_topo = datainfo.z_ifc[-1, nn_sel[nnidx]] + print(station_name[iloc]) + if station_name[iloc] not in mountain_stations: + print('not a mountain station') + target_asl[nnidx] = model_topo + elev[iloc] + else: + print('mountain station') + target_asl[nnidx] = asl[iloc] + elev[iloc] + for i_mc, mc in enumerate(datainfo.z_mc[:, nn_sel[nnidx]]): + # if mc>=asl[iloc]: + if mc >= target_asl[nnidx]: + idx_above[nnidx] = i_mc + else: + idx_below[nnidx] = i_mc + break + + #in case of data point below lowest midlevel: + for nnidx in np.arange(len(nn_sel)): + if idx_below[nnidx] == -1: + idx_below[nnidx] = idx_above[nnidx] + if any(np.ravel([idx_above, idx_below]) < 0): + print("At least one nearest neigbor has no valid height levels") + u[:] = np.nan #-999. + return np.zeros((gridinfo.nn)), np.zeros( + (gridinfo.nn)).astype(int), np.zeros( + (gridinfo.nn)).astype(int), nn_sel[:], u + + vert_scaling_fact = np.zeros((len(nn_sel))) + for nnidx in np.arange(len(nn_sel)): + if idx_below[nnidx] != idx_above[nnidx]: + vert_scaling_fact[nnidx] = ( + target_asl[nnidx] - + datainfo.z_mc[idx_below[nnidx], nn_sel[nnidx]]) / ( + datainfo.z_mc[idx_above[nnidx], nn_sel[nnidx]] - + datainfo.z_mc[idx_below[nnidx], nn_sel[nnidx]]) + else: + vert_scaling_fact[nnidx] = 0. + + print('---idx above:', idx_above) + print('---idx below:', idx_below) + print('---vert_scaling_fact:', vert_scaling_fact) + + #% + + return vert_scaling_fact, idx_below, idx_above, nn_sel[:], u + + +#%% +"""Paths, Variables and Constants""" + +ICON_grid = '{cfg.dynamics_grid_filename}' +fname_base = '{fname_base}' +DATA_path = '{icon_path}' +starttime = '{start}' +enddtime = '{end}' +obs_dir = '{stationdir}' +nlev = {cfg.nvlev} #number of vertical levels +nneighb = {nneighb} #number of nearest neighbors to consider + +#%% +"""Variables to extract""" + +n_member = {cfg.ctdas_optimizer_nmembers} + +meta = {{ + 'TRCO2_A': {{ + 'offset': 0. + }}, + 'TRCO2_BG': {{ + 'offset': 0. + }}, + 'CO2_RA': {{ + 'offset': 0. + }}, + 'CO2_GPP': {{ + 'offset': 0. + }}, + 'u': {{ + 'offset': 0. + }}, + 'v': {{ + 'offset': 0. + }}, + 'qv': {{ + 'offset': 0. + }}, + 'temp': {{ + 'offset': 0. + }}, + 'biosink_chemtr': {{ + 'offset': 0. + }}, + 'biosource_all_chemtr': {{ + 'offset': 0. + }}, + 'TRCO2_A-ENS': {{ + 'offset': 0, + 'ensemble': n_member + }}, +}} + +#%% +"""Get grid data""" + +fh_grid = Dataset(ICON_grid, 'r') + + +class gridinfo: + clon_vertices = np.array(fh_grid.variables['clon_vertices']) + clat_vertices = np.array(fh_grid.variables['clat_vertices']) + cells_of_vertex = np.array(fh_grid.variables['cells_of_vertex']) + vertex_of_cell = np.array(fh_grid.variables['vertex_of_cell']) + neighbor_cell_index = np.array(fh_grid.variables['neighbor_cell_index']) + vlon = np.array(fh_grid.variables['vlon']) + vlat = np.array(fh_grid.variables['vlat']) + clon = np.array(fh_grid.variables['clon']) + clat = np.array(fh_grid.variables['clat']) + ncells = len(fh_grid.dimensions['cell']) + nn = nneighb + + +#%% +"""Times""" + +firstfile = True +startdate = datetime.datetime.strptime(starttime, '%Y-%m-%d %H:%M:%S') +enddate = datetime.datetime.strptime(enddtime, '%Y-%m-%d %H:%M:%S') +delta = datetime.timedelta(hours=1) +looptime = startdate + +#%% +"""Get locations of measurement stations""" + +longitudes = [] +latitudes = [] +obsnames = [] +# stationnames = [] +asl = [] +elev = [] +#for station in stationlist: +for ncfile in os.listdir(obs_dir): + if not ncfile.endswith('.nc'): continue + + infile = os.path.join(obs_dir, ncfile) + + f = xr.open_dataset(infile) + stationnames = f.Stations_names.values + st_ind = np.arange(len(stationnames)) + for x in st_ind: + latitudes.append(f.Lat[x].values) + obsnames.append(unidecode(str(f.Stations_names[x].values))) + longitudes.append(f.Lon[x].values) + asl.append(f.Stations_masl[x]) + elev.append(f.Stations_elev[x]) +print("Found %i locations." % (len(latitudes))) + +# """Add 5 missing stations""" +# missing_longitudes = [-1.15, 4.93, 0.23, 8.4, 8.18] +# missing_latitudes = [54.36, 51.97, 50.98, 47.48, 47.19] +# missing_stationnames = ['bsd', 'cbw', 'hea', 'lae', 'beo'] +# missing_asl = [628.,200.,250.,872.,1009.] +# for imiss in np.arange(len(missing_longitudes)): +# latitudes.append(missing_latitudes[imiss]) +# longitudes.append(missing_longitudes[imiss]) +# stationnames.append(missing_stationnames[imiss]) +# asl.append(missing_asl[imiss]) +# print("Added %i missing locations."%(len(missing_longitudes))) + +#%% +"""Initialize output variables""" +n_det = int( + np.nansum([1 for var in meta.keys() if 'ensemble' not in meta[var]])) +n_ens = int(np.nansum([1 for var in meta.keys() if 'ensemble' in meta[var]])) +intp_ICON_data_det = np.zeros((n_det, len(latitudes), 0)) +maxmem = 0 +for var in meta.keys(): + if 'ensemble' in meta[var]: + if meta[var]['ensemble'] > maxmem: + maxmem = meta[var]['ensemble'] +maxmem = int(maxmem) +intp_ICON_data_ens = np.zeros((n_ens, maxmem, len(latitudes), 0)) +#%% +"""Loop over Data Files (=timesteps)""" + + +def process_data(index, gridinfo, datainfo, latitudes, longitudes, asl, elev, + obsnames, results): + result = intp_icon_data(index, gridinfo, datainfo, latitudes, longitudes, + asl, elev, obsnames) + results.append((index, result)) + + +datetime_list = [] +print('======================================') +date_idx = 0 +while looptime <= enddate: + + intp_ICON_data_det = np.concatenate( + (intp_ICON_data_det, np.zeros((n_det, len(latitudes), 1))), axis=2) + intp_ICON_data_ens = np.concatenate( + (intp_ICON_data_ens, np.zeros((n_ens, maxmem, len(latitudes), 1))), + axis=3) + + timestring = datetime.datetime.strftime(looptime, '%Y-%m-%dT%H') + datetime_list.append(timestring) + DATA_file = os.path.join(DATA_path, + '%s_%s:00:00.000.nc' % (fname_base, timestring)) + + print('extracting from %s' % (DATA_file), flush=True) + + fh_data = Dataset(DATA_file, 'r') + + class datainfo: + z_mc = np.array(fh_data.variables['z_mc']) + z_ifc = np.array(fh_data.variables['z_ifc']) + + ICON_data_det = np.zeros((n_det, nlev, gridinfo.ncells)) + ICON_data_ens = np.zeros((n_ens, maxmem, nlev, gridinfo.ncells)) + ivar = 0 + for var in meta.keys(): + if not 'ensemble' in meta[var]: + ICON_data_det[ivar, ...] = np.array( + fh_data.variables[var]) - meta[var]['offset'] + ivar += 1 + ivar = 0 + for var in meta.keys(): + if 'ensemble' in meta[var]: + # for iens in np.arange(n_member): + for iens in np.arange(meta[var]['ensemble']): + varnc = var.split('-')[0] + '-%.3i' % (iens + 1) + ICON_data_ens[ivar, iens, ...] = np.array( + fh_data.variables[varnc]) - meta[var]['offset'] + ivar += 1 +#%% + """Since the stations don't walk around, I only call the function at the first timestep""" + + if looptime == startdate: + manager = Manager() + results = manager.list() + processes = [] + + for i in range(len(latitudes)): + p = Process(target=process_data, + args=(i, gridinfo, datainfo, latitudes, longitudes, + asl, elev, obsnames, results)) + processes.append(p) + p.start() + + for p in processes: + p.join() + + # Sort the results based on the index + results = sorted(results, key=lambda x: x[0]) + + # Extract the sorted results + sorted_results = [result for _, result in results] + vsf, idxb, idxa, neighbours, u_ret = zip(*sorted_results) + + vsf = np.array(vsf) + idxb = np.array(idxb, dtype=int) + idxa = np.array(idxa, dtype=int) + neighbours = np.array(neighbours, dtype=int) + u_ret = np.array(u_ret) + + #Do the interpolation + for iloc in np.arange(len(latitudes)): + + ###First, the deterministic values: + ##First, the vertical interpolation: + vert_intp_data = np.zeros((n_det, len(idxb[iloc]))) + for nn in np.arange(len(idxb[iloc])): + vert_intp_data[:,nn] = ICON_data_det[:,idxb[iloc,nn],neighbours[iloc,nn]] + vsf[iloc,nn]*(ICON_data_det[:,idxa[iloc,nn],neighbours[iloc,nn]] \ + -ICON_data_det[:,idxb[iloc,nn],neighbours[iloc,nn]]) + ##Now the horizontal interpolation: + #intp_ICON_data_det[:,iloc,date_idx] = u_ret[iloc,0]*vert_intp_data[:,0]+u_ret[iloc,1]*vert_intp_data[:,1] \ + # +u_ret[iloc,2]*vert_intp_data[:,2] + # + intp_ICON_data_det[:, iloc, date_idx] = np.nansum( + [w * vert_intp_data[:, i] for i, w in enumerate(u_ret[iloc, :])], + axis=0) / np.nansum(u_ret[iloc, :]) + ###Second, the ensemble values: + vert_intp_data = np.zeros((n_ens, maxmem, len(idxb[iloc]))) + for nn in np.arange(len(idxb[iloc])): + vert_intp_data[:,:,nn] = ICON_data_ens[:,:,idxb[iloc,nn],neighbours[iloc,nn]] + vsf[iloc,nn]*(ICON_data_ens[:,:,idxa[iloc,nn],neighbours[iloc,nn]] \ + -ICON_data_ens[:,:,idxb[iloc,nn],neighbours[iloc,nn]]) + #intp_ICON_data_ens[:,:,iloc,date_idx] = u_ret[iloc,0]*vert_intp_data[:,:,0]+u_ret[iloc,1]*vert_intp_data[:,:,1] \ + # +u_ret[iloc,2]*vert_intp_data[:,:,2] + intp_ICON_data_ens[:, :, iloc, date_idx] = np.nansum( + [ + w * vert_intp_data[:, :, i] + for i, w in enumerate(u_ret[iloc, :]) + ], + axis=0) / np.nansum(u_ret[iloc, :]) +#%% + """Update time""" + looptime += delta + date_idx += 1 +#%% +"""Save as netcdf""" +with Dataset("{extracted_file}", mode='w') as ofile: + + osites = ofile.createDimension('sites', len(latitudes)) + otime = ofile.createDimension('time', (date_idx)) + + oname = ofile.createVariable('site_name', str, ('sites')) + otimes = ofile.createVariable('time', np.unicode_, ('time')) + + ivar = 0 + for var in meta.keys(): + if 'ensemble' not in meta[var]: + ovar = ofile.createVariable(var, np.float32, ('sites', 'time')) + ovar[:, :] = intp_ICON_data_det[ivar, :, :] + ivar += 1 + ivar = 0 + for var in meta.keys(): + if 'ensemble' in meta[var]: + oens = ofile.createDimension('ens_%.2i' % (ivar + 1), + meta[var]['ensemble']) + varnc = var.split('-')[0] + '_ENS' + ovar = ofile.createVariable( + varnc, np.float32, ('ens_%.2i' % (ivar + 1), 'sites', 'time')) + ovar[:, :, :] = intp_ICON_data_ens[ivar, + 0:meta[var]['ensemble'], :, :] + ivar += 1 + + oname[:] = np.array(stationnames[:]) + otimes[:] = np.array(datetime_list) diff --git a/cases/ctdas-test/runscript_template_restart_icos_EU b/cases/ctdas-test/runscript_template_restart_icos_EU new file mode 100644 index 00000000..d22542c3 --- /dev/null +++ b/cases/ctdas-test/runscript_template_restart_icos_EU @@ -0,0 +1,423 @@ +#!/usr/bin/env bash +#SBATCH --job-name="{cfg.casename}_{cfg.inidate_yyyymmddhh}_{cfg.forecasttime}" +#SBATCH --account={cfg.compute_account} +#SBATCH --time="00:50:00" +#SBATCH --nodes={cfg.icon_np_tot} +#SBATCH --ntasks-per-core=1 +#SBATCH --ntasks-per-node={cfg.ntasks_per_node} +#SBATCH --cpus-per-task=1 +#SBATCH --partition={cfg.compute_queue} +#SBATCH --constraint={cfg.constraint} +#SBATCH --open-mode=append +#SBATCH --chdir={cfg.icon_work} + +set -x +export OMP_NUM_THREADS=$SLURM_CPUS_PER_TASK + +#export ECCODES_DEFINITION_PATH=/users/msteiner/eccodes_definitions/definitions.edzw-2.12.5-2:/store/empa/em05/easybuild/software/ecCodes/2.12.5-CrayGNU-19.10/share/eccodes/definitions +export ECCODES_DEFINITION_PATH=/store/empa/em05/eccodes_definitions/definitions.edzw-2.12.5-2 +export ECCODES_DEFINITION_PATH=$ECCODES_DEFINITION_PATH:/store/empa/em05/easybuild.backup/software/ecCodes/2.12.5-CrayGNU-20.08/share/eccodes/definitions +export ECCODES_DEFINITION_PATH=$ECCODES_DEFINITION_PATH:/project/g110/spack-install/daint/eccodes/2.19.0/pgi/6skdmw5lsn6mjv4esxkyalf6xogllshi/share/eccodes/definitions/ +echo $ECCODES_DEFINITION_PATH + + +# ---------------------------------------------------------------------------- +# Link radiation input files +# ---------------------------------------------------------------------------- +ln -sf {cfg.art_input_folder}/runctrl_examples/photo_ctrl/* . +ln -sf {cfg.art_input_folder}/runctrl_examples/init_ctrl/* . + +# ---------------------------------------------------------------------------- +# create ICON master namelist +# ---------------------------------------------------------------------------- + +cat > icon_master.namelist << EOF +! master_nml: ---------------------------------------------------------------- +&master_nml + lrestart = .FALSE. ! .TRUE.=current experiment is resumed +/ + +! master_model_nml: repeated for each model ---------------------------------- +&master_model_nml + model_type = 1 ! identifies which component to run (atmosphere,ocean,...) + model_name = "ATMO" ! character string for naming this component. + model_namelist_filename = "NAMELIST_NWP" ! file name containing the model namelists + model_min_rank = 1 ! start MPI rank for this model + model_max_rank = 65536 ! end MPI rank for this model + model_inc_rank = 1 ! stride of MPI ranks +/ + +! time_nml: specification of date and time------------------------------------ +&time_nml + ini_datetime_string = "{ini_restart_string}" ! initial date and time of the simulation + end_datetime_string = "{ini_restart_end_string}" ! end date and time of the simulation 10T00 +/ +EOF + +# ---------------------------------------------------------------------- +# model namelists +# ---------------------------------------------------------------------- + +cat > NAMELIST_NWP << EOF +! parallel_nml: MPI parallelization ------------------------------------------- +¶llel_nml + nproma = 16 ! loop chunk length + p_test_run = .FALSE. ! .TRUE. means verification run for MPI parallelization + num_io_procs = 3 ! number of I/O processors + num_restart_procs = 0 ! number of restart processors + num_prefetch_proc = 1 ! number of processors for LBC prefetching + iorder_sendrecv = 3 ! sequence of MPI send/receive calls +/ + + +! run_nml: general switches --------------------------------------------------- +&run_nml + ltestcase = .FALSE. ! real case run + num_lev = 60 ! number of full levels (atm.) for each domain + lvert_nest = .FALSE. ! no vertical nesting + dtime = 120 ! timestep in seconds + ldynamics = .TRUE. ! compute adiabatic dynamic tendencies + ltransport = .TRUE. ! compute large-scale tracer transport + ntracer = 0 ! number of advected tracers + iforcing = 3 ! forcing of dynamics and transport by parameterized processes + msg_level = 7 ! detailed report during integration + ltimer = .TRUE. ! timer for monitoring the runtime of specific routines + timers_level = 10 ! performance timer granularity + check_uuid_gracefully = .TRUE. ! give only warnings for non-matching uuids + output = "nml" ! main switch for enabling/disabling components of the model output + lart = .TRUE. ! main switch for ART + debug_check_level = 10 +/ + +! art_nml: Aerosols and Reactive Trace gases extension------------------------------------------------- +&art_nml + lart_chem = .TRUE. ! enables chemistry + lart_pntSrc = .FALSE. ! enables point sources + lart_aerosol = .FALSE. ! main switch for the treatment of atmospheric aerosol + lart_chemtracer = .TRUE. ! main switch for the treatment of chemical tracer + lart_diag_out = .TRUE. ! If this switch is set to .TRUE., diagnostic + ! ... output elds are available. Set it to + ! ... .FALSE. when facing memory problems. + cart_chemtracer_xml = '{cfg.chemtracer_xml_filename_scratch}' ! path to xml file for passive tracers + cart_input_folder = '{cfg.art_input_folder}' ! absolute Path to ART source code +/ + +! oem_nml: online emission module --------------------------------------------- +&oemctrl_nml + gridded_emissions_nc = '{emissionsgrid}' + vertical_profile_nc = '{cfg.oae_vertical_profiles_nc_scratch}' + hour_of_day_nc = '{cfg.oae_hourofday_nc_scratch}' + day_of_week_nc = '{cfg.oae_dayofweek_nc_scratch}' + month_of_year_nc = '{cfg.oae_monthofyear_nc_scratch}' + hour_of_year_nc = '{cfg.oae_hourofyear_nc_scratch}' + chem_init_nc = '{inifile}' + ens_reg_nc = '{cfg.vprm_regions_synth_nc_scratch}' + ens_lambda_nc = '{lambda_f}' + chem_restart_nc = '{restart_file}' + restart_init_time = {cfg.ctdas_restart_init_time} ![s] + vegetation_indices_nc = '{cfg.vprm_coeffs_nc_scratch}' + boundary_lambda_nc = '{bg_lambda}' + boundary_regions_nc = '{cfg.ctdas_boundary_mask_file}' + vprm_par = 3.139508E+02, 3.132859E+02, 5.149856E+02, 1.009878E+02, 6.820000E+02, 1.322951E+03, 5.794354E+02, 0.0E+00 + vprm_lambda = -1.968102E-01, -1.808074E-01, -1.428753E-01, -2.019100E-01, -1.141000E-01, -8.061275E-02, -1.704107E-01, 0.0E+00 + vprm_alpha = 2.237246E-01, 1.273649E-01, 1.714644E-01, 5.186480E-02, 4.900000E-03, 6.764283E-02, 8.623752E-02, 0.0E+00 + vprm_beta = -6.411373E-01, 1.140385E+00, 1.072367E-02, -1.675250E-01, 0.000000E+00, 5.772793E-01, 3.629423E-01, 0.0E+00 + vprm_tmin = 0.0, 0.0, 0.0, 2.0, 2.0, 5.0, 2.0, 0.0 + vprm_tmax = 40.0, 40.0, 40.0, 40.0, 40.0, 40.0, 40.0, 0.0 + vprm_topt = 20.0, 20.0, 20.0, 20.0, 20.0, 22.0, 18.0, 0.0 + vprm_tlow = 4.0, 0.0, 2.0, 4.0, 0.0, 0.0, 0.0, 0.0 + lcut_area = .FALSE. + lon_cut_start = 5.0 + lon_cut_end = 15.0 + lat_cut_start = 45.0 + lat_cut_end = 55.0 +/ + + +! diffusion_nml: horizontal (numerical) diffusion ---------------------------- +&diffusion_nml + lhdiff_vn = .TRUE. ! diffusion on the horizontal wind field + lhdiff_temp = .TRUE. ! diffusion on the temperature field + lhdiff_w = .TRUE. ! diffusion on the vertical wind field + hdiff_order = 5 ! order of nabla operator for diffusion + itype_vn_diffu = 1 ! reconstruction method used for Smagorinsky diffusion + itype_t_diffu = 2 ! discretization of temperature diffusion + hdiff_efdt_ratio = 24.0 ! ratio of e-folding time to time step + hdiff_smag_fac = 0.025 ! scaling factor for Smagorinsky diffusion +/ + +! dynamics_nml: dynamical core ----------------------------------------------- +&dynamics_nml + iequations = 3 ! type of equations and prognostic variables + idiv_method = 1 ! method for divergence computation + divavg_cntrwgt = 0.50 ! weight of central cell for divergence averaging + lcoriolis = .TRUE. ! Coriolis force +/ + +! extpar_nml: external data -------------------------------------------------- +&extpar_nml + itopo = 1 ! topography (0:analytical) + extpar_filename = '{cfg.extpar_filename_scratch}' ! filename of external parameter input file + n_iter_smooth_topo = 1,1 ! iterations of topography smoother + heightdiff_threshold = 3000. ! height difference between neighb. grid points + hgtdiff_max_smooth_topo = 750.,750. ! see Namelist doc + heightdiff_threshold = 2250.,1500. +/ + +! initicon_nml: specify read-in of initial state ------------------------------ +&initicon_nml + init_mode = 2 ! 7: start from DWD fg with subsequent vertical remapping + lread_ana = .FALSE. ! no analysis data will be read + ifs2icon_filename = "{inifile}" ! initial data filename + ana_varnames_map_file = "{cfg.map_file_ana_scratch}" ! dictionary mapping internal names onto GRIB2 shortNames + ltile_coldstart = .TRUE. ! coldstart for surface tiles + ltile_init = .FALSE. ! set it to .TRUE. if FG data originate from run without tiles + +/ + +! grid_nml: horizontal grid -------------------------------------------------- +&grid_nml + dynamics_grid_filename = "{cfg.dynamics_grid_filename_scratch}" ! array of the grid filenames for the dycore + radiation_grid_filename = "{cfg.radiation_grid_filename_scratch}" ! array of the grid filenames for the radiation model + dynamics_parent_grid_id = 0 ! array of the indexes of the parent grid filenames + lredgrid_phys = .TRUE. ! .true.=radiation is calculated on a reduced grid + lfeedback = .TRUE. ! specifies if feedback to parent grid is performed + l_limited_area = .TRUE. ! .TRUE. performs limited area run + ifeedback_type = 2 ! feedback type (incremental/relaxation-based) + start_time = 0. ! Time when a nested domain starts to be active [s] +/ + +! gridref_nml: grid refinement settings -------------------------------------- +&gridref_nml + denom_diffu_v = 150. ! denominator for lateral boundary diffusion of velocity +/ + +! interpol_nml: settings for internal interpolation methods ------------------ +&interpol_nml + nudge_zone_width = 10 ! width of lateral boundary nudging zone + support_baryctr_intp = .FALSE. ! barycentric interpolation support for output + nudge_max_coeff = 0.069 + nudge_efold_width = 2.0 +/ + +! io_nml: general switches for model I/O ------------------------------------- +&io_nml + itype_pres_msl = 5 ! method for computation of mean sea level pressure + itype_rh = 1 ! method for computation of relative humidity + lmask_boundary = .TRUE. ! mask out interpolation zone in output +/ + +! limarea_nml: settings for limited area mode --------------------------------- +&limarea_nml + itype_latbc = 1 ! 1: time-dependent lateral boundary conditions + dtime_latbc = 10800 ! time difference between 2 consecutive boundary data + nlev_latbc = 137 ! Number of vertical levels in boundary data + latbc_boundary_grid = "{cfg.lateral_boundary_grid_scratch}" ! Grid file defining the lateral boundary + latbc_path = "{cfg.icon_input_icbc}" ! Absolute path to boundary data + latbc_varnames_map_file = "{cfg.map_file_latbc_scratch}" + latbc_filename = "{cfg.latbc_filename}" ! boundary data input filename + init_latbc_from_fg = .FALSE. ! .TRUE.: take lbc for initial time from first guess +/ + +! lnd_nml: land scheme switches ----------------------------------------------- +&lnd_nml + ntiles = 3 ! number of tiles + nlev_snow = 3 ! number of snow layers + lmulti_snow = .FALSE. ! .TRUE. for use of multi-layer snow model + idiag_snowfrac = 20 ! type of snow-fraction diagnosis + lsnowtile = .TRUE. ! .TRUE.=consider snow-covered and snow-free separately + itype_root = 2 ! root density distribution + itype_heatcond = 3 ! type of soil heat conductivity + itype_lndtbl = 4 ! table for associating surface parameters + itype_evsl = 4 ! type of bare soil evaporation + itype_root = 2 ! root density distribution + cwimax_ml = 5.e-4 ! scaling parameter for max. interception storage + c_soil = 1.75 ! surface area density of the evaporative soil surface + c_soil_urb = 0.5 ! same for urban areas + lseaice = .TRUE. ! .TRUE. for use of sea-ice model + llake = .TRUE. ! .TRUE. for use of lake model +/ + +! nonhydrostatic_nml: nonhydrostatic model ----------------------------------- +&nonhydrostatic_nml + iadv_rhotheta = 2 ! advection method for rho and rhotheta + ivctype = 2 ! type of vertical coordinate + itime_scheme = 4 ! time integration scheme + ndyn_substeps = 5 ! number of dynamics steps per fast-physics step + exner_expol = 0.333 ! temporal extrapolation of Exner function + vwind_offctr = 0.2 ! off-centering in vertical wind solver + damp_height = 12500.0 ! height at which Rayleigh damping of vertical wind starts + rayleigh_coeff = 1.5 ! Rayleigh damping coefficient + divdamp_order = 24 ! order of divergence damping + divdamp_type = 3 ! type of divergence damping + divdamp_fac = 0.004 ! scaling factor for divergence damping + l_open_ubc = .FALSE. ! .TRUE.=use open upper boundary condition + igradp_method = 3 ! discretization of horizontal pressure gradient + l_zdiffu_t = .TRUE. ! specifies computation of Smagorinsky temperature diffusion + thslp_zdiffu = 0.02 ! slope threshold (temperature diffusion) + thhgtd_zdiffu = 125.0 ! threshold of height difference (temperature diffusion) + htop_moist_proc = 22500.0 ! max. height for moist physics + hbot_qvsubstep = 22500.0 ! height above which QV is advected with substepping scheme +/ + +! nwp_phy_nml: switches for the physics schemes ------------------------------ +&nwp_phy_nml + inwp_gscp = 2 ! cloud microphysics and precipitation + inwp_convection = 1 ! convection + lshallowconv_only = .FALSE. ! only shallow convection + inwp_radiation = 1 ! radiation + inwp_cldcover = 1 ! cloud cover scheme for radiation + inwp_turb = 1 ! vertical diffusion and transfer + inwp_satad = 1 ! saturation adjustment + inwp_sso = 1 ! subgrid scale orographic drag + inwp_gwd = 0 ! non-orographic gravity wave drag + inwp_surface = 1 ! surface scheme + latm_above_top = .TRUE. ! take into account atmosphere above model top for radiation computation + ldetrain_conv_prec = .TRUE. + efdt_min_raylfric = 7200. ! minimum e-folding time of Rayleigh friction + itype_z0 = 2 ! type of roughness length data + icapdcycl = 3 ! apply CAPE modification to improve diurnalcycle over tropical land + icpl_aero_conv = 1 ! coupling between autoconversion and Tegen aerosol climatology + icpl_aero_gscp = 1 ! coupling between autoconversion and Tegen aerosol climatology + lrtm_filename = '{cfg.lrtm_filename_scratch}' ! longwave absorption coefficients for RRTM_LW + cldopt_filename = '{cfg.cldopt_filename_scratch}' ! RRTM cloud optical properties + dt_rad = 720. ! time step for radiation in s + dt_conv = 240. ! time step for convection in s (domain specific) + dt_sso = 240. ! time step for SSO parameterization + dt_gwd = 360. ! time step for gravity wave drag parameterization +/ + +! nwp_tuning_nml: additional tuning parameters ---------------------------------- +&nwp_tuning_nml + itune_albedo = 1 ! reduced albedo (w.r.t. MODIS data) over Sahara + tune_gkwake = 1.8 + tune_gkdrag = 0.01 + tune_minsnowfrac = 0.3 +/ + +! output_nml: specifies an output stream -------------------------------------- +&output_nml + filetype = 4 ! output format: 2=GRIB2, 4=NETCDFv2 + dom = 1 ! write domain 1 only + output_bounds = 0., 10000000., 43200. ! start, end, increment + steps_per_file = 1 ! number of steps per file + mode = 1 ! 1: forecast mode (relative t-axis), 2: climate mode (absolute t-axis) + include_last = .TRUE. + output_filename = 'ICON-ART-OEM' + filename_format = '{output_directory}/_' ! file name base + steps_per_file_inclfirst = .FALSE. + output_grid = .TRUE. + remap = 1 ! 1: remap to lat-lon grid + !north_pole = -170.,40. ! definition of north_pole for rotated lat-lon grid + reg_lon_def = -8.25,0.05,17.65 ! + reg_lat_def = 40.75,0.05,58.85 ! + ml_varlist = 'group:ART_CHEMISTRY', + 'group:ART_PASSIVE', + 'group:ATMO_ML_VARS', + 'z_mc', 'z_ifc', + +/ + +! output_nml: specifies an output stream -------------------------------------- +&output_nml + filetype = 4 ! output format: 2=GRIB2, 4=NETCDFv2 + dom = 1 ! write domain 1 only + output_bounds = 0., 10000000., 3600. ! start, end, increment + steps_per_file = 1 ! number of steps per file + mode = 1 ! 1: forecast mode (relative t-axis), 2: climate mode (absolute t-axis) + include_last = .TRUE. + output_filename = 'ICON-ART-UNSTR' + filename_format = '{output_directory}/_' ! file name base + steps_per_file_inclfirst = .FALSE. + output_grid = .TRUE. + remap = 0 ! 1: remap to lat-lon grid + !north_pole = -170.,40. ! definition of north_pole for rotated lat-lon grid + reg_lon_def = -8.25,0.05,17.65 ! + reg_lat_def = 40.75,0.05,58.85 ! + ml_varlist = 'group:PBL_VARS', + 'group:ATMO_ML_VARS', + 'group:precip_vars', + 'group:land_vars', + 'group:nh_prog_vars', + 'group:ART_CHEMISTRY', + 'z_mc', 'z_ifc', + 'topography_c', + 'group:ART_PASSIVE', + 'group:ART_AEROSOL' +/ + +! output_nml: specifies an output stream -------------------------------------- +&output_nml + filetype = 4 ! output format: 2=GRIB2, 4=NETCDFv2 + dom = 1 ! write domain 1 only + output_bounds = 604800., 10000000., 604800. ! start, end, increment 518400 + steps_per_file = 1 ! number of steps per file + mode = 1 ! 1: forecast mode (relative t-axis), 2: climate mode (absolute t-axis) + include_last = .TRUE. + output_filename = 'ICON-ART-OEM-INIT' + filename_format = '{output_directory}/_' ! file name base + steps_per_file_inclfirst = .FALSE. + output_grid = .TRUE. + remap = 0 ! 1: remap to lat-lon grid + ml_varlist = 'group:ART_CHEMISTRY' +/ + +! radiation_nml: radiation scheme --------------------------------------------- +&radiation_nml + irad_o3 = 7 ! ozone climatology + irad_aero = 6 ! aerosols + albedo_type = 2 ! type of surface albedo + vmr_co2 = 390.e-06 + vmr_ch4 = 1800.e-09 + vmr_n2o = 322.0e-09 + vmr_o2 = 0.20946 + vmr_cfc11 = 240.e-12 + vmr_cfc12 = 532.e-12 + ecrad_data_path = '/scratch/snx3000/nponomar/icon-vprm/externals/ecrad/data' +/ + +! sleve_nml: vertical level specification ------------------------------------- +&sleve_nml + min_lay_thckn = 20.0 ! layer thickness of lowermost layer + top_height = 23000.0 ! height of model top + stretch_fac = 0.65 ! stretching factor to vary distribution of model levels + decay_scale_1 = 4000.0 ! decay scale of large-scale topography component + decay_scale_2 = 2500.0 ! decay scale of small-scale topography component + decay_exp = 1.2 ! exponent of decay function + flat_height = 16000.0 ! height above which the coordinate surfaces are flat +/ + +! transport_nml: tracer transport --------------------------------------------- +&transport_nml + npassive_tracer = 0 ! number of additional passive tracers + ivadv_tracer = 3, 3, 3, 3, 3, 3 ! tracer specific method to compute vertical advection + itype_hlimit = 3, 4, 4, 4, 4, 4 ! type of limiter for horizontal transport + ihadv_tracer = 52, 2, 2, 2, 2, 2 ! tracer specific method to compute horizontal advection + llsq_svd = .TRUE. ! use SV decomposition for least squares design matrix +/ + +! turbdiff_nml: turbulent diffusion ------------------------------------------- +&turbdiff_nml + tkhmin = 0.75 ! scaling factor for minimum vertical diffusion coefficient + tkmmin = 0.75 ! scaling factor for minimum vertical diffusion coefficient + pat_len = 750.0 ! effective length scale of thermal surface patterns + c_diff = 0.2 ! length scale factor for vertical diffusion of TKE + rat_sea = 0.8 ! controls laminar resistance for sea surface + ltkesso = .TRUE. ! consider TKE-production by sub-grid SSO wakes + frcsmot = 0.2 ! these 2 switches together apply vertical smoothing of the TKE source terms + imode_frcsmot = 2 ! in the tropics (only), which reduces the moist bias in the tropical lower troposphere + itype_sher = 3 ! type of shear forcing used in turbulence + ltkeshs = .TRUE. ! include correction term for coarse grids in hor. shear production term + a_hshr = 2.0 ! length scale factor for separated horizontal shear mode + icldm_turb = 1 ! mode of cloud water representation in turbulence + ldiff_qi = .TRUE. +/ + +EOF + +# ---------------------------------------------------------------------- +# run the model! +# ---------------------------------------------------------------------- +cp -p /scratch/snx3000/nponomar/icon-vprm-try2/icon-vprm/bin/icon_bg_vprm_ens icon.exe + +srun ./icon.exe diff --git a/cases/ctdas-test/sbatch_extract_template b/cases/ctdas-test/sbatch_extract_template new file mode 100644 index 00000000..9e8fc606 --- /dev/null +++ b/cases/ctdas-test/sbatch_extract_template @@ -0,0 +1,12 @@ +#!/usr/bin/env bash +#SBATCH --account=s1152 +#SBATCH --job-name="extract_synth" +#SBATCH --partition=normal +#SBATCH --time=00:40:00 +#SBATCH --constraint=mc +#SBATCH --nodes=1 +#SBATCH --ntasks-per-core=1 +#SBATCH --ntasks-per-node=36 +#SBATCH --cpus-per-task=1 + +python {extract_script} diff --git a/jobs/ctdas.py b/jobs/ctdas.py new file mode 100644 index 00000000..fc0d421f --- /dev/null +++ b/jobs/ctdas.py @@ -0,0 +1,70 @@ +#!/usr/bin/env python +# -*- coding: utf-8 -*- +# + +import logging +import os +import subprocess +from . import tools + + +def main(starttime, hstart, hstop, cfg): + + logging.info("Setup the namelist for an CTDAS ICON run and " + "submit the jobs to the queue") + + # Copy icon executable + execname = 'icon.exe' + tools.copy_file(cfg.icon_bin, os.path.join(cfg.icon_work, execname)) + + cfg.ctdas_exec = os.path.join(cfg.ctdas_root, 'exec') + cfg.ctdas_rc = os.path.join(cfg.ctdas_exec, 'ctdas-icon.rc') + cfg.ctdas_job = os.path.join(cfg.ctdas_exec, 'ctdas-icon.jb') + cfg.ctdas_sys_rc = os.path.join(cfg.ctdas_exec, 'da', 'rc', 'cteco2', + 'carbontracker_icon.rc') + cfg.ctdas_case_rc = os.path.join( + cfg.ctdas_exec, 'ctdas-icon-' + cfg.ini_datetime_string + '.rc') + cfg.ctdas_case_sys_rc = os.path.join( + cfg.ctdas_exec, 'da', 'rc', 'cteco2', + 'carbontracker_icon_' + cfg.ini_datetime_string + '.rc') + cfg.ctdas_rc = os.path.join(cfg.ctdas_exec, 'ctdas-icon.rc') + cfg.ctdas_ini_datetime_string = cfg.ini_datetime_string.replace( + 'T', ' ').replace('Z', '') + cfg.ctdas_end_datetime_string = cfg.end_datetime_string.replace( + 'T', ' ').replace('Z', '') + + # Write ctdas_job.rc file + with open(cfg.ctdas_rc) as input_file: + to_write = input_file.read() + output_file = os.path.join(cfg.ctdas_case_rc) + with open(output_file, "w") as outf: + outf.write(to_write.format(cfg=cfg)) + + # Write ctdas_system.rc file + with open(cfg.ctdas_sys_rc) as input_file: + to_write = input_file.read() + output_file = os.path.join(cfg.ctdas_case_sys_rc) + with open(output_file, "w") as outf: + outf.write(to_write.format(cfg=cfg)) + + # Write ctdas.job file + with open(cfg.ctdas_job) as input_file: + to_write = input_file.read() + output_file = os.path.join(cfg.ctdas_root, 'exec', + 'ctdas-icon-' + cfg.ini_datetime_string + '.jb') + with open(output_file, "w") as outf: + outf.write(to_write.format(cfg=cfg)) + + #Create dir for extracted data + tools.create_dir(os.path.join(cfg.icon_base, "extracted"), + "extracted data") + + # Run CTDAS + exitcode = subprocess.call([ + "sbatch", "--wait", + os.path.join(cfg.ctdas_root, 'exec', + 'ctdas-icon-' + cfg.ini_datetime_string + '.jb') + ]) + + if exitcode != 0: + raise RuntimeError("sbatch returned exitcode {}".format(exitcode)) diff --git a/jobs/ctdas_prepare_scripts.py b/jobs/ctdas_prepare_scripts.py new file mode 100644 index 00000000..d52b9d69 --- /dev/null +++ b/jobs/ctdas_prepare_scripts.py @@ -0,0 +1,135 @@ +#!/usr/bin/env python +# -*- coding: utf-8 -*- +# + +import os +from re import A +from . import tools +from datetime import timedelta +import shutil + + +def main(starttime, hstart, hstop, cfg): + + # Set up SBATCH extraction template + tools.create_dir(os.path.join(cfg.icon_work, 'templates'), + "ctdas_icon_templates") + shutil.copyfile( + cfg.ctdas_sbatch_extract_template, + os.path.join(cfg.icon_work, 'templates', 'sbatch_extract_template')) + cfg.ctdas_sbatch_extract_template_scratch = os.path.join( + cfg.icon_work, 'templates', 'sbatch_extract_template') + + # Link the first simulation + restart_first_sim_init = os.path.join( + cfg.ctdas_first_restart_init, + 'ICON-ART-OEM-INIT_' + '%sT00:00:00.000.nc' % + ((starttime + timedelta(hours=24)).strftime('%Y-%m-%d'))) + link_restart_dir = os.path.join( + cfg.icon_base, 'output_%s_opt' % + ((starttime - + timedelta(hours=cfg.ctdas_cycle * 24)).strftime('%Y%m%d%H'))) + restart_file = os.path.join( + link_restart_dir, 'ICON-ART-OEM-INIT_%sT00:00:00.000.nc' % + ((starttime + timedelta(hours=24)).strftime('%Y-%m-%d'))) + tools.create_dir(link_restart_dir, + "restart ini dir for the first CTDAS simulation") + os.system('ln -sf ' + restart_first_sim_init + ' ' + restart_file) + + # Write script files for each ctdas lag + for time in tools.iter_hours(starttime, hstart, hstop, + cfg.ctdas_cycle * 24): + # Set the starting date of the current cycle (filenames are based on this) + current_cycle_filename = time.strftime('%Y%m%d%H') + + # Go back [ctdas-cycle * 24] hours to find the ICON restart file + restart_t = ( + time - timedelta(hours=cfg.ctdas_cycle * 24)).strftime('%Y%m%d%H') + + # Allow for spinup + spinup_time = time + timedelta(seconds=cfg.ctdas_restart_init_time) + + # Limit the simulation end-time if needed + a = spinup_time + timedelta(days=cfg.ctdas_cycle) + if a < starttime + timedelta(hours=hstop - hstart): + end_time = a + else: + end_time = starttime + timedelta(hours=hstop - hstart) + + # Name of the ICBC file + ICBC_filename = os.path.join( + cfg.icon_input_icbc, + time.strftime(cfg.meteo_nameformat) + '.nc') + + # File name for the extraction output + extraction_output_file_name = os.path.join( + cfg.icon_base, 'extracted', + 'output_%s_' % (current_cycle_filename)) + + # Lambdas for emissions & background + lambda_fs = os.path.join(cfg.icon_base, 'input', 'oae', + 'lambda_%s_' % (current_cycle_filename)) + bg_lambda_fs = os.path.join(cfg.icon_base, 'input', 'oae', + 'bg_lambda_%s_' % (current_cycle_filename)) + + # Restart folders & filenames + restart_folder_opt = os.path.join(cfg.icon_base, + 'output_%s_opt' % (restart_t)) + restart_folder_prior = os.path.join( + cfg.icon_base, 'output_%s_priorcycle1' % (restart_t)) + restart_file = 'ICON-ART-OEM-INIT_%sT00:00:00.000.nc' % ( + (spinup_time).strftime('%Y-%m-%d')) + + # Create extraction & runscripts + with open(cfg.ctdas_extract_template) as input_file: + write_extraction = input_file.read() + extract_file = os.path.join(cfg.icon_work, + 'extract_' + current_cycle_filename) + + with open(cfg.ctdas_ICON_template) as input_file: + write_runscript = input_file.read() + runscript_file = os.path.join(cfg.icon_work, + 'runscript_' + current_cycle_filename) + + for case in ['prior', 'priorcycle1', 'opt']: + # Set the case + setattr(cfg, 'suffix', case) + + # Write the extraction script + with open(extract_file + cfg.suffix, "w") as outf: + outf.write( + write_extraction.format( + cfg=cfg, + icon_path=cfg.icon_output + '_' + + time.strftime('%Y%m%d%H') + '_' + cfg.suffix, + fname_base='ICON-ART-UNSTR', + stationdir=cfg.ctdas_datadir, + nneighb=5, + start=spinup_time, + end=end_time, + extracted_file=extraction_output_file_name + + cfg.suffix)) + + # Write the runscript + ICON_output_folder = cfg.icon_output + '_' + time.strftime( + '%Y%m%d%H') + '_' + cfg.suffix + restart_folder = restart_folder_opt if ( + case != 'prior') else restart_folder_prior + with open(runscript_file + cfg.suffix, "w") as outf: + outf.write( + write_runscript.format( + cfg=cfg, + ini_restart_string=time.strftime('%Y-%m-%dT%H:%M:%SZ'), + ini_restart_end_string=end_time.strftime( + '%Y-%m-%dT%H:%M:%SZ'), + inifile=ICBC_filename, + lambda_f=lambda_fs + cfg.suffix + '.nc', + bg_lambda=bg_lambda_fs + cfg.suffix + '.nc', + emissionsgrid=cfg.ctdas_oae_grid.format( + time.strftime(cfg.ctdas_oae_grid_fname)) + if "{}" in cfg.ctdas_oae_grid else cfg.ctdas_oae_grid, + restart_file=os.path.join(restart_folder, + restart_file), + output_directory=ICON_output_folder)) + tools.create_dir(ICON_output_folder, + "creating ICON output folder...")