diff --git a/cases/icon-art-full-chem/config.yaml b/cases/icon-art-full-chem/config.yaml new file mode 100644 index 00000000..b666a6eb --- /dev/null +++ b/cases/icon-art-full-chem/config.yaml @@ -0,0 +1,83 @@ +# 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 +compute_account: s1317 +ntasks_per_node: 128 +restart_step: PT03H +startdate: 2019-01-01T00:00:00Z +enddate: 2019-01-01T09:00:00Z + +eccodes_dir: /capstor/store/cscs/empa/em05/ckeller/spack-c2sm/spack/opt/spack/linux-sles15-zen2/gcc-12.3.0/eccodes-2.25.0-hrevcu72qaiav7csmdhgth323uozzutd/share/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: /capstor/store/cscs/userlab/s1317/ckeller/ICON-ART/Running_model/proc_chain_data/input/test-proc-chain-full-chem/art + +walltime: + prepare_icon: '00:20:00' + prepare_art_full_chem: '00:30:00' + icon: '01:00:00' + +meteo: + dir: /capstor/store/cscs/userlab/s1317/ckeller/ICON-ART/Running_model/proc_chain_data/input/test-proc-chain-full-chem/meteo + prefix: era5_ + nameformat: '%Y%m%d%H' + suffix: .grb + inc: 3 + +chem: + dir: /capstor/store/cscs/userlab/s1317/ckeller/ICON-ART/Running_model/proc_chain_data/input/test-proc-chain-full-chem/chem + prefix: camchem_ + nameformat: '%Y%m%d%H' + suffix: .nc + inc: 3 + icbc_filename: camchem-20190101-20190201.nc + ref_date: '2000-03-01 00:00:00' + + +input_files: + radiation_grid_filename: /capstor/store/cscs/userlab/s1317/ckeller/ICON-ART/Running_model/proc_chain_data/input/test-proc-chain-full-chem/grid/base_grid.nc + dynamics_grid_filename: /capstor/store/cscs/userlab/s1317/ckeller/ICON-ART/Running_model/proc_chain_data/input/test-proc-chain-full-chem/grid/icon_R03B07_DOM01.nc + map_file_latbc: /capstor/store/cscs/userlab/s1317/ckeller/ICON-ART/Running_model/proc_chain_data/input/test-proc-chain-full-chem/grid/map_file.latbc + lateral_boundary_grid: /capstor/store/cscs/userlab/s1317/ckeller/ICON-ART/Running_model/proc_chain_data/input/test-proc-chain-full-chem/grid/lateral_boundary_DOM01.grid.nc + extpar_filename: /capstor/store/cscs/userlab/s1317/ckeller/ICON-ART/Running_model/proc_chain_data/input/test-proc-chain-full-chem/grid/extpar_DOM01.nc + cldopt_filename: /capstor/store/cscs/userlab/s1317/ckeller/ICON-ART/Running_model/proc_chain_data/input/test-proc-chain-full-chem/rad/rrtm_cldopt.nc + lrtm_filename: /capstor/store/cscs/userlab/s1317/ckeller/ICON-ART/Running_model/proc_chain_data/input/test-proc-chain-full-chem/rad/rrtmg_lw.nc + map_file_ana: /capstor/store/cscs/userlab/s1317/ckeller/ICON-ART/Running_model/proc_chain_data/input/test-proc-chain-full-chem/mapping/map_file.ana + meccatracer_xml_filename: /capstor/store/cscs/userlab/s1317/ckeller/ICON-ART/Running_model/proc_chain_data/input/test-proc-chain-full-chem/xml/mecca_tracers.xml + oem_vertical_profiles_nc: /capstor/store/cscs/userlab/s1317/ckeller/ICON-ART/Running_model/proc_chain_data/input/test-proc-chain-full-chem/oem/vertical_profiles.nc + oem_gridded_emissions_nc: /capstor/store/cscs/userlab/s1317/ckeller/ICON-ART/Running_model/proc_chain_data/input/test-proc-chain-full-chem/oem/oem_gridded_emissions.nc + oem_hourofday_nc: /capstor/store/cscs/userlab/s1317/ckeller/ICON-ART/Running_model/proc_chain_data/input/test-proc-chain-full-chem/oem/hourofday.nc + oem_dayofweek_nc: /capstor/store/cscs/userlab/s1317/ckeller/ICON-ART/Running_model/proc_chain_data/input/test-proc-chain-full-chem/oem/dayofweek.nc + oem_monthofyear_nc: /capstor/store/cscs/userlab/s1317/ckeller/ICON-ART/Running_model/proc_chain_data/input/test-proc-chain-full-chem/oem/monthofyear.nc + aerodyn_tracers: /capstor/store/cscs/userlab/s1317/ckeller/ICON-ART/Running_model/proc_chain_data/input/test-proc-chain-full-chem/aerodyn/tracers_aerosol.xml + aerodyn_modes: /capstor/store/cscs/userlab/s1317/ckeller/ICON-ART/Running_model/proc_chain_data/input/test-proc-chain-full-chem/aerodyn/modes.xml + aerodyn_coag: /capstor/store/cscs/userlab/s1317/ckeller/ICON-ART/Running_model/proc_chain_data/input/test-proc-chain-full-chem/aerodyn/coagulate.xml + aerodyn_emiss: /capstor/store/cscs/userlab/s1317/ckeller/ICON-ART/Running_model/proc_chain_data/input/test-proc-chain-full-chem/aerodyn/emiss.xml + aerodyn_diag: /capstor/store/cscs/userlab/s1317/ckeller/ICON-ART/Running_model/proc_chain_data/input/test-proc-chain-full-chem/aerodyn/diag.xml + +icon: + binary_file: /capstor/store/cscs/userlab/s1317/ckeller/ICON-ART/Running_model/proc_chain_data/input/test-proc-chain-full-chem/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: '01:00:00' + np_tot: 8 + 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..b8686c7f --- /dev/null +++ b/cases/icon-art-full-chem/icon_runjob.cfg @@ -0,0 +1,383 @@ +#!/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 + cart_diagnostics_xml = '{cfg.input_files_scratch_aerodyn_diag}' + 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 = 95 ! 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 = 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 = 480. ! time step for radiation in s + dt_conv = 120. ! time step for convection in s (domain specific) + dt_sso = 240. ! time step for SSO parameterization + dt_gwd = 240. ! 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., 100000000., 3600. ! start, end, increment + steps_per_file = 1 ! number of steps per file + include_last = .TRUE. + steps_per_file_inclfirst = .FALSE. + remap = 0 + output_filename = 'icon-art-full-chem' + filename_format = '{cfg.icon_output}/_' ! 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', + 'group:ATMO_ML_VARS' , + 'group:PBL_VARS', + 'group:PRECIP_VARS', + 'group:NH_PROG_VARS', +/ + +! 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/config.py b/config.py index 8540c43d..8e1a445f 100644 --- a/config.py +++ b/config.py @@ -71,6 +71,8 @@ def __init__(self, casename): # Specific settings based on the node type ('gpu' or 'mc') if self.machine == 'daint': self.set_node_info() + elif self.machine == 'eiger': + self.set_node_info() def load_config_file(self): """Load configuration settings from a YAML file and set them as attributes. @@ -152,6 +154,8 @@ def set_machine(self): self.machine = 'daint' elif hostname.startswith('eu-'): self.machine = 'euler' + elif hostname.startswith('eiger'): + self.machine = 'eiger' else: raise ValueError(f"Unsupported hostname: {hostname}") print(f"You are on the {self.machine} machine.") @@ -191,7 +195,7 @@ def set_node_info(self): 'export MPICH_G2G_PIPELINE=256\n' 'export CRAY_CUDA_MPS=1\n') elif self.constraint == 'mc': - self.ntasks_per_node = 36 + self.ntasks_per_node = 128 self.mpich_cuda = '' else: raise ValueError( @@ -472,6 +476,22 @@ def submit_basic_python(self, job_name): f'./run_chain.py {self.casename} -j {job_name} -c {self.chunk_id} -f -s --no-logging', '', ] + elif self.machine == 'eiger': + script_lines = [ + '#!/usr/bin/env bash', + f'#SBATCH --job-name={job_name}', + '#SBATCH --nodes=1', + f'#SBATCH --time={walltime}', + f'#SBATCH --output={self.logfile}', + '#SBATCH --open-mode=append', + f'#SBATCH --account={self.compute_account}', + f'#SBATCH --partition={self.compute_queue}', + f'#SBATCH --constraint={self.constraint}', + '', + f'cd {self.chain_src_dir}', + f'./run_chain.py {self.casename} -j {job_name} -c {self.chunk_id} -f -s --no-logging', + '', + ] job_path = self.chain_root / 'job_scripts' job_path.mkdir(parents=True, exist_ok=True) @@ -516,6 +536,17 @@ def wait_for_previous(self): f'#SBATCH --dependency=afterany:{dep_str}', '', '# Do nothing', 'exit 0' ] + elif self.machine == 'eiger': + script_lines = [ + '#!/usr/bin/env bash', '#SBATCH --job-name="wait"', + '#SBATCH --nodes=1', '#SBATCH --time=00:01:00', + f'#SBATCH --output={log_file}', + f'#SBATCH --account={self.compute_account}', + f'#SBATCH --partition={self.compute_queue}', + f'#SBATCH --constraint={self.constraint}', + f'#SBATCH --dependency=afterany:{dep_str}', '', + '# Do nothing', 'exit 0' + ] with open(job_file, mode='w') as wait_job: wait_job.write('\n'.join(script_lines)) diff --git a/docs/config.rst b/docs/config.rst index f38ef329..c821eedb 100644 --- a/docs/config.rst +++ b/docs/config.rst @@ -154,13 +154,13 @@ files within the test cases. +-----------------------+-------------------------------------------------------------------------------------------------------------------------------------+ | Dictionary variable | Used in job | +=======================+=====================================================================================================================================+ -| ``meteo`` | ``prepare_cosmo``, ``prepare_icon``, ``icontools``, ``int2lm``, ``icon`` | +| ``meteo`` | ``prepare_cosmo``, ``prepare_icon``, ``prepare_art_full_chem``, ``icontools``, ``int2lm``, ``icon`` | +-----------------------+-------------------------------------------------------------------------------------------------------------------------------------+ | ``icontools_runjobs`` | ``icontools`` | +-----------------------+-------------------------------------------------------------------------------------------------------------------------------------+ | ``input_files`` | ``prepare_icon`` | +-----------------------+-------------------------------------------------------------------------------------------------------------------------------------+ -| ``chem`` | ``prepare_icon`` | +| ``chem`` | ``prepare_icon``, ``prepare_art_full_chem`` | +-----------------------+-------------------------------------------------------------------------------------------------------------------------------------+ | ``era5`` | ``prepare_icon`` | +-----------------------+-------------------------------------------------------------------------------------------------------------------------------------+ @@ -186,5 +186,5 @@ files within the test cases. +-----------------------+-------------------------------------------------------------------------------------------------------------------------------------+ | ``verify_chain`` | ``verify_chain`` | +-----------------------+-------------------------------------------------------------------------------------------------------------------------------------+ -| ``icon`` | ``oem``, ``prepare_icon``, ``icon`` | +| ``icon`` | ``oem``, ``prepare_icon``, ``prepare_art_full_chem``, ``icon`` | +-----------------------+-------------------------------------------------------------------------------------------------------------------------------------+ diff --git a/docs/functions.rst b/docs/functions.rst index 4305a932..c36b9d86 100644 --- a/docs/functions.rst +++ b/docs/functions.rst @@ -21,6 +21,7 @@ Jobs * :func:`jobs.prepare_icon.main` * :func:`jobs.reduce_output.main` * :func:`jobs.verify_chain.main` +* :func:`jobs.prepare_art_full_chem.main` ------------------------------------------- @@ -94,6 +95,10 @@ Jobs .. autofunction:: jobs.verify_chain.main +------------------------------------------- + +.. autofunction:: jobs.prepare_art_full_chem.main + Tools ----- @@ -112,6 +117,12 @@ a look into ``jobs/tools`` directly. * :func:`jobs.tools.vprmsplit.main` * :func:`jobs.tools.write_cosmo_input_ghg.main` * :func:`jobs.tools.write_int2lm_input_art.main` +* :func:`jobs.tools.camchem_ic_lbc_reggrid.process_ic` +* :func:`jobs.tools.camchem_ic_lbc_reggrid.process_lbc` +* :func:`jobs.tools.camchem_interpolation.time_intpl` +* :func:`jobs.tools.camchem_interpolation.vert_intpl` +* :func:`jobs.tools.camchem_interpolation.extract_timeslice` + ------------------------------------------- @@ -153,3 +164,23 @@ a look into ``jobs/tools`` directly. ------------------------------------------- .. autofunction:: jobs.tools.write_int2lm_input_art.main + +------------------------------------------- + +.. autofunction:: jobs.tools.camchem_ic_lbc_reggrid.process_ic + +------------------------------------------- + +.. autofunction:: jobs.tools.camchem_ic_lbc_reggrid.process_lbc + +------------------------------------------- + +.. autofunction:: jobs.tools.camchem_interpolation.time_intpl + +------------------------------------------- + +.. autofunction:: jobs.tools.camchem_interpolation.vert_intpl + +------------------------------------------- + +.. autofunction:: jobs.tools.camchem_interpolation.extract_timeslice diff --git a/docs/requirements.txt b/docs/requirements.txt index 96e1f7f9..1608630e 100644 --- a/docs/requirements.txt +++ b/docs/requirements.txt @@ -13,3 +13,4 @@ scikit-learn sphinx sphinx_rtd_theme sphinx-copybutton +setuptools \ No newline at end of file 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/icon.py b/jobs/icon.py index 250df513..7ed50cf9 100644 --- a/jobs/icon.py +++ b/jobs/icon.py @@ -44,6 +44,12 @@ def main(cfg): # Get name of initial file if hasattr(cfg, 'inicond_filename'): inidata_filename = cfg.icon_input_icbc / cfg.inicond_filename + elif (hasattr(cfg, 'inidata_prefix') and hasattr(cfg, 'inidata_nameformat') + and hasattr(cfg, 'inidata_filename_suffix')): + inidata_filename = cfg.icon_input_icbc / str( + cfg.startdate.strftime(cfg.inidata_prefix + + cfg.inidata_nameformat + + cfg.inidata_filename_suffix)) else: inidata_filename = cfg.icon_input_icbc / str( cfg.startdate_sim.strftime(cfg.meteo['prefix'] + diff --git a/jobs/prepare_art_full_chem.py b/jobs/prepare_art_full_chem.py new file mode 100644 index 00000000..57dcf841 --- /dev/null +++ b/jobs/prepare_art_full_chem.py @@ -0,0 +1,461 @@ +#!/usr/bin/env python +# -*- coding: utf-8 -*- + +import os +import logging +import numpy as np +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 all meteorological and chemistry initial and boundary conditions + (IC/BC) required for an ICON-ART full-chemistry simulation. + + The workflow includes: + - Initializing configuration and logging. + - Generating meteorological IC/LBC from ERA5 using template bash scripts. + - Processing CAM-Chem fields (vertical interpolation, temporal slicing, + temporal interpolation). + - Generating chemical IC/LBC from CAM-Chem using template bash scripts. + - Merging meteorological and chemistry fields into final ICON IC/LBC files. + - Ensuring required variables (e.g., PS, Q) are present and consistent. + + Note + ---- + Chemical tracers, molar weights, aerosol species, and aerosol–mode partitioning + are **hardcoded** in this script. + To add/remove species or alter aerosol mode partitioning, this script must + be modified accordingly. + + Parameters + ---------- + cfg : Config + Object holding all user-configuration parameters as attributes. + + Returns + ------- + None + All outputs are written as NetCDF files in the configured IC/BC directory. + """ + 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 + MW_NH4 = 18.04 + + aero_mode_dict = { + 'pom_a4': { + 'pom_mixed_ait': 0.9, # 90%-->mixed, 10%-->sol/insol + 'pom_insol_ait': 0.1 + }, + 'pom_a1': { + 'pom_mixed_acc': 0.9, + 'pom_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, + 'nh4_mixed_ait': 0.9 * MW_NH4 / MW_NH4HSO4, + 'nh4_sol_ait': 0.1 * MW_NH4 / MW_NH4HSO4 + }, + 'so4_a1': { + 'so4_mixed_acc': 0.9 * MW_SO4 / MW_NH4HSO4, + 'so4_sol_acc': 0.1 * MW_SO4 / MW_NH4HSO4, + 'nh4_mixed_acc': 0.9 * MW_NH4 / MW_NH4HSO4, + 'nh4_sol_acc': 0.1 * MW_NH4 / MW_NH4HSO4 + }, + 'so4_a3': { + 'so4_mixed_coa': 0.9 * MW_SO4 / MW_NH4HSO4, + 'so4_sol_coa': 0.1 * MW_SO4 / MW_NH4HSO4, + 'nh4_mixed_coa': 0.9 * MW_NH4 / MW_NH4HSO4, + 'nh4_sol_coa': 0.1 * MW_NH4 / MW_NH4HSO4 + }, + '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 + }, + } + + # -- 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) + # Replace "PS" from CAM-Chem by ERA5 value + ds_chem["PS"] = np.exp(ds_meteo["LNPS"]) + ds_chem["PS"] = ds_chem["PS"].squeeze(dim="lev_2") + #if 'Q' not in ds_chem: + ds_chem['Q'] = ds_meteo['QV'] + logging.info(f"Added PS and Q to file {merged_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) + + if 'GEOP_ML' not in ds_meteo.variables: + logging.warning( + f"'GEOP_ML' missing in {meteo_file}. Attempting to copy from 00:00 UTC file." + ) + zero_hour = time.replace(hour=0, minute=0, second=0) + zero_hour_file = os.path.join( + cfg.icon_input_icbc, + zero_hour.strftime(cfg.meteo_prefix + cfg.meteo_nameformat) + + '_lbc.nc') + + try: + ds_zero_hour = xr.open_dataset(zero_hour_file) + if 'GEOP_ML' in ds_zero_hour.variables: + ds_meteo['GEOP_ML'] = ds_zero_hour['GEOP_ML'] + logging.info( + f"Copied 'GEOP_ML' from {zero_hour_file} to {meteo_file}." + ) + else: + logging.error(f"'GEOP_ML' not found in {zero_hour_file}.") + except FileNotFoundError: + logging.error(f"00:00 UTC file {zero_hour_file} not found.") + except Exception as e: + logging.error( + f"Failed to process 00:00 UTC file {zero_hour_file}: {e}") + + 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) and PS=(exp(LNPS)) to IC file -- + ## ------------------------------------------------------ + + logging.info('Add Q (copy of QV) and PS (exp(LNPS)) 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 "PS" not in ds: + merging = True + ds["PS"] = np.exp(ds["LNPS"]) + ds_chem["PS"] = ds_chem["PS"].squeeze(dim="lev_2") + for var in ds.data_vars: + ds[var].encoding = {} + logging.info(f"Added PS to file {ic_file}") + logging.info(f"Added PS to file {ic_file}") + if 'Q' not in ds: + merging = True + ds['Q'] = ds['QV'] + # -- Bug fix for: + # "Variable has conflicting _FillValue (nan) and missing_value (-9e+33)" + for var in ds.data_vars: + ds[var].encoding = {} + logging.info(f"Added Q to file {ic_file}") + # add surface pressure compatible with ERA5 data + if merging: + ds.to_netcdf(merged_file) + tools.rename_file(merged_file, ic_file) diff --git a/jobs/prepare_icon.py b/jobs/prepare_icon.py index 9c0ae7ec..41fc4d25 100644 --- a/jobs/prepare_icon.py +++ b/jobs/prepare_icon.py @@ -34,6 +34,11 @@ def set_cfg_variables(cfg): cfg.restart_filename = 'restart_atm_DOM01.nc' cfg.restart_file = cfg.icon_restart_in / cfg.restart_filename cfg.restart_file_scratch = cfg.icon_work / cfg.restart_filename + cfg.ini_LBC_filename = cfg.startdate.strftime( + cfg.lbcdata_prefix + cfg.lbcdata_nameformat + + cfg.lbcdata_filename_suffix) + cfg.ini_LBC_file = cfg.icon_input_icbc_prev / cfg.ini_LBC_filename + cfg.ini_LBC_file_scratch = cfg.icon_input_icbc / cfg.ini_LBC_filename # Nudge type (global or nothing) cfg.nudge_type = 2 if hasattr(cfg, diff --git a/jobs/tools/camchem_ic_lbc_reggrid.py b/jobs/tools/camchem_ic_lbc_reggrid.py new file mode 100644 index 00000000..085dab75 --- /dev/null +++ b/jobs/tools/camchem_ic_lbc_reggrid.py @@ -0,0 +1,233 @@ +#!/usr/bin/env python3 + +import os +from os.path import join +from netCDF4 import Dataset +from .nc_operations import VariableCreator, VariableCopier, DimensionCopier +from .. import tools + + +def process_ic(file_path_in, file_path_out, ic_datetime, aero_dict, prefix): + """ + Generate initial condition (IC) NetCDF files from CAM-Chem output. + + This function copies chemical and aerosol fields for a specified datetime + from an input CAM-Chem NetCDF file. Aerosol fields can be split into + different modes. + No unit transformation is performed. + + Parameters + ---------- + file_path_in : str + Path to the input directory containing CAM-Chem NetCDF files. + file_path_out : str + Path to the output directory where IC files will be saved. + ic_datetime : datetime.datetime + Datetime for the IC to process. + aero_dict : dict + Dictionary specifying how aerosol variables should be split into modes + and their corresponding fractions. + Format: `{aerosol_var_name: {mode_name: split_fraction, ...}, ...}`. + prefix : str + Filename prefix for input and output NetCDF files. + + Returns + ------- + None + + Notes + ----- + - Input files are expected to follow the naming convention: `{prefix}%Y%m%d%H.nc`. + - Output files are saved as `{prefix}%Y%m%d%H_ic.nc`. + - Variables not present in `aero_dict` are copied directly. + - Aerosol variables are split according to the fractions defined in `aero_dict`. + """ + + 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): + """ + Generate lateral boundary condition (LBC) NetCDF files from CAM-Chem output. + + This function processes CAM-Chem NetCDF files and produces LBC files with + appropriate unit transformations. + Chemical fields are converted from [mol/mol] to [kg/kg], while aerosol fields + are split into modes and converted from [kg/kg] to [μg/kg]. + + Parameters + ---------- + file_path_in : str + Path to the input directory containing CAM-Chem NetCDF files. + file_path_out : str + Path to the output directory where LBC files will be saved. + start_time : datetime.datetime + Start datetime for the LBC processing range. + end_time : datetime.datetime + End datetime for the LBC processing range. + inc : int + Time increment in hours between consecutive files to process. + prefix : str + Filename prefix used for both input and output files. + suffix : str + Filename suffix used for both input and output files (e.g., '.nc'). + nameformat : str + Datetime format string used in the filenames (e.g., '%Y%m%d%H'). + var_mw_dict : dict + Dictionary mapping chemical variable names to their molecular weights + for unit conversion. Format: `{var_name: molecular_weight, ...}`. + aero_dict : dict + Dictionary specifying how aerosol variables should be split into modes + and their corresponding fractions. + Format: `{aerosol_var_name: {mode_name: split_fraction, ...}, ...}`. + + Returns + ------- + None + + Notes + ----- + - Input files are iterated over using specified increments between `start_time` + and `end_time`. + - Chemical fields are converted from [mol/mol] to [kg/kg] using the molecular + weight of air (28.96 g/mol) and the variable-specific molecular weight. + - Aerosol fields are split into modes and converted from [kg/kg] to [μg/kg]. + - Essential variables like coordinates and pressure are copied directly + - Output files are saved with `_lbc` appended to the original filename. + """ + + for time in tools.iter_hours(start_time, end_time, inc): + + # -- Extract datetime information from the filename + 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..12def493 --- /dev/null +++ b/jobs/tools/camchem_interpolation.py @@ -0,0 +1,680 @@ +#!/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): + """ + Interpolate chemistry fields with respect to time between consecutive NetCDF files. + + This function generates new NetCDF files at intermediate times between existing + CAM-Chem output files. Takes the arithmetic mean of chemical and surface pressure + fields and creates interpolated 'time' and 'datesec' variables. + + Parameters + ---------- + file_path : str + Path to the directory containing the input NetCDF files. + prefix : str + Filename prefix for identifying the input NetCDF files to interpolate. + + Returns + ------- + None + + """ + + 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 single time slice from a NetCDF file and adjust longitudes. + + Copies selected species and metadata from an input NetCDF file into a new + NetCDF file for each time step. + Longitudes are adjusted from [0, 360) to [-180, 180). + + Parameters + ---------- + in_file : str + Path to the input NetCDF file. + out_file_template : str + Template string for the output filename. + spec_intpl : list of str + List of chemical species to extract from the input file. + ref_date : str + Reference date (YYYY-MM-DD HH:MM:SS) used to convert time variables. + + 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_intpl_timeslice(pres_m, pres, var_data): + """ + Perform vertical interpolation of chemical data using log-pressure weighting. + + For each chemistry level, the two nearest meteorological levels are identified + and a logarithmic interpolation is applied. + + Parameters + ---------- + pres_m : ndarray + Pressure on meteorological levels (Pa). + pres : ndarray + Pressure on chemical levels (Pa). + var_data : ndarray + Chemical field data to interpolate. + + Returns + ------- + ndarray + Interpolated chemical field on meteorological levels. + """ + pres_m_flipped = np.flip(pres_m, axis=0) + pres_flipped = np.flip(pres, axis=0) + var_data_flipped = np.flip(var_data, axis=0) + + # -- Find meteo levels lying above given chem level + # -- (that is, meteo levels with lower pressure than given chem level) + above_target = pres_flipped[:, np.newaxis] <= pres_m_flipped + + # -- Calculate number of levels below given chem level + # -- (that is, number of meteo levels with higher pressire than given chem level) + lev_count = (pres_flipped[:, np.newaxis] >= pres_m_flipped).sum(axis=0) + + # -- Find indices of first meteo levels lying above given chem level + # -- The index is set to -1 if no meteo levels are found + first_above_target = np.where(lev_count == len(pres_flipped), -1, + np.argmax(above_target, axis=0)) + + # -- Select indices of 2 closest neighbouring levels for given chem level + vertical_indices = np.stack([first_above_target, first_above_target - 1], + axis=0) + + # -- Set both indices to highest/lower chem level if given chem level is above/below all meteo levels + vertical_indices[:, first_above_target == 0] = 0 + vertical_indices[:, first_above_target == -1] = len(pres_flipped) - 1 + + # -- Select pressure values of 2 closest neighbouring levels + pk_below = np.take_along_axis( + (pres_flipped), vertical_indices[1], + axis=0) # this is the pressure of the lower level (higher pressure) + pk_up = np.take_along_axis( + (pres_flipped), vertical_indices[0], + axis=0) # this is the pressure of the upper level (lower pressure) + + # -- Compute interpolation weights + alpha = (np.log(pres_m_flipped / pk_up)) / (np.log(pk_below / pk_up)) + weights = np.stack([alpha, 1 - alpha], axis=0) + + # -- Weights for chem levels lying above/below all meteo levels + weights[:, first_above_target == 0] = 0.5 + weights[:, first_above_target == -1] = 0.5 + + # Interpolate + var_data_intpl_flipped = np.take_along_axis(var_data_flipped, vertical_indices[1], axis=0) * weights[0] + \ + np.take_along_axis(var_data_flipped, vertical_indices[0], axis=0) * weights[1] + + var_data_intpl = np.flip(var_data_intpl_flipped, axis=0) + return var_data_intpl + + +def hybrid_pressure_interpolation(in_ds, out_ds, var_name, time_indices): + """ + Interpolate chemical variable vertically using log-pressure interpolation. + + Parameters + ---------- + in_ds : netCDF4.Dataset + Input dataset containing the chemical variable. + out_ds : netCDF4.Dataset + Output dataset with meteorological levels. + var_name : str + Name of the chemical variable to interpolate. + time_indices : list of int + Time indices to process in the datasets. + + Returns + ------- + ndarray + Array of vertically interpolated variable on meteorological levels. + """ + # 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)) + + for dt in range(len(out_ds['time'][:])): + var_arr[dt, :, :, :] = log_intpl_timeslice(pres_m[dt, :, :, :], + pres[dt, :, :, :], + var_data[dt, :, :, :]) + + return var_arr + + +def vert_intpl(chem_filename, meteo_filename, out_filename, spec, start_chunk, + end_chunk, ref_date): + """ + Perform vertical interpolation of chemical fields using meteorological data. + + This function interpolates chemical species from their CAM-Chem vertical levels + onto the hybrid sigma-pressure levels of the meteorological dataset. + + Parameters + ---------- + chem_filename : str + Path to the NetCDF file containing chemical data. + meteo_filename : str + Path to the NetCDF file containing meteorological data. + out_filename : str + Path for the output NetCDF file with interpolated results. + spec : list of str + List of chemical species to interpolate. + start_chunk : datetime.datetime + Start of the simulation time window to process. + end_chunk : datetime.datetime + End of the simulation time window to process. + ref_date : str + Reference date (YYYY-MM-DD HH:MM:SS) used to convert time variables. + + 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) + + 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) + + 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: + 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