diff --git a/cases/icon-art-full-chem/config.yaml b/cases/icon-art-full-chem/config.yaml new file mode 100644 index 00000000..3e34efe1 --- /dev/null +++ b/cases/icon-art-full-chem/config.yaml @@ -0,0 +1,91 @@ +# Configuration file for the 'icon-art-full-chem-test' case with ICON + +workflow: icon-art-full-chem +constraint: mc +run_on: cpu +compute_queue: normal +ntasks_per_node: 36 +restart_step: PT6H +startdate: 2019-02-12T00:00:00Z +enddate: 2019-02-26T00:00:00Z + +eccodes_dir: /users/icontest/pool/data/ICON/mch/eccodes_definitions +latbc_filename: LBC_.nc +inidata_prefix: IC_ +inidata_nameformat: '%Y%m%d%H' +inidata_filename_suffix: .nc +lbcdata_prefix: LBC_ +lbcdata_nameformat: '%Y%m%d%H' +lbcdata_filename_suffix: .nc +output_filename: icon-full-chem-test +filename_format: _ +lateral_boundary_grid_order: lateral_boundary +art_input_folder: /scratch/snx3000/ckeller/proc_chain_data/input/icon-art-full-chem/art + +walltime: + prepare_icon: '00:10:00' + prepare_art_full_chem: '01:00:00' + icon: '02:00:00' + +meteo: + dir: /scratch/snx3000/ckeller/proc_chain_data/input/meteo + prefix: era5_ + nameformat: '%Y%m%d%H' + suffix: .grb + inc: 3 + +chem: + dir: /scratch/snx3000/ckeller/proc_chain_data/input/icon-art-full-chem/chem + prefix: camchem_ + nameformat: '%Y%m%d%H' + suffix: .nc + inc: 3 + icbc_filename: camchem-20190212-20190226.nc + ref_date: '2000-03-01 00:00:00' + + +input_files: + radiation_grid_filename: /scratch/snx3000/ckeller/proc_chain_data/input/icon-art-full-chem/grid/base_grid.nc + #radiation_grid_filename: /scratch/snx3000/ckeller/proc_chain_data/input/icon-art-full-chem/grid/icon_R03B07_DOM01.nc + dynamics_grid_filename: /scratch/snx3000/ckeller/proc_chain_data/input/icon-art-full-chem/grid/icon_R03B07_DOM01.nc + #dynamics_grid_filename: /scratch/snx3000/ckeller/proc_chain_data/input/icon-art-full-chem/grid/icon_R03B08_DOM02.nc + map_file_latbc: /scratch/snx3000/ckeller/proc_chain_data/input/icon-art-full-chem/grid/map_file.latbc + lateral_boundary_grid: /scratch/snx3000/ckeller/proc_chain_data/input/icon-art-full-chem/grid/lateral_boundary_DOM01.grid.nc + #lateral_boundary_grid: /scratch/snx3000/ckeller/proc_chain_data/input/icon-art-full-chem/grid/lateral_boundary_DOM02.grid.nc + extpar_filename: /scratch/snx3000/ckeller/proc_chain_data/input/icon-art-full-chem/grid/extpar_DOM01.nc + #extpar_filename: /scratch/snx3000/ckeller/proc_chain_data/input/icon-art-full-chem/grid/extpar_DOM02.nc + cldopt_filename: /scratch/snx3000/ckeller/proc_chain_data/input/icon-art-full-chem/rad/rrtm_cldopt.nc + lrtm_filename: /scratch/snx3000/ckeller/proc_chain_data/input/icon-art-full-chem/rad/rrtmg_lw.nc + map_file_ana: /scratch/snx3000/ckeller/proc_chain_data/input/icon-art-full-chem/mapping/map_file.ana + meccatracer_xml_filename: /scratch/snx3000/ckeller/proc_chain_data/input/icon-art-full-chem/xml/mecca_tracers.xml + oem_gridded_emissions_nc: /scratch/snx3000/ckeller/proc_chain_data/input/icon-art-full-chem/oem/oem_gridded_emissions.nc + oem_vertical_profiles_nc: /scratch/snx3000/ckeller/proc_chain_data/input/icon-art-full-chem/oem/vertical_profiles.nc + oem_hourofday_nc: /scratch/snx3000/ckeller/proc_chain_data/input/icon-art-full-chem/oem/hourofday.nc + oem_dayofweek_nc: /scratch/snx3000/ckeller/proc_chain_data/input/icon-art-full-chem/oem/dayofweek.nc + oem_monthofyear_nc: /scratch/snx3000/ckeller/proc_chain_data/input/icon-art-full-chem/oem/monthofyear.nc + #oem_gridded_emissions_nc: /scratch/snx3000/ckeller/proc_chain_data/input/icon-art-full-chem/oem/oem_gridded_emissions_DOM02.nc + #oem_vertical_profiles_nc: /scratch/snx3000/ckeller/proc_chain_data/input/icon-art-full-chem/oem/vertical_profiles_DOM02.nc + #oem_hourofday_nc: /scratch/snx3000/ckeller/proc_chain_data/input/icon-art-full-chem/oem/hourofday_DOM02.nc + #oem_dayofweek_nc: /scratch/snx3000/ckeller/proc_chain_data/input/icon-art-full-chem/oem/dayofweek_DOM02.nc + #oem_monthofyear_nc: /scratch/snx3000/ckeller/proc_chain_data/input/icon-art-full-chem/oem/monthofyear_DOM02.nc + aerodyn_tracers: /scratch/snx3000/ckeller/proc_chain_data/input/icon-art-full-chem/aerodyn/tracers_aerosol.xml + aerodyn_modes: /scratch/snx3000/ckeller/proc_chain_data/input/icon-art-full-chem/aerodyn/modes.xml + aerodyn_coag: /scratch/snx3000/ckeller/proc_chain_data/input/icon-art-full-chem/aerodyn/coagulate.xml + aerodyn_emiss: /scratch/snx3000/ckeller/proc_chain_data/input/icon-art-full-chem/aerodyn/emiss.xml + + +icon: + binary_file: /scratch/snx3000/ckeller/proc_chain_data/ext/icon-art/bin/icon + runjob_filename: icon_runjob.cfg + era5_icjob: icon_era5_ic.sh + era5_lbcjob: icon_era5_lbc.sh + chem_icjob: icon_camchem_ic.sh + chem_lbcjob: icon_camchem_lbc.sh + compute_queue: normal + walltime: '02:00:00' + np_tot: 32 + np_io: 3 + np_restart: 1 + np_prefetch: 1 + timestep: 120. + diff --git a/cases/icon-art-full-chem/icon_camchem_ic.sh b/cases/icon-art-full-chem/icon_camchem_ic.sh new file mode 100755 index 00000000..05280572 --- /dev/null +++ b/cases/icon-art-full-chem/icon_camchem_ic.sh @@ -0,0 +1,40 @@ +#!/bin/bash + +cd {cfg.icon_input_icbc} + +# --------------------------------- +# -- Pre-processing +# --------------------------------- + +rm -f {cfg.icon_input_icbc}/{cfg.chem_prefix}{cfg.startdate_sim_yyyymmddhh}_ic_icon{cfg.inidata_filename_suffix} + +# -- Change variable and coordinates names to be consistent with MECCA nomenclature +cdo setpartabn,partab_chem,convert {cfg.icon_input_icbc}/{cfg.chem_prefix}{cfg.startdate_sim_yyyymmddhh}_ic{cfg.chem_suffix} data_in.nc + +# --------------------------------- +# -- Re-mapping +# --------------------------------- + +# -- Retrieve the triangular horizontal grid +cdo -s selgrid,2 {cfg.input_files_scratch_dynamics_grid_filename} triangular-grid.nc + +# -- Create the weights for remapping CAM-Chem (lat,lon) grid onto the triangular grid +echo "creating weights" +cdo genbil,triangular-grid.nc data_in.nc weights.nc + +# -- Remap +cdo -s remap,triangular-grid.nc,weights.nc data_in.nc chem_final.nc +rm data_in.nc triangular-grid.nc + +# --------------------------------- +# -- Post-processing +# --------------------------------- + +# -- Rename dimensions and order alphabetically +ncrename -h -d cell,ncells chem_final.nc +ncrename -h -d nv,vertices chem_final.nc +ncks chem_final.nc {cfg.icon_input_icbc}/{cfg.chem_prefix}{cfg.startdate_sim_yyyymmddhh}_ic_icon{cfg.inidata_filename_suffix} +rm chem_final.nc weights.nc + +# -- Clean up +rm {cfg.icon_input_icbc}/{cfg.chem_prefix}{cfg.startdate_sim_yyyymmddhh}_ic{cfg.chem_suffix} diff --git a/cases/icon-art-full-chem/icon_camchem_lbc.sh b/cases/icon-art-full-chem/icon_camchem_lbc.sh new file mode 100755 index 00000000..4f4ff65c --- /dev/null +++ b/cases/icon-art-full-chem/icon_camchem_lbc.sh @@ -0,0 +1,53 @@ +#!/bin/bash + +cd {cfg.icon_input_icbc} + +# -- Loop over file list +i=0 +echo "DATAFILELIST is {datafile_list_chem}" +for datafilename in {datafile_list_chem} ; do + datafile="${{datafilename##*/}}" # get filename without path + outdatafile=${{datafile%.*}} # get filename without suffix + ((i++)) + + # --------------------------------- + # -- Pre-processing + # --------------------------------- + + rm -f {cfg.icon_input_icbc}/${{outdatafile}}_icon{cfg.inidata_filename_suffix} + + # -- Change variable and coordinates names to be consistent with ICON nomenclature + cdo setpartabn,partab_chem,convert $datafilename data_in.nc + + # --------------------------------- + # -- Re-mapping + # --------------------------------- + + # -- Retrieve the lateral boundary grid + cdo -s selgrid,2 {cfg.input_files_scratch_lateral_boundary_grid} triangular-grid.nc + + # -- Create the weights for remapping CAM-Chem data from latlon grid onto the triangular grid + if [[ $i == "1" ]] ; then + echo "creating weights" + cdo genbil,triangular-grid.nc data_in.nc weights.nc + fi + + # -- Remap + cdo -s remap,triangular-grid.nc,weights.nc data_in.nc chem_final.nc + rm data_in.nc triangular-grid.nc + + # --------------------------------- + # -- Post-processing + # --------------------------------- + + # -- Rename dimensions and order alphabetically + ncrename -h -d cell,ncells chem_final.nc + ncrename -h -d nv,vertices chem_final.nc + ncks chem_final.nc {cfg.icon_input_icbc}/${{outdatafile}}_icon{cfg.inidata_filename_suffix} + rm chem_final.nc + + # -- Clean up + rm $datafilename +done +# -- Clean up +rm weights.nc diff --git a/cases/icon-art-full-chem/icon_era5_ic.sh b/cases/icon-art-full-chem/icon_era5_ic.sh new file mode 100644 index 00000000..8744a89d --- /dev/null +++ b/cases/icon-art-full-chem/icon_era5_ic.sh @@ -0,0 +1,132 @@ +#!/bin/bash + +cd {cfg.icon_input_icbc} + +# --------------------------------- +# -- Pre-processing +# --------------------------------- + +rm -f {cfg.icon_input_icbc}/{cfg.meteo_prefix}{cfg.startdate_sim_yyyymmddhh}_ic{cfg.inidata_filename_suffix} + +# -- Convert the GRIB files to NetCDF +cdo -t ecmwf -f nc copy {cfg.meteo_dir}/{cfg.meteo_prefix}{cfg.startdate_sim_yyyymmddhh}{cfg.meteo_suffix} era5_ori.nc + +# -- Change variable and coordinates names to be consistent with ICON nomenclature +cdo setpartabn,partab_meteo,convert era5_ori.nc tmp.nc + +# -- Order the variables alphabetically +ncks tmp.nc data_in.nc +rm tmp.nc era5_ori.nc + +# --------------------------------- +# -- Re-mapping +# --------------------------------- + +# -- Retrieve the triangular grid +cdo -s selgrid,2 {cfg.input_files_scratch_dynamics_grid_filename} triangular-grid.nc + +# -- Create the weights for remapping era5 latlon grid onto the triangular grid +cdo gendis,triangular-grid.nc data_in.nc weights.nc + +# -- Extract the land-sea mask variable +cdo selname,LSM data_in.nc LSM_in.nc +ncrename -h -v LSM,FR_LAND LSM_in.nc +cdo selname,FR_LAND {cfg.input_files_scratch_extpar_filename} LSM_out_tmp.nc + +# -- Add time dimension to LSM_out.nc +ncecat -O -u time LSM_out_tmp.nc LSM_out_tmp.nc +ncks -h -A -v time LSM_in.nc LSM_out_tmp.nc + +# -- Create two different files for land- and sea-mask +# -------------------------------------------------- +# -- Comments: +# -- setctomoss,0. = setting missing values to 0. +# -- gtc = greater than constant (o(t,x) = 1 if i(t,x) > const, o(t,x) = 0 else) +# -- ltc = lower than constant (o(t,x) = 1 if i(t,x) < const, o(t,x) = 0 else) +# -------------------------------------------------- +cdo -L setctomiss,0. -ltc,0.5 LSM_in.nc oceanmask_in.nc +cdo -L setctomiss,0. -gec,0.5 LSM_in.nc landmask_in.nc +cdo -L setctomiss,0. -ltc,0.5 LSM_out_tmp.nc oceanmask_out.nc +cdo -L setctomiss,0. -gec,0.5 LSM_out_tmp.nc landmask_out.nc +cdo setrtoc2,0.5,1.0,1,0 LSM_out_tmp.nc LSM_out.nc +rm LSM_in.nc LSM_out_tmp.nc + +# -- Select surface sea variables defined only on sea +ncks -h -v SST,CI data_in.nc datasea_in.nc + +# -- Select surface variables defined on both that must be remapped differently on sea and on land +ncks -h -v SKT,STL1,STL2,STL3,STL4,ALB_SNOW,W_SNOW,T_SNOW data_in.nc dataland_in.nc + +# ----------------------------------------------------------------------------- +# -- Remap land and ocean area differently for variables +# ----------------------------------------------------------------------------- + +# -- Ocean part +# ----------------- + +# -- Apply the ocean mask (by dividing) +cdo div dataland_in.nc oceanmask_in.nc tmp1_land.nc +cdo div datasea_in.nc oceanmask_in.nc tmp1_sea.nc + +# -- Set missing values to a distance-weighted average +cdo setmisstodis tmp1_land.nc tmp2_land.nc +cdo setmisstodis tmp1_sea.nc tmp2_sea.nc + +# -- Remap +cdo remap,triangular-grid.nc,weights.nc tmp2_land.nc tmp3_land.nc +cdo remap,triangular-grid.nc,weights.nc tmp2_sea.nc tmp3_sea.nc + + +# -- Apply the ocean mask to remapped variables (by dividing) +cdo div tmp3_land.nc oceanmask_out.nc dataland_ocean_out.nc +cdo div tmp3_sea.nc oceanmask_out.nc datasea_ocean_out.nc + +# -- Clean the repository +rm tmp*.nc oceanmask*.nc + +# -- Land part +# ----------------- + +cdo div dataland_in.nc landmask_in.nc tmp1.nc +cdo setmisstodis tmp1.nc tmp2.nc +cdo remap,triangular-grid.nc,weights.nc tmp2.nc tmp3.nc +cdo div tmp3.nc landmask_out.nc dataland_land_out.nc +rm tmp*.nc landmask*.nc dataland_in.nc datasea_in.nc + +# -- merge remapped land and ocean part +# -------------------------------------- + +cdo ifthenelse LSM_out.nc dataland_land_out.nc dataland_ocean_out.nc dataland_out.nc +rm dataland_ocean_out.nc dataland_land_out.nc + +# remap the rest and merge all files +# -------------------------------------- + +# -- Select all variables apart from these ones +ncks -h -x -v SKT,STL1,STL2,STL3,STL4,ALB_SNOW,W_SNOW,T_SNOW,SST,CI,LSM data_in.nc datarest_in.nc +# -- Remap +cdo -s remap,triangular-grid.nc,weights.nc datarest_in.nc era5_final.nc +rm datarest_in.nc data_in.nc + +# -- Fill NaN values for SST and CI +cdo setmisstodis -selname,SST,CI datasea_ocean_out.nc dataland_ocean_out_filled.nc +rm datasea_ocean_out.nc + +# -- Merge remapped files plus land-sea mask from EXTPAR +ncks -h -A dataland_out.nc era5_final.nc +ncks -h -A dataland_ocean_out_filled.nc era5_final.nc +ncks -h -A -v FR_LAND LSM_out.nc era5_final.nc +ncrename -h -v FR_LAND,LSM era5_final.nc +rm LSM_out.nc dataland_out.nc dataland_ocean_out_filled.nc + +# --------------------------------- +# -- Post-processing +# --------------------------------- + +ncrename -h -d cell,ncells era5_final.nc +ncrename -h -d nv,vertices era5_final.nc +ncks era5_final.nc {cfg.icon_input_icbc}/{cfg.meteo_prefix}{cfg.startdate_sim_yyyymmddhh}_ic{cfg.inidata_filename_suffix} +rm era5_final.nc + +# -- Clean up +rm triangular-grid.nc weights.nc diff --git a/cases/icon-art-full-chem/icon_era5_lbc.sh b/cases/icon-art-full-chem/icon_era5_lbc.sh new file mode 100644 index 00000000..fbeb9772 --- /dev/null +++ b/cases/icon-art-full-chem/icon_era5_lbc.sh @@ -0,0 +1,54 @@ +#!/bin/bash + +cd {cfg.icon_input_icbc} + +# -- Loop over file list +i=0 +for datafilename in {datafile_list_meteo} ; do + datafile="${{datafilename##*/}}" # get filename without path + outdatafile=${{datafile%.*}} # get filename without suffix + ((i++)) + + # --------------------------------- + # -- Pre-processing + # --------------------------------- + + rm -f {cfg.icon_input_icbc}/${{outdatafile}}_lbc{cfg.inidata_filename_suffix} + + # -- Convert the GRIB files to NetCDF + cdo -t ecmwf -f nc copy $datafilename era5_ori.nc + + # -- Change variable and coordinates names to be consistent with ICON nomenclature + cdo setpartabn,partab_meteo,convert era5_ori.nc tmp.nc + # -- Order the variables alphabetically + ncks tmp.nc data_in.nc + rm tmp.nc era5_ori.nc + + # --------------------------------- + # -- Re-mapping + # --------------------------------- + + # -- Retrieve the lateral boundary grid + cdo -s selgrid,2 {cfg.input_files_scratch_lateral_boundary_grid} triangular-grid.nc + + # -- Create the weights for remapping era5 latlon grid onto the triangular grid + if [[ $i == "1" ]] ; then + echo "creating weights" + cdo gendis,triangular-grid.nc data_in.nc weights.nc + fi + + # -- Remap + cdo -s remap,triangular-grid.nc,weights.nc data_in.nc era5_final.nc + rm data_in.nc triangular-grid.nc + + # --------------------------------- + # -- Post-processing + # --------------------------------- + + ncrename -h -d cell,ncells era5_final.nc + ncrename -h -d nv,vertices era5_final.nc + ncks era5_final.nc {cfg.icon_input_icbc}/${{outdatafile}}_lbc{cfg.inidata_filename_suffix} + rm era5_final.nc +done +# -- Clean up +rm weights.nc \ No newline at end of file diff --git a/cases/icon-art-full-chem/icon_runjob.cfg b/cases/icon-art-full-chem/icon_runjob.cfg new file mode 100644 index 00000000..e01813e1 --- /dev/null +++ b/cases/icon-art-full-chem/icon_runjob.cfg @@ -0,0 +1,412 @@ +#!/usr/bin/env bash +#SBATCH --job-name=icon-full-chem +#SBATCH --account={cfg.compute_account} +#SBATCH --time={cfg.walltime_icon} +#SBATCH --nodes={cfg.icon_np_tot} +#SBATCH --ntasks-per-node={cfg.ntasks_per_node} +#SBATCH --partition={cfg.compute_queue} +#SBATCH --constraint={cfg.constraint} +#SBATCH --hint=nomultithread +#SBATCH --output={cfg.logfile} +#SBATCH --chdir={cfg.icon_work} + +# OpenMP environment variables +# ---------------------------- +export OMP_NUM_THREADS=1 +export ICON_THREADS=1 +export OMP_SCHEDULE=static,12 +export OMP_DYNAMIC="false" +export OMP_STACKSIZE=200M + +set -e -x + +# -- ECCODES path +export ECCODES_DEFINITION_PATH={cfg.eccodes_dir} + +# ---------------------------------------------------------------------------- +# 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 = {cfg.lrestart} ! .TRUE.=current experiment is resumed + read_restart_namelists = .TRUE. +/ + +! master_time_control_nml: --------------------------------------------------- +&master_time_control_nml + calendar = 'proleptic gregorian' + restartTimeIntval = '{cfg.restart_step}' + checkpointTimeIntval = '{cfg.restart_step}' + experimentStartDate = '{cfg.ini_datetime_string}' + experimentStopDate = '{cfg.end_datetime_string}' +/ + +! 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_{cfg.casename}" ! 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 +/ +EOF + +# ---------------------------------------------------------------------- +# model namelists +# ---------------------------------------------------------------------- + +cat > NAMELIST_{cfg.casename} << 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 = {cfg.icon_np_io} ! number of I/O processors + num_restart_procs = {cfg.icon_np_restart} ! number of restart processors + num_prefetch_proc = {cfg.icon_np_prefetch} ! 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 = {cfg.icon_timestep} ! timestep in seconds + ldynamics = .TRUE. ! compute adiabatic dynamic tendencies + ltransport = .TRUE. ! compute large-scale tracer transport + ntracer = 6 ! 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 + restart_filename = "{cfg.icon_restart_out}/{cfg.output_filename}_.nc" +/ + +! art_nml: Aerosols and Reactive Trace gases extension------------------------------------------------- +&art_nml + lart_chem = .TRUE. ! enables chemistry + lart_mecca = .TRUE. ! enables MECCA chemistry + lart_aerosol = .TRUE. ! main switch for the treatment of atmospheric aerosol + lart_chemtracer = .FALSE. ! main switch for the treatment of chemical tracer + lart_diag_out = .FALSE. ! If this switch is set to .TRUE., diagnostic + ! ... output elds are available. Set it to + ! ... .FALSE. when facing memory problems. + iart_init_gas = 4 + cart_cheminit_file = '{inidata_filename}' + cart_cheminit_type = 'EMAC' + cart_cheminit_coord = '{inidata_filename}' + cart_mecca_xml = '{cfg.input_files_scratch_meccatracer_xml_filename}' ! path to xml file for passive tracers + cart_input_folder = '{cfg.art_input_folder}' ! absolute Path to ART source code + iart_init_aero = 4 + cart_aero_emiss_xml = '{cfg.input_files_scratch_aerodyn_emiss}' + cart_aerosol_xml = '{cfg.input_files_scratch_aerodyn_tracers}' ! Path to XML file for aerosol tracers + cart_modes_xml = '{cfg.input_files_scratch_aerodyn_modes}' ! Path to XML file for modes + cart_coag_xml = '{cfg.input_files_scratch_aerodyn_coag}' ! Path to XML file for coagulation processes + iart_modeshift = 1 ! Shift2larger + iart_aero_washout = 1 ! Washout of aerosols + iart_isorropia = 1 ! Gas-aerosol partitioning +/ + +! oem_nml: online emission module --------------------------------------------- +&oemctrl_nml + gridded_emissions_nc = '{cfg.input_files_scratch_oem_gridded_emissions_nc}' + vertical_profile_nc = '{cfg.input_files_scratch_oem_vertical_profiles_nc}' + hour_of_day_nc = '{cfg.input_files_scratch_oem_hourofday_nc}' + day_of_week_nc = '{cfg.input_files_scratch_oem_dayofweek_nc}' + month_of_year_nc = '{cfg.input_files_scratch_oem_monthofyear_nc}' +/ + +! 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.input_files_scratch_extpar_filename}' ! 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 = "{inidata_filename}" ! initial data filename + ana_varnames_map_file = "{cfg.input_files_scratch_map_file_ana}" ! 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.input_files_scratch_dynamics_grid_filename}" ! array of the grid filenames for the dycore + radiation_grid_filename = "{cfg.input_files_scratch_radiation_grid_filename}" ! 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 = 8 ! width of lateral boundary nudging zone + support_baryctr_intp = .FALSE. ! barycentric interpolation support for output + nudge_max_coeff = 0.07 + 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 = 60 ! Number of vertical levels in boundary data + latbc_boundary_grid = "{cfg.input_files_scratch_lateral_boundary_grid}" ! Grid file defining the lateral boundary + latbc_path = "{cfg.icon_input_icbc}" ! Absolute path to boundary data + latbc_varnames_map_file = "{cfg.input_files_scratch_map_file_latbc}" + 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.input_files_scratch_lrtm_filename}' ! longwave absorption coefficients for RRTM_LW + cldopt_filename = '{cfg.input_files_scratch_cldopt_filename}' ! RRTM cloud optical properties + dt_rad = 720. ! time step for radiation in s + dt_conv = 120.,90.,90. ! time step for convection in s (domain specific) + dt_sso = 120.,360.,360. ! time step for SSO parameterization + dt_gwd = 360.,360.,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 = 5 ! 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 + include_last = .TRUE. + remap = 0 + output_filename = 'icon-art-full-chem_EU' + filename_format = '{cfg.icon_output}/_LBC_' ! file name base + output_grid = .TRUE. + ml_varlist = 'U','V','W', + 'PRES','TEMP', + 'QV','QC','QI', + 'QR','QS', + 'z_ifc','z_mc', + 'group:ART_CHEMISTRY', + 'group:ART_AEROSOL', + 'group:LAND_VARS', + 'geopot', + 'topography_c', + 'group:ATMO_ML_VARS' , + 'group:PBL_VARS', + 'group:PRECIP_VARS', + 'group:NH_PROG_VARS' +/ + +&output_nml + filetype = 5 ! output format: 2=GRIB2, 4=NETCDFv2 + dom = 1 ! write domain 1 only + output_bounds = 0., 0., 10000. ! start, end, increment + steps_per_file = 1 ! number of steps per file + include_last = .TRUE. + remap = 0 + output_filename = 'icon-art-full-chem_EU' + filename_format = '{cfg.icon_output}/_IC_' ! file name base + output_grid = .TRUE. + ml_varlist = 'u','v','w', + 'pres','temp', + 'qv','qc','qi', + 'qr','qs','tke', + 't_g','t_ice', + 'h_ice','t_mnw_lk', + 't_wml_lk','h_ml_lk', + 't_bot_lk','c_t_lk', + 'qv_s','w_i','t_so', + 'gz0','alb_si', + 'smi','w_so_ice', + 't_snow','w_snow', + 'rho_snow','h_snow', + 'freshsnow','fr_seaice', + 'fr_land','z_ifc', + 'group:ART_CHEMISTRY', + 'group:ART_AEROSOL' +/ + +! 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 +/ + +! 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, 22 ! 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! +# ---------------------------------------------------------------------- +handle_error(){{ + # Check for invalid pointer error at the end of icon-art + if grep -q "free(): invalid pointer" {cfg.logfile} && grep -q "clean-up finished" {cfg.logfile}; then + exit 0 + else + exit 1 + fi +}} +srun ./{cfg.icon_execname} || handle_error diff --git a/cases/icon-art-full-chem/partab_chem b/cases/icon-art-full-chem/partab_chem new file mode 100644 index 00000000..a7f70058 --- /dev/null +++ b/cases/icon-art-full-chem/partab_chem @@ -0,0 +1,507 @@ +¶meter +name = ALKNIT +out_name = ALKNO3_full +units = "mol mol-1" +/ +¶meter +name = BCARY +out_name = BCARY_full +units = "mol mol-1" +/ +¶meter +name = BENZENE +out_name = BENZENE_full +units = "mol mol-1" +/ +¶meter +name = BIGALD +out_name = C5H6O2_full +units = "mol mol-1" +/ +¶meter +name = BIGALD1 +out_name = BIGALD1_full +units = "mol mol-1" +/ +¶meter +name = BIGALD2 +out_name = BIGALD2_full +units = "mol mol-1" +/ +¶meter +name = BIGALD3 +out_name = BIGALD3_full +units = "mol mol-1" +/ +¶meter +name = BIGALD4 +out_name = BIGALD4_full +units = "mol mol-1" +/ +¶meter + name = BIGALK + out_name = C5H12_full + units = "mol mol-1" +/ +¶meter + name = BIGENE + out_name = C4H8_full + units = "mol mol-1" +/ +¶meter +name = C2H2 +out_name = C2H2_full +units = "mol mol-1" +/ +¶meter +name = C2H4 +out_name = C2H4_full +units = "mol mol-1" +/ +¶meter +name = C2H5OH +out_name = C2H5OH_full +units = "mol mol-1" +/ +¶meter + name = C2H6 + out_name = C2H6_full + units = "mol mol-1" +/ +¶meter + name = C3H6 + out_name = C3H6_full + units = "mol mol-1" +/ +¶meter + name = C3H8 + out_name = C3H8_full + units = "mol mol-1" +/ +¶meter + name = CH2O + out_name = HCHO_full + units = "mol mol-1" +/ +¶meter + name = CH3CHO + out_name = CH3CHO_full + units = "mol mol-1" +/ +¶meter + name = CH3COCH3 + out_name = CH3COCH3_full + units = "mol mol-1" +/ +¶meter + name = CH3COCHO + out_name = MGLYOX_full + units = "mol mol-1" +/ +¶meter + name = CH3COOH + out_name = CH3CO2H_full + units = "mol mol-1" +/ +¶meter + name = CH3OH + out_name = CH3OH_full + units = "mol mol-1" +/ +¶meter + name = CH3OOH + out_name = CH3OOH_full + units = "mol mol-1" +/ +¶meter + name = CH4 + out_name = CH4_full + units = "mol mol-1" +/ +¶meter + name = CO + out_name = CO_full + units = "mol mol-1" +/ +¶meter + name = CRESOL + out_name = CRESOL_full + units = "mol mol-1" +/ +¶meter + name = DMS + out_name = DMS_full + units = "mol mol-1" +/ +¶meter + name = GLYOXAL + out_name = GLYOX_full + units = "mol mol-1" +/ +¶meter + name = H2O + out_name = H2O_full + units = "mol mol-1" +/ +¶meter + name = H2O2 + out_name = H2O2_full + units = "mol mol-1" +/ +¶meter + name = HCOOH + out_name = HCOOH_full + units = "mol mol-1" +/ +¶meter + name = HNO3 + out_name = HNO3_full + units = "mol mol-1" +/ +¶meter + name = HO2 + out_name = HO2_full + units = "mol mol-1" +/ +¶meter + name = HO2NO2 + out_name = HNO4_full + units = "mol mol-1" +/ +¶meter + name = HONITR + out_name = LHONITR_full + units = "mol mol-1" +/ +¶meter + name = HYAC + out_name = ACETOL_full + units = "mol mol-1" +/ +¶meter + name = ISOP + out_name = C5H8_full + units = "mol mol-1" +/ +¶meter + name = ISOPNITA + out_name = LISOPBDNO3_full + units = "mol mol-1" +/ +¶meter + name = ISOPNITB + out_name = LISOPACNO3_full + units = "mol mol-1" +/ +¶meter + name = MACR + out_name = MACR_full + units = "mol mol-1" +/ +¶meter + name = MEK + out_name = MEK_full + units = "mol mol-1" +/ +¶meter + name = MPAN + out_name = MPAN_full + units = "mol mol-1" +/ +¶meter + name = MTERP + out_name = LTERP_full + units = "mol mol-1" +/ +¶meter + name = MVK + out_name = MVK_full + units = "mol mol-1" +/ +¶meter + name = N2O + out_name = N2O_full + units = "mol mol-1" +/ +¶meter + name = N2O5 + out_name = N2O5_full + units = "mol mol-1" +/ +¶meter + name = NH3 + out_name = NH3_full + units = "mol mol-1" +/ +¶meter + name = NO2 + out_name = NO2_full + units = "mol mol-1" +/ +¶meter + name = NO3 + out_name = NO3_full + units = "mol mol-1" +/ +¶meter + name = NO + out_name = NO_full + units = "mol mol-1" +/ +¶meter + name = NOA + out_name = NOA_full + units = "mol mol-1" +/ +¶meter + name = O3 + out_name = O3_full + units = "mol mol-1" +/ +¶meter + name = OH + out_name = OH_full + units = "mol mol-1" +/ +¶meter + name = ONITR + out_name = LONITR_full + units = "mol mol-1" +/ +¶meter + name = PAN + out_name = PAN_full + units = "mol mol-1" +/ +¶meter + name = PBZNIT + out_name = PBZNIT_full + units = "mol mol-1" +/ +¶meter + name = PHENOL + out_name = PHENOL_full + units = "mol mol-1" +/ +¶meter + name = SO2 + out_name = SO2_full + units = "mol mol-1" +/ +¶meter + name = TERPNIT + out_name = BPINANO3_full + units = "mol mol-1" +/ +¶meter + name = TOLUENE + out_name = TOLUENE_full + units = "mol mol-1" +/ +¶meter + name = XYLENES + out_name = LXYL_full + units = "mol mol-1" +/ +¶meter + name = poa_mixed_ait + out_name = poa_mixed_ait + units = "kg kg-1" +/ +¶meter + name = poa_insol_ait + out_name = poa_insol_ait + units = "kg kg-1" +/ +¶meter + name = poa_mixed_acc + out_name = poa_mixed_acc + units = "kg kg-1" +/ +¶meter + name = poa_insol_acc + out_name = poa_insol_acc + units = "kg kg-1" +/ +¶meter + name = soot_mixed_ait + out_name = soot_mixed_ait + units = "kg kg-1" +/ +¶meter + name = soot_insol_ait + out_name = soot_insol_ait + units = "kg kg-1" +/ +¶meter + name = soot_mixed_acc + out_name = soot_mixed_acc + units = "kg kg-1" +/ +¶meter + name = soot_insol_acc + out_name = soot_insol_acc + units = "kg kg-1" +/ +¶meter + name = so4_mixed_ait + out_name = so4_mixed_ait + units = "kg kg-1" +/ +¶meter + name = so4_sol_ait + out_name = so4_sol_ait + units = "kg kg-1" +/ +¶meter + name = nh4_mixed_ait + out_name = nh4_mixed_ait + units = "kg kg-1" +/ +¶meter + name = nh4_sol_ait + out_name = nh4_sol_ait + units = "kg kg-1" +/ +¶meter + name = so4_mixed_acc + out_name = so4_mixed_acc + units = "kg kg-1" +/ +¶meter + name = so4_sol_acc + out_name = so4_sol_acc + units = "kg kg-1" +/ +¶meter + name = nh4_mixed_acc + out_name = nh4_mixed_acc + units = "kg kg-1" +/ +¶meter + name = nh4_sol_acc + out_name = nh4_sol_acc + units = "kg kg-1" +/ +¶meter + name = so4_mixed_coa + out_name = so4_mixed_coa + units = "kg kg-1" +/ +¶meter + name = so4_sol_coa + out_name = so4_sol_coa + units = "kg kg-1" +/ +¶meter + name = nh4_mixed_coa + out_name = nh4_mixed_coa + units = "kg kg-1" +/ +¶meter + name = nh4_sol_coa + out_name = nh4_sol_coa + units = "kg kg-1" +/ +¶meter + name = na_mixed_ait + out_name = na_mixed_ait + units = "kg kg-1" +/ +¶meter + name = na_sol_ait + out_name = na_sol_ait + units = "kg kg-1" +/ +¶meter + name = cl_mixed_ait + out_name = cl_mixed_ait + units = "kg kg-1" +/ +¶meter + name = cl_sol_ait + out_name = cl_sol_ait + units = "kg kg-1" +/ +¶meter + name = na_mixed_acc + out_name = na_mixed_acc + units = "kg kg-1" +/ +¶meter + name = na_sol_acc + out_name = na_sol_acc + units = "kg kg-1" +/ +¶meter + name = cl_mixed_acc + out_name = cl_mixed_acc + units = "kg kg-1" +/ +¶meter + name = cl_sol_acc + out_name = cl_sol_acc + units = "kg kg-1" +/ +¶meter + name = na_mixed_coa + out_name = na_mixed_coa + units = "kg kg-1" +/ +¶meter + name = na_sol_coa + out_name = na_sol_coa + units = "kg kg-1" +/ +¶meter + name = cl_mixed_coa + out_name = cl_mixed_coa + units = "kg kg-1" +/ +¶meter + name = cl_sol_coa + out_name = cl_sol_coa + units = "kg kg-1" +/ +¶meter + name = dust_mixed_ait + out_name = dust_mixed_ait + units = "kg kg-1" +/ +¶meter + name = dust_insol_ait + out_name = dust_insol_ait + units = "kg kg-1" +/ +¶meter + name = dust_mixed_acc + out_name = dust_mixed_acc + units = "kg kg-1" +/ +¶meter + name = dust_insol_acc + out_name = dust_insol_acc + units = "kg kg-1" +/ +¶meter + name = dust_mixed_coa + out_name = dust_mixed_coa + units = "kg kg-1" +/ +¶meter + name = dust_insol_coa + out_name = dust_insol_coa + units = "kg kg-1" +/ +¶meter + name = PS + out_name = PS + units = "Pa" + long_name = "Surface pressure" +/ +¶meter + name = Q + out_name = Q + units = "kg kg-1" + long_name = "Specific humidity" +/ \ No newline at end of file diff --git a/cases/icon-art-full-chem/partab_meteo b/cases/icon-art-full-chem/partab_meteo new file mode 100644 index 00000000..4d6986ea --- /dev/null +++ b/cases/icon-art-full-chem/partab_meteo @@ -0,0 +1,185 @@ +¶meter + name = u + out_name = U + standard_name = eastward_wind + long_name = "horizontal velocity component" + units = "m s-1" +/ +¶meter + name = v + out_name = V + standard_name = northward_wind + long_name = "horizontal velocity component" + units = "m s-1" +/ +¶meter + name = w + out_name = W + standard_name = lagrangian_tendency_of_air_pressure + long_name = "vertical velocity" + units = "Pa s-1" +/ +¶meter + name = t + out_name = T + standard_name = air_temperature + long_name = "temperature" + units = "K" +/ +¶meter + name = q + out_name = QV + standard_name = specific_humidity + long_name = "specific humidity" + units = "kg kg-1" +/ +¶meter + name = clwc + out_name = QC + long_name = "specific cloud liquid water content" + units = "kg kg-1" +/ +¶meter + name = ciwc + out_name = QI + long_name = "specific cloud ice water content" + units = "kg kg-1" +/ +¶meter + name = crwc + out_name = QR + long_name = "specific rain water content" + units = "kg kg-1" +/ +¶meter + name = cswc + out_name = QS + long_name = "specific snow water content" + units = "kg kg-1" +/ +¶meter + name = lnsp + out_name = LNPS + long_name = "logarithm of surface pressure" +/ +¶meter + name = Z + out_name = GEOP_SFC + standard_name = geopotential + long_name = "surface geopotential" + units = "m2 s-2" +/ +¶meter + name = z + out_name = GEOP_ML + standard_name = geopotential + long_name = "model level geopotential" + units = "m2 s-2" +/ +¶meter + name = TSN + out_name = T_SNOW + standard_name = temperature_in_surface_snow + long_name = "temperature of snow layer" + units = "K" +/ +¶meter + name = RSN + out_name = RHO_SNOW + long_name = "snow density" + unit = "kg m-3" +/ +¶meter + name = ASN + out_name = ALB_SNOW + long_name = "snow albedo" + unit = "(0 - 1)" +/ +¶meter + name = SD + out_name = W_SNOW + standard_name = lwe_thickness_of_surface_snow_amount + long_name = "snow depth" + unit = "m of water equivalent" +/ +¶meter + name = SKT + out_name = SKT + long_name = "skin temperature" + unit = "K" +/ +¶meter + name = SSTK + out_name = SST + long_name = "sea surface temperature" + unit = "K" +/ +¶meter + name = SRC + out_name = W_I + long_name = "skin reservoir content" + unit = "m of water equivalent" +/ +¶meter + name = LSM + out_name = LSM + standard_name = land_binary_mask + long_name = "land-sea mask" + unit = "(0 - 1)" +/ +¶meter + name = CI + out_name = CI + standard_name = sea_ice_area_fraction + long_name = "sea ice cover" + unit = "(0 - 1)" +/ +¶meter + name = STL1 + out_name = STL1 + standard_name = surface_temperature + long_name = "soil temperature level 1" + unit = "K" +/ +¶meter + name = STL2 + out_name = STL2 + long_name = "soil temperature level 2" + unit = "K" +/ +¶meter + name = STL3 + out_name = STL3 + long_name = "soil temperature level 3" + unit = "K" +/ +¶meter + name = STL4 + out_name = STL4 + long_name = "soil temperature level 4" + unit = "K" +/ +¶meter + name = SWVL1 + out_name = SMIL1 + long_name = "soil moisture index layer 1" + unit = "m3 m-3" +/ +¶meter + name = SWVL2 + out_name = SMIL2 + long_name = "soil moisture index layer 2" + unit = "m3 m-3" +/ +¶meter + name = SWVL3 + out_name = SMIL3 + long_name = "soil moisture index layer 3" + unit = "m3 m-3" +/ +¶meter + name = SWVL4 + out_name = SMIL4 + long_name = "soil moisture index layer 4" + unit = "m3 m-3" +/ diff --git a/jobs/__init__.py b/jobs/__init__.py index 332f34a8..a3d75541 100644 --- a/jobs/__init__.py +++ b/jobs/__init__.py @@ -22,3 +22,4 @@ from . import prepare_icon from . import reduce_output from . import verify_chain +from . import prepare_art_full_chem diff --git a/jobs/prepare_art_full_chem.py b/jobs/prepare_art_full_chem.py new file mode 100644 index 00000000..ce9a194c --- /dev/null +++ b/jobs/prepare_art_full_chem.py @@ -0,0 +1,439 @@ +#!/usr/bin/env python +# -*- coding: utf-8 -*- + +import os +import logging +import xarray as xr +from . import tools, prepare_icon +import shutil +import subprocess +from datetime import datetime +from .tools.camchem_interpolation import vert_intpl, time_intpl, extract_timeslice +from .tools.camchem_ic_lbc_reggrid import process_lbc, process_ic + +BASIC_PYTHON_JOB = True + +def main(cfg): + """ + Prepare ICON-ART simulations with full chemistry. + + Parameters + ---------- + cfg : Config + Object holding all user-configuration parameters as attributes. + """ + prepare_icon.set_cfg_variables(cfg) + tools.change_logfile(cfg.logfile) + logging.info("Prepare ICON-ART for full chemistry simulation.") + + ## --------------------------- + ## -- Meteorological fields -- + ## --------------------------- + + # -- Copy partab_meteo in workdir + shutil.copy( + os.path.join(cfg.case_path, 'partab_meteo'), + os.path.join(cfg.icon_input_icbc, 'partab_meteo')) + + # -- Create initial conditions on ICON grid + if cfg.lrestart == '.FALSE.': + + # -- Copy ERA5 processing script in workdir + with open( + os.path.join(cfg.case_path, 'icon_era5_ic.sh') + ) as inf: + to_write = inf.read() + output_file = os.path.join(cfg.icon_input_icbc, 'icon_era5_ic.sh') + with open(output_file, "w") as outf: + outf.write(to_write.format(cfg=cfg)) + + # -- Run ERA5 processing script + process = subprocess.Popen([ + "bash", + os.path.join(cfg.icon_input_icbc, 'icon_era5_ic.sh') + ], + stdout=subprocess.PIPE) + process.communicate() + + # -- Create lateral boundary conditions on ICON grid + + # -- Collect file list + datafile_list_meteo = [] + for time in tools.iter_hours(cfg.startdate_sim, cfg.enddate_sim, + cfg.meteo_inc): + meteo_file = os.path.join(cfg.meteo_dir, + cfg.meteo_prefix + time.strftime(cfg.meteo_nameformat)) + datafile_list_meteo.append(str(meteo_file) + cfg.meteo_suffix) + datafile_list_meteo = ' '.join([str(v) for v in datafile_list_meteo]) + + # -- Copy ERA5 processing script in workdir + with open( + os.path.join(cfg.case_path, 'icon_era5_lbc.sh') + ) as inf: + to_write = inf.read() + output_file = os.path.join(cfg.icon_input_icbc, 'icon_era5_lbc.sh') + with open(output_file, "w") as outf: + outf.write(to_write.format(cfg=cfg, + datafile_list_meteo = datafile_list_meteo)) + + # -- Run ERA5 processing script + process = subprocess.Popen([ + "bash", + os.path.join(cfg.icon_input_icbc, 'icon_era5_lbc.sh') + ], + stdout=subprocess.PIPE) + process.communicate() + + + ## ---------------------- + ## -- Chemistry fields -- + ## ---------------------- + + # -- Define CAM-Chem (CESM2.1) fields for interpolation + + # -- Define aerosol tracers and split of aerosols into modes + MW_NH4HSO4 = 115.11 + MW_SO4 = 96.06 + MW_CL = 35.453 + MW_NA = 23.0 + MW_NACL = 58.453 + + aero_mode_dict = { + 'pom_a4':{ + 'poa_mixed_ait':0.9, # 90%-->mixed, 10%-->sol/insol + 'poa_insol_ait':0.1 + }, + 'pom_a1':{ + 'poa_mixed_acc':0.9, + 'poa_insol_acc':0.1 + }, + 'bc_a4':{ + 'soot_mixed_ait':0.9, + 'soot_insol_ait':0.1 + }, + 'bc_a1':{ + 'soot_mixed_acc':0.9, + 'soot_insol_acc':0.1 + }, + 'so4_a2':{ + 'so4_mixed_ait':0.9*MW_SO4/MW_NH4HSO4, # Sulfat aerosol in CAM-Chem has composition NH4HSO4 + 'so4_sol_ait':0.1*MW_SO4/MW_NH4HSO4 + }, + 'so4_a1':{ + 'so4_mixed_acc':0.9*MW_SO4/MW_NH4HSO4, + 'so4_sol_acc':0.1*MW_SO4/MW_NH4HSO4 + }, + 'so4_a3':{ + 'so4_mixed_coa':0.9*MW_SO4/MW_NH4HSO4, + 'so4_sol_coa':0.1*MW_SO4/MW_NH4HSO4 + }, + 'ncl_a2':{ + 'na_mixed_ait':0.9*MW_NA/MW_NACL, + 'na_sol_ait':0.1*MW_NA/MW_NACL, + 'cl_mixed_ait':0.9*MW_CL/MW_NACL, + 'cl_sol_ait':0.1*MW_CL/MW_NACL + }, + 'ncl_a1':{ + 'na_mixed_acc':0.9*MW_NA/MW_NACL, + 'na_sol_acc':0.1*MW_NA/MW_NACL, + 'cl_mixed_acc':0.9*MW_CL/MW_NACL, + 'cl_sol_acc':0.1*MW_CL/MW_NACL}, + 'ncl_a3':{ + 'na_mixed_coa':0.9*MW_NA/MW_NACL, + 'na_sol_coa':0.1*MW_NA/MW_NACL, + 'cl_mixed_coa':0.9*MW_CL/MW_NACL, + 'cl_sol_coa':0.1*MW_CL/MW_NACL}, + 'dst_a2':{ + 'dust_mixed_ait':0.9, + 'dust_insol_ait':0.1 + }, + 'dst_a1':{ + 'dust_mixed_acc':0.9, + 'dust_insol_acc':0.1 + }, + 'dst_a3':{ + 'dust_mixed_coa':0.9, + 'dust_insol_coa':0.1 + }, + 'NH4':{ + 'nh4_mixed_acc':0.9*0.9*0.9, # Ammonium aerosol 90%-->Aitken/Accumulation, + 'nh4_sol_acc':0.9*0.9*0.1, # 10%-->Coarse, 90%-->Accumulation, 10%-->Aitken + 'nh4_mixed_ait':0.9*0.1*0.9, + 'nh4_sol_ait':0.9*0.1*0.1, + 'nh4_mixed_coa':0.1*0.9, + 'nh4_sol_coa':0.1*0.1 + } + } + + # -- Define chemical tracers and molar weights for VMR --> MMR + chem_mw_dict = { + 'ALKNIT':133.1460, + 'BENZENE':78.1121, + 'BIGALD':98.1001, + 'BIGALD1':84.0735, + 'BIGALD2':98.1001, + 'BIGALD3':98.1001, + 'BIGALD4':112.1268, + 'BIGALK':72.1490, + 'BIGENE':56.1065, + 'C2H2':26.0374, + 'C2H4':28.0532, + 'C2H5OH':46.0685, + 'C2H6':30.0691, + 'C3H6':42.0799, + 'C3H8':44.0957, + 'CH2O':30.0260, + 'CH3CHO':44.0526, + 'CH3COCH3':58.0793, + 'CH3COCHO':72.0628, + 'CH3COOH':60.0520, + 'CH3OH':32.0419, + 'CH3OOH':48.0413, + 'CH4':16.0425, + 'CO':28.0101, + 'CRESOL':108.1381, + 'DMS':62.1340, + 'GLYOXAL':26.0374, + 'H2O':18.0153, + 'H2O2':34.0147, + 'HCOOH':46.0254, + 'HNO3':63.0128, + 'HO2':33.0068, + 'HO2NO2':79.0123, + 'HONITR':135.1188, + 'HYAC':74.0787, + 'ISOP':68.1172, + 'ISOPNITA':147.1295, + 'ISOPNITB':147.1295, + 'MACR':70.0900, + 'MEK':72.1059, + 'MPAN':147.0864, + 'MVK':70.0900, + 'N2O':44.0129, + 'N2O5':108.0105, + 'NH3':17.0306, + 'NO':30.0061, + 'NO2':46.0056, + 'NO3':62.0050, + 'NOA':119.0762, + 'O3':47.9982, + 'OH':17.0073, + 'ONITR':133.1029, + 'PAN':121.0491, + 'PBZNIT':183.1186, + 'PHENOL':94.1115, + 'SO2':64.0648, + 'TERPNIT':215.2467, + 'TOLUENE':92.1387, + 'XYLENES':106.1653 + } + + camchem_spec = list(aero_mode_dict.keys()) + list(chem_mw_dict.keys()) + + chem_filename = os.path.join(cfg.chem_dir, + cfg.chem_icbc_filename) + meteo_filename = os.path.join(cfg.icon_input_icbc, + cfg.meteo_prefix + + cfg.startdate_sim.strftime(cfg.meteo_nameformat) + + '_lbc' + + cfg.inidata_filename_suffix) + tmp_filename = os.path.join(cfg.icon_input_icbc, + 'camchem_temp.nc') + + # -- Vertically interpolate CAM-Chem fields + vert_intpl(chem_filename, + meteo_filename, + tmp_filename, + camchem_spec, + cfg.startdate_sim, + cfg.enddate_sim, + cfg.chem_ref_date) + + # -- Create temporary folder for time slices of chem data + tmp_dir = os.path.join(cfg.icon_input_icbc,'chem_data_tmp') + if os.path.exists(tmp_dir): + os.rmdir(tmp_dir) + os.makedirs(tmp_dir) + + # -- Extract one NetCDF file per time step + out_template = os.path.join(tmp_dir, + cfg.chem_prefix + + cfg.chem_nameformat + + cfg.chem_suffix) + + extract_timeslice(tmp_filename, + out_template, + camchem_spec, + cfg.chem_ref_date) + + # -- Remove tmp file with vertically interpolated fields + os.remove(tmp_filename) + + # -- Linear interpolation with respect to time + time_intpl(tmp_dir, + cfg.chem_prefix) + + # -- Prepare data for initial conditions + process_ic(tmp_dir, + cfg.icon_input_icbc, + cfg.startdate_sim, + aero_mode_dict, + cfg.chem_prefix) + + # -- Prepare data for lateral boundary conditions + process_lbc(tmp_dir, + cfg.icon_input_icbc, + cfg.startdate_sim, + cfg.enddate_sim, + cfg.chem_inc, + cfg.chem_prefix, + cfg.chem_suffix, + cfg.chem_nameformat, + chem_mw_dict, + aero_mode_dict) + + # -- Remove tmp folder with chem data + shutil.rmtree(tmp_dir) + + # -- Copy partab_chem in workdir + shutil.copy( + os.path.join(cfg.case_path, 'partab_chem'), + os.path.join(cfg.icon_input_icbc, 'partab_chem')) + + # -- Create initial conditions on ICON grid + if cfg.lrestart == '.FALSE.': + # -- Copy CAM-Chem processing script in workdir + with open( + os.path.join(cfg.case_path, 'icon_camchem_ic.sh') + ) as inf: + to_write = inf.read() + output_file = os.path.join(cfg.icon_input_icbc, 'icon_camchem_ic.sh') + with open(output_file, "w") as outf: + outf.write(to_write.format(cfg=cfg)) + + # -- Run CAM-Chem processing script + process = subprocess.Popen([ + "bash", + os.path.join(cfg.icon_input_icbc, 'icon_camchem_ic.sh') + ], + stdout=subprocess.PIPE) + process.communicate() + + # -- Create lateral boundary conditions on ICON grid + + # -- Collect file list + datafile_list_chem = [] + for time in tools.iter_hours(cfg.startdate_sim, cfg.enddate_sim, + cfg.meteo_inc): + chem_file = os.path.join(cfg.icon_input_icbc, + cfg.chem_prefix + time.strftime(cfg.chem_nameformat)) + datafile_list_chem.append(str(chem_file) + '_lbc' + cfg.chem_suffix) + datafile_list_chem = ' '.join([str(v) for v in datafile_list_chem]) + + # -- Copy CAM-Chem processing script in workdir + with open( + os.path.join(cfg.case_path, 'icon_camchem_lbc.sh') + ) as inf: + to_write = inf.read() + output_file = os.path.join(cfg.icon_input_icbc, 'icon_camchem_lbc.sh') + with open(output_file, "w") as outf: + outf.write(to_write.format(cfg=cfg, + datafile_list_chem = datafile_list_chem)) + + # -- Run CAM-Chem processing script + process = subprocess.Popen([ + "bash", + os.path.join(cfg.icon_input_icbc, 'icon_camchem_lbc.sh') + ], + stdout=subprocess.PIPE) + process.communicate() + + # -- Create symbolic link to LBC file at experiment start + if cfg.lrestart == '.TRUE.': + tools.symlink_file(cfg.ini_LBC_file, cfg.ini_LBC_file_scratch) + + ## ---------------------------------------------- + ## -- Merging IC/LBC for meteo and chem fields -- + ## ---------------------------------------------- + + logging.info('Merging IC and LBC') + + for time in tools.iter_hours(cfg.startdate_sim, cfg.enddate_sim, + cfg.meteo_inc): + if time == cfg.startdate_sim: + # -- Merge IC + meteo_file = os.path.join( + cfg.icon_input_icbc, + time.strftime(cfg.meteo_prefix + cfg.meteo_nameformat) + + '_ic.nc') + if os.path.isfile(meteo_file): + chem_file = os.path.join( + cfg.icon_input_icbc, + time.strftime(cfg.chem_prefix + cfg.chem_nameformat) + + '_ic_icon.nc') + merged_file = os.path.join( + cfg.icon_input_icbc, + time.strftime(cfg.inidata_prefix + cfg.inidata_nameformat) + + cfg.inidata_filename_suffix) + ds_meteo = xr.open_dataset(meteo_file) + ds_chem = xr.open_dataset(chem_file) + ds_merged = xr.merge([ds_meteo, ds_chem], + compat="override") + ds_merged.to_netcdf(merged_file) + # Rename file to get original file name + tools.remove_file(meteo_file) + tools.remove_file(chem_file) + logging.info( + "Added MECCA tracers to file {}".format(merged_file)) + + # -- Merge LBC + meteo_file = os.path.join( + cfg.icon_input_icbc, + time.strftime(cfg.meteo_prefix + cfg.meteo_nameformat) + + '_lbc.nc') + chem_file = os.path.join( + cfg.icon_input_icbc, + time.strftime(cfg.chem_prefix + cfg.chem_nameformat) + + '_lbc_icon.nc') + merged_file = os.path.join( + cfg.icon_input_icbc, + time.strftime(cfg.lbcdata_prefix + cfg.lbcdata_nameformat) + + cfg.lbcdata_filename_suffix) + ds_meteo = xr.open_dataset(meteo_file) + ds_chem = xr.open_dataset(chem_file) + ds_merged = xr.merge([ds_meteo, ds_chem], + compat="override") + ds_merged.to_netcdf(merged_file) + tools.remove_file(meteo_file) + tools.remove_file(chem_file) + logging.info( + "Added MECCA tracers to file {}".format(merged_file)) + + ## ----------------------------------- + ## -- Add Q (copy of QV) to IC file -- + ## ----------------------------------- + + logging.info('Add Q (copy of QV) to initial file') + ic_file = os.path.join( + cfg.icon_input_icbc, + cfg.startdate_sim.strftime(cfg.inidata_prefix + + cfg.inidata_nameformat + + cfg.inidata_filename_suffix) + ) + if os.path.isfile(ic_file): + merged_file = os.path.join( + cfg.icon_input_icbc, + cfg.startdate_sim.strftime(cfg.inidata_prefix + + cfg.inidata_nameformat + + '_merged' + + cfg.inidata_filename_suffix) + ) + ds = xr.open_dataset(ic_file) + merging = False + if 'Q' not in ds: + merging = True + ds['Q'] = ds['QV'] + logging.info(f"Added Q to file {ic_file}") + if merging: + ds.to_netcdf(merged_file) + tools.rename_file(merged_file, ic_file) \ No newline at end of file diff --git a/jobs/tools/camchem_ic_lbc_reggrid.py b/jobs/tools/camchem_ic_lbc_reggrid.py new file mode 100644 index 00000000..7f05aec9 --- /dev/null +++ b/jobs/tools/camchem_ic_lbc_reggrid.py @@ -0,0 +1,153 @@ +#!/usr/bin/env python3 + +import os +from os.path import join +from netCDF4 import Dataset +from .nc_operations import VariableCreator, VariableCopier, DimensionCopier +from datetime import datetime +from .. import tools + +####################################################################### +## Script for roducing IC and LBC on regular grid from CAM-Chem output. +## +## Author: Corina Keller (Empa). +####################################################################### + +def process_ic(file_path_in, file_path_out, ic_datetime, aero_dict, prefix): + """ + Copy chem and aerosol fields for specified dates (no unit transformation). + Splitting of aerosols into modes. + """ + #ic_datetime = datetime.strptime(ic_date, '%Y-%m-%d %H:%M:%S') + + in_template = prefix + '%Y%m%d%H.nc' + outname_template = prefix + '%Y%m%d%H_ic.nc' + in_filename = ic_datetime.strftime(in_template) + out_filename = ic_datetime.strftime(outname_template) + in_path = join(file_path_in, in_filename) + + with Dataset(in_path, 'r') as ds: + # -- Create output filename + out_path = join(file_path_out, out_filename) + # -- Write output netCDF file containing IC + with Dataset(out_path, 'w') as out_ds: + # -- Copy global attributes + out_ds.setncatts(ds.__dict__) + + # -- Copy dimensions + dimpairs = [(dim_name, dim_name) for dim_name, dim in ds.dimensions.items()] + dim_copier = [ + DimensionCopier(src_name, dst_name) for src_name, dst_name in dimpairs + ] + + # -- Copy variables + var_to_copy = [ + var_name for var_name, v in ds.variables.items() if var_name not in aero_dict + ] + dict_var = [{'src_names': var, + 'dst_name': var} + for var in var_to_copy] + var_copier = [VariableCopier(**kwargs) for kwargs in dict_var] + + for op in dim_copier + var_copier: + op.apply_to(ds,out_ds) + + # -- Split aerosols into modes + for var_name, modes_dict in aero_dict.items(): + for mode_name, split in modes_dict.items(): + var_mmr = ds[var_name][:] * split + (VariableCreator( + var_args={ + 'varname': mode_name, + 'datatype': 'f8', + 'dimensions': ('time', 'lev_m', 'lat', 'lon') + }, + var_attrs={ + 'long_name' : mode_name + ' mass mixing ratio', + 'units': 'kg/kg' + }, + var_vals=var_mmr).apply_to(out_ds)) + +######################################## + +def process_lbc(file_path_in, file_path_out, + start_time, end_time, inc, + prefix, suffix, nameformat, + var_mw_dict, aero_dict): + """ + Unit transformation [mol/mol] --> [kg/kg] for chem fields. + Splitting, renaming and unit transformation [kg/kg] --> [ug/kg] for aerosols. + """ + # Files in the input folder + #files = [f for f in os.listdir(file_path) if f.startswith('camchem_') and f.endswith('.nc')] + #for file in files: + for time in tools.iter_hours(start_time, end_time, + inc): + + # -- Extract datetime information from the filename + #datetime_str = file.split('_')[1].split('.')[0] + in_file = os.path.join(file_path_in, + time.strftime(prefix + nameformat + suffix)) + with Dataset(in_file, 'r') as ds: + # -- Create output filename + out_filename = time.strftime(prefix + nameformat + '_lbc' + suffix) + out_file = join(file_path_out, out_filename) + + # -- Write output netCDF file containing LBC + with Dataset(out_file, 'w') as out_ds: + ## Global attributes + out_ds.setncatts(ds.__dict__) + + # -- Copy dimensions + dimpairs = [(dim_name, dim_name) for dim_name, dim in ds.dimensions.items()] + dim_copier = [ + DimensionCopier(src_name, dst_name) for src_name, dst_name in dimpairs + ] + + # -- Copy variables + var_to_copy = ['lev_s', 'lev_m', 'lon', 'lat', 'hyam', 'hybm', 'hyai', 'hybi', 'time', 'date', 'datesec', 'PS'] + dict_var = [{'src_names': var, + 'dst_name': var} + for var in var_to_copy] + var_copier = [VariableCopier(**kwargs) for kwargs in dict_var] + + for op in dim_copier + var_copier: + op.apply_to(ds,out_ds) + + # -- Transform units [mol/mol] --> [kg/kg] of chemical fields + # Molar weight air + mw_air = 28.96 + for var_name, mw in var_mw_dict.items(): + # [mol/mol] --> [kg/kg] + var_mmr = ds[var_name][:] * mw/mw_air + (VariableCreator( + var_args={ + 'varname': var_name, + 'datatype': 'f8', + 'dimensions': ('time', 'lev_m', 'lat', 'lon') + }, + var_attrs={ + 'long_name' : var_name + ' mass mixing ratio', + 'units': 'kg/kg' + }, + var_vals=var_mmr).apply_to(out_ds)) + + # -- Split aerosols into modes and transform units [kg/kg] --> [ug/kg] + kg_to_mug = 1.e9 + for var_name, modes_dict in aero_dict.items(): + for mode_name, split in modes_dict.items(): + # [kg/kg] --> [ug/kg] + var_mmr = ds[var_name][:] * split * kg_to_mug + (VariableCreator( + var_args={ + 'varname': mode_name, + 'datatype': 'f8', + 'dimensions': ('time', 'lev_m', 'lat', 'lon') + }, + var_attrs={ + 'long_name' : mode_name + ' mass mixing ratio', + 'units': 'mug/kg' + }, + var_vals=var_mmr).apply_to(out_ds)) + +######################################## diff --git a/jobs/tools/camchem_interpolation.py b/jobs/tools/camchem_interpolation.py new file mode 100644 index 00000000..dbafbe10 --- /dev/null +++ b/jobs/tools/camchem_interpolation.py @@ -0,0 +1,537 @@ +#!/usr/bin/env python3 + +import numpy as np +import glob +import os +from netCDF4 import Dataset +from os.path import join +import multiprocessing +from .nc_operations import VariableCreator, VariableCopier, DimensionCopier +from datetime import datetime, timedelta, timezone + +########################################################################################## +## - Vertical interpolation of chemistry fields from CAM-Chem (CESM2.1) output. +## - Extraction of one NetCDF file per time step. +## - Interpolation of chemistry fields with respect to time variable. +## +## Based on 'mozart2int2m.py', with additions by Louise Aubet and Corina Keller (Empa). +########################################################################################## + +def time_intpl(file_path, prefix): + """ Interpolation of chemistry fields with respect to time variable. + """ + files = sorted([f for f in os.listdir(file_path) if f.startswith(prefix)]) + + for i in range(len(files) - 1): + + in_filename1 = join(file_path, files[i]) + in_filename2 = join(file_path, files[i + 1]) + + date_str1 = files[i].split('_')[1].split('.')[0] + date_str2 = files[i + 1].split('_')[1].split('.')[0] + + # Time difference between two files + date1 = datetime.strptime(date_str1, '%Y%m%d%H') + date2 = datetime.strptime(date_str2, '%Y%m%d%H') + time_difference = (date2 - date1).total_seconds() / 3600.0 + + # Date of the new file + delta_t = time_difference / 2. + interp_date = date1 + timedelta(hours=delta_t) + + outname_template = join(file_path, prefix + '%Y%m%d%H.nc') + out_name = interp_date.strftime(outname_template) + + with Dataset(in_filename1, 'r') as ds1, Dataset(in_filename2, 'r') as ds2: + if os.path.exists(out_name): + os.remove(out_name) + + with Dataset(out_name, 'w') as out_ds: + # -- Copy global attributes + out_ds.setncatts(ds1.__dict__) + + # -- Copy dimensions and variables + dimpairs = [(dim_name, dim_name) for dim_name, dim in ds1.dimensions.items()] + dim_copier = [ + DimensionCopier(src_name, dst_name) for src_name, dst_name in dimpairs + ] + var_to_copy = ['lev_s', 'lev_m', 'lon', 'lat', 'hyam', 'hybm', 'hyai', 'hybi', 'date'] + dict_var = [{'src_names': var, + 'dst_name': var} + for var in var_to_copy] + var_copier = [VariableCopier(**kwargs) for kwargs in dict_var] + for op in dim_copier + var_copier: + op.apply_to(ds1, out_ds) + + # -- Create time and date variables + (VariableCreator( + var_args={ + 'varname': 'time', + 'datatype': 'f8', + 'dimensions': ('time',) + }, + var_attrs={ + 'axis': 'T', + 'calendar': 'proleptic_gregorian', + 'long_name': 'simulation_time', + 'standard_name': 'time', + 'units': + interp_date.strftime('hours since %Y-%m-%d %H:%M:%S') + }, + var_vals=0).apply_to(out_ds)) + + datesec = ds1['datesec'][:] + delta_t * 3600 + (VariableCreator( + var_args={ + 'varname': 'datesec', + 'datatype': 'f8', + 'dimensions': ('time', ) + }, + var_attrs={ + 'long_name':'current seconds of current date', + 'units':'s' + }, + var_vals=datesec).apply_to(out_ds)) + + # -- Interpolate surface pressure + mean_pres = np.mean([ds1['PS'][:], ds2['PS'][:]], axis=0) + (VariableCreator( + var_args={ + 'varname': 'PS', + 'datatype': 'f8', + 'dimensions': ('time', 'lat', 'lon') + }, + var_attrs={ + 'long_name':'surface pressure', + 'units':'Pa' + }, + var_vals=mean_pres).apply_to(out_ds)) + + # -- Interpolate chemical fields + var_to_copy.extend(('time','datesec', 'PS')) + for var_name, var in ds1.variables.items(): + if var_name not in var_to_copy: + mean_field = np.mean([ds1[var_name][:], ds2[var_name][:]],axis=0) + (VariableCreator( + var_args={ + 'varname': var_name, + 'datatype': 'f8', + 'dimensions': ('time', 'lev_m', 'lat', 'lon') + }, + var_attrs=var.__dict__, + var_vals=mean_field).apply_to(out_ds)) + +################################## + +def date_from_days_since_ref(days_since_ref, ref_date): + """Convert a float counting days since a reference date to a datetime object. + + Parameters: + - days_since_ref: float, the number of days since the reference date. + - ref_date: str, the reference date in the format 'yyyy-mm-dd hh:mm:ss'. + + Returns: + - datetime object corresponding to the given number of days since the reference date. + """ + timestamp_date = datetime.strptime(ref_date, '%Y-%m-%d %H:%M:%S').date() + timestamp_time = timedelta(days=days_since_ref) + + return datetime.combine(timestamp_date, datetime.min.time()) + timestamp_time + +################################## + +def extract_timeslice(in_file, out_file_template, spec_intpl, ref_date): + """ + - Extract a time slice from the input NetCDF file. + - Adjust longitude from [0, 360) to (-180, 180]. + + Parameters: + - in_file (str): Path to the input NetCDF file containing data for each time step. + - out_file_template (str): Template for the output NetCDF file name with strftime placeholder. + - spec_intpl (list): List with species to be copied from in_file. + + Returns: + None + """ + + # -- Create list of DimensionCopier objects for copying dimensions. + dimpairs = [('lat','lat'), + ('lon','lon'), + ('lev_m','lev_m'), + ('lev_s','lev_s'), + ('nhym','nhym'), + ('nhyi','nhyi') + ] + + dim_copier = [ + DimensionCopier(src_name, dst_name) for src_name, dst_name in dimpairs + ] + + # -- Create list of VariableCopier objects for copying variables for each time step. + varpairs_to_copy = [(spc, spc) for spc in spec_intpl] + + for time_index in range(Dataset(in_file).dimensions['time'].size): + + # -- Specify the slices of the variable values to be copied. + sliced_variable_options = {'var_args': {'dimensions': ('time', + 'lev_m', + 'lat', + 'lon')}, + 'var_val_indices': np.s_[time_index, :]} + + var_opts = [{'src_names': src, + 'dst_name': dst, + **sliced_variable_options} + for src, dst in varpairs_to_copy] + + var_opts += [{'src_names': 'lat', + 'dst_name': 'lat'}, + {'src_names': 'lon', + 'dst_name': 'lon'}, + {'src_names': 'PS', + 'dst_name': 'PS', + 'var_args': {'dimensions': ('time', 'lat', 'lon')}, + 'var_val_indices': np.s_[time_index, :]}, + {'src_names': 'hyam', + 'dst_name': 'hyam', + 'var_args': {'dimensions': ('nhym', )}}, + {'src_names': 'hybm', + 'dst_name': 'hybm', + 'var_args': {'dimensions': ('nhym', )}}, + {'src_names': 'lev_m', + 'dst_name': 'lev_m'}, + {'src_names': 'lev_s', + 'dst_name': 'lev_s'}, + {'src_names': 'hyai', + 'dst_name': 'hyai', + 'var_args': {'dimensions': ('nhyi', )}}, + {'src_names': 'hybi', + 'dst_name': 'hybi', + 'var_args': {'dimensions': ('nhyi', )}}] + + var_copier = [VariableCopier(**kwargs) + for kwargs in var_opts] + + # -- Extract data for given time step. + with Dataset(in_file) as ds: + days_since_ref = float(ds.variables['time'][time_index]) + timestamp = date_from_days_since_ref(days_since_ref,ref_date) + out_path = timestamp.strftime(out_file_template) + + with Dataset(out_path, 'w') as out_ds: + # -- Copy global attributes + out_ds.setncatts(ds.__dict__) + + # -- Create new time dimension. + out_ds.createDimension(dimname='time', size=1) + + # -- Create time variable. + (VariableCreator( + var_args={ + 'varname': 'time', + 'datatype': 'f8', + 'dimensions': ('time', ) + }, + var_attrs={ + 'axis': 'T', + 'calendar': 'proleptic_gregorian', + 'long_name': 'simulation_time', + 'standard_name': 'time', + 'units': + timestamp.strftime('hours since %Y-%m-%d %H:%M:%S') + }, + var_vals=0).apply_to(out_ds)) + + # -- Create date variable. + (VariableCreator( + var_args={ + 'varname': 'date', + 'datatype': 'i4', + 'dimensions': ('time', ) + }, + var_attrs={ + 'long_name': ('current date as 8 digit' + ' integer (YYYYMMDD)') + }, + var_vals=int(timestamp.strftime('%Y%m%d'))).apply_to(out_ds)) + + # -- Create datesec variable. + sec_from_midnight = ((timestamp - timestamp.replace( + hour=0, minute=0, second=0,microsecond=0)).total_seconds()) + (VariableCreator( + var_args={ + 'varname': 'datesec', + 'datatype': 'i4', + 'dimensions': ('time', ) + }, + var_attrs={ + 'long_name': ('seconds to complete','current date'), + 'units': 's' + }, + var_vals=sec_from_midnight).apply_to(out_ds)) + + # -- Write specified dimensions and variables to output. + for op in dim_copier + var_copier: + op.apply_to(ds, out_ds) + + # -- Transform longitude: [0,360) -> [-180,180). + for i in range(len(out_ds['lon'])): + if out_ds['lon'][i] > 180: + out_ds['lon'][i] -= 360 + +################################## + +def log_interpolate_1d(dt,dlat,dlon,pres_m,pres,var_data): + + p_interp = pres_m[dt, :, dlat, dlon] + p_in = pres[dt, :, dlat, dlon] + y_in = var_data[dt, :, dlat, dlon] + + index = np.searchsorted(p_in, p_interp, side='left', sorter=np.argsort(p_in)) + y_interp = np.zeros_like(p_interp, dtype=float) + alpha = np.zeros_like(p_interp, dtype=float) + + valid_indices = (0 < index) & (index <= len(p_in) - 1) + alpha[valid_indices] = np.log(p_interp[valid_indices] / p_in[index[valid_indices] - 1]) / \ + np.log(p_in[index[valid_indices]] / p_in[index[valid_indices] - 1]) + + y_interp[valid_indices] = alpha[valid_indices] * y_in[index[valid_indices]] + \ + (1 - alpha[valid_indices]) * y_in[index[valid_indices] - 1] + + y_interp[index == 0] = y_in[0] + y_interp[index == len(p_in)] = y_in[-1] + + return dt, dlat, dlon, y_interp + +def interpolate_at_point(dt, dlat, dlon, pres_m, pres, var_data): + var_col = log_interpolate_1d(pres_m[dt, :, dlat, dlon], pres[dt, :, dlat, dlon], var_data[dt, :, dlat, dlon]) + return dt, dlat, dlon, var_col + +def hybrid_pressure_interpolation(in_ds, out_ds, var_name, time_indices): + """Perform vertical interpolation of 'var_name'. + The interpolation is linear with the log pressure. + """ + # Pressure on vertical levels of meteo data + pres_m = out_ds['pres_m'][:] + # Pressure on vertical levels of chemistry data + pres = out_ds['pres'][:] + var_data = in_ds[var_name][np.s_[time_indices,:]] + + var_arr = np.zeros((out_ds.dimensions['time'].size, + out_ds.dimensions['lev_m'].size, + out_ds.dimensions['lat'].size, + out_ds.dimensions['lon'].size)) + + args_list = [(dt, dlat, dlon, pres_m, pres, var_data) + for dt in range(len(out_ds['time'][:])) + for dlat in range(len(out_ds['lat'][:])) + for dlon in range(len(out_ds['lon'][:]))] + + num_proc = multiprocessing.cpu_count() + + with multiprocessing.Pool(num_proc) as pool: + results = pool.starmap(log_interpolate_1d, args_list) + + for result in results: + dt, dlat, dlon, var_col = result + var_arr[dt, :, dlat, dlon] = var_col + return var_arr + +def vert_intpl(chem_filename, meteo_filename, out_filename, spec, start_chunk, end_chunk, ref_date): + """Perform vertical interpolation of atmospheric chemical fields using meteorological data. + + Parameters: + - chem_filename (str): Path to the netCDF file containing atmospheric chemical data. + - meteo_filename (str): Path to the netCDF file containing meteorological data. + - out_filename (str): Path to the output netCDF file for interpolated data. + - species (list): List of chemical species names to be interpolated. + - start_chunk (datetime object): Start of simulation time. + - end_chunk (datetime object): End of simulation time. + + Returns: + None + + """ + + with Dataset(chem_filename, 'r') as c_ds, Dataset(meteo_filename, 'r') as m_ds: + + # -- Extract indices for times between simulation start and end + time_index_lst = [] + + for time_index in range(c_ds.dimensions['time'].size): + days_since_ref_now = float(c_ds.variables['time'][time_index]) + timestamp_now = date_from_days_since_ref(days_since_ref_now,ref_date) + timestamp_now = timestamp_now.replace(tzinfo=timezone.utc) + + if start_chunk <= timestamp_now <= end_chunk: + if time_index not in time_index_lst: + time_index_lst.append(time_index) + + if time_index != c_ds.dimensions['time'].size-1: + days_since_ref_next = float(c_ds.variables['time'][time_index+1]) + timestamp_next = date_from_days_since_ref(days_since_ref_next,ref_date) + timestamp_next = timestamp_next.replace(tzinfo=timezone.utc) + + # Time difference of two time steps + time_difference = (timestamp_next - timestamp_now).total_seconds() / 3600.0 + + # Date between two time steps + delta_t = time_difference / 2. + timestamp_interp = timestamp_now + timedelta(hours=delta_t) + print(time_index) + print(timestamp_interp) + + if start_chunk <= timestamp_interp <= end_chunk: + if time_index not in time_index_lst: + time_index_lst.append(time_index) + time_index_lst.append(time_index+1) + + print(time_index_lst) + + with Dataset(out_filename, 'w') as out_ds: + # -- Set global attributes + out_ds.setncatts(c_ds.__dict__) + + # -- Create new dimensions + out_ds.createDimension("lev_m", len(m_ds['lev'][:])) + out_ds.createDimension("lev_s", 1) + out_ds.createDimension('time', size=len(time_index_lst)) + + # -- Copy dimensions + dimpairs_c = [('lev', 'lev'), + ('lat', 'lat'), + ('lon', 'lon'), + ('ilev', 'ilev'), + ('nbnd', 'nbnd')] + dim_copier_c = [ + DimensionCopier(src_name, dst_name) for src_name, dst_name in dimpairs_c + ] + dimpairs_m = [('nhym', 'nhym'), + ('nhyi', 'nhyi')] + dim_copier_m = [ + DimensionCopier(src_name, dst_name) for src_name, dst_name in dimpairs_m + ] + + # -- Copy variables + copy_from_chem = ['lat', 'lon'] + dict_c = [{'src_names': var, + 'dst_name': var} + for var in copy_from_chem] + + dict_c += [{'src_names': 'PS', + 'dst_name': 'PS', + 'var_args': {'dimensions': ('time', 'lat', 'lon')}, + 'var_val_indices': np.s_[time_index_lst, :]}] + + sliced_variable_options = {'var_args': {'dimensions': ('time')}, + 'var_val_indices': np.s_[time_index_lst]} + dict_c += [{'src_names': src, + 'dst_name': dst, + **sliced_variable_options} + for src, dst in [('time','time'), + ('date', 'date'), + ('datesec','datesec')]] + + var_copier_c = [VariableCopier(**kwargs) for kwargs in dict_c] + for op in dim_copier_c + var_copier_c: + op.apply_to(c_ds, out_ds) + + copy_from_meteo = ['hyam', 'hybm', 'hyai', 'hybi'] + dict_m = [{'src_names': var, + 'dst_name': var} + for var in copy_from_meteo] + + var_copier_m = [VariableCopier(**kwargs) for kwargs in dict_m] + + for op in dim_copier_m + var_copier_m: + op.apply_to(m_ds, out_ds) + + # -- Create vertical levels and pressure arrays for vertial interpolation + + # Vertical levels from meteo data + hybrid_levels = m_ds['hyam'][:]*1e-2 + m_ds['hybm'][:]*1e3 + (VariableCreator( + var_args={ + 'varname': 'lev_m', + 'datatype': 'f8', + 'dimensions': ('lev_m', ) + }, + var_attrs={ + 'long_name':'hybrid sigma levels', + 'units':'1' + }, + var_vals=hybrid_levels).apply_to(out_ds)) + + # Surface level + (VariableCreator( + var_args={ + 'varname': 'lev_s', + 'datatype': 'f8', + 'dimensions': ('lev_s', ) + }, + var_attrs={ + 'long_name':'surface level', + 'units':'1' + }, + var_vals=1).apply_to(out_ds)) + + # Pressure array for chem data + pres_c_arr = np.zeros((len(time_index_lst), + c_ds.dimensions['lev'].size, + c_ds.dimensions['lat'].size, + c_ds.dimensions['lon'].size)) + + for ilev in range(len(c_ds['ilev'][:])-1): + pres_c_col = c_ds['hyam'][ilev]*c_ds['P0'] + c_ds['hybm'][ilev]*c_ds['PS'][np.s_[time_index_lst,:]] + pres_c_arr[:, ilev, :, :] = pres_c_col + + # Pressure on vertical levels of chem data + (VariableCreator( + var_args={ + 'varname': 'pres', + 'datatype': 'f8', + 'dimensions': ('time', 'lev', 'lat', 'lon') + }, + var_attrs={ + 'long_name':'pressure on vertical levels of chem data', + 'units':'Pa' + }, + var_vals=pres_c_arr).apply_to(out_ds)) + + # Pressure array for meteo data + pres_m_arr = np.zeros((len(time_index_lst), + m_ds.dimensions['lev'].size, + c_ds.dimensions['lat'].size, + c_ds.dimensions['lon'].size)) + + for ilev in range(len(m_ds['lev'][:])): + pres_m_col = m_ds['hyam'][ilev] + m_ds['hybm'][ilev]*c_ds['PS'][np.s_[time_index_lst,:]] + pres_m_arr[:, ilev, :, :] = pres_m_col + + # Pressure on vertical levels of meteo data + (VariableCreator( + var_args={ + 'varname': 'pres_m', + 'datatype': 'f8', + 'dimensions': ('time', 'lev_m', 'lat', 'lon') + }, + var_attrs={ + 'long_name':'pressure on vertical levels of meteo data', + 'units':'Pa' + }, + var_vals=pres_m_arr).apply_to(out_ds)) + + # -- Interpolate fields + for var_name in spec: + print('Vertical interpolation for ' + var_name) + arr = hybrid_pressure_interpolation(c_ds, out_ds, var_name, time_index_lst) + (VariableCreator( + var_args={ + 'varname': var_name, + 'datatype': 'f8', + 'dimensions': ('time', 'lev_m', 'lat', 'lon') + }, + var_attrs=c_ds[var_name].__dict__, + var_vals=arr).apply_to(out_ds)) + + diff --git a/workflows.yaml b/workflows.yaml index f87c22b5..abe778d7 100644 --- a/workflows.yaml +++ b/workflows.yaml @@ -235,3 +235,21 @@ icon-art-oem: - prepare_art_oem previous: - icon + +icon-art-full-chem: + features: + - restart + jobs: + - prepare_icon + - prepare_art_full_chem + - icon + dependencies: + prepare_art_full_chem: + current: + - prepare_icon + icon: + current: + - prepare_icon + - prepare_art_full_chem + previous: + - icon