diff --git a/run_scripts/extpar_mch.copernicus.sh b/run_scripts/extpar_mch.copernicus.sh new file mode 100755 index 00000000..ff4a66fd --- /dev/null +++ b/run_scripts/extpar_mch.copernicus.sh @@ -0,0 +1,666 @@ +#!/bin/bash + +#SBATCH --job-name="extpar" +#SBATCH --nodes=1 +#SBATCH --output="job.out" +#SBATCH --time=00:23:00 +#SBATCH --partition=postproc +#SBATCH --cpus-per-task=64 + + +#-------------------------------------------------------------------------------- +# variables to define by user +#-------------------------------------------------------------------------------- + +# define model for which Extpar should run: c1, c2, i05, i1, i2, i1_dev, i2_dev +model="i05" + +# Sandbox (make sure you have enough disk place at that location)! +sandboxdir=$SCRATCH/ExtPar/output/${model} + +#-------------------------------------------------------------------------------- +# environment +#-------------------------------------------------------------------------------- + +source ../modules.env +source ../.venv/bin/activate +export OMP_NUM_THREADS=$SLURM_CPUS_PER_TASK # manually set to 1 if run as ./script.sh + +export NETCDF_OUTPUT_FILETYPE=NETCDF4 + +# import functions to launch Extpar executables +. ../test/testsuite/bin/runcontrol_functions.sh + +ulimit -s unlimited +ulimit -c unlimited + +# get hostname +hostname="`echo $HOSTNAME`" +logfile="extpar_runscript.log" + + +# auto-set paths + +# directory of runscripts => need to be as in original repository +scriptdir=`pwd` +src_python=${scriptdir}/../python/lib + +# change dir to src_python to get absolute path +cd $src_python +export PYTHONPATH=$PYTHONPATH:$(pwd) +cd - > /dev/null 2>&1 + +# directory of compiled extpar executables +exedir=$scriptdir/../bin + +# define host-dependent paths and variables +# CSCS-machines +#if [[ $hostname == tsa* || $hostname == arolla* || $hostname == daint*]]; then +if [[ $hostname == tsa* || $hostname == arolla* ]]; then + + # NetCDF raw data for external parameter + data_dir=/store/c2sm/extpar_raw_data/linked_data + +# balfrin and tasna +elif [[ $hostname == nid* ]]; then + + # NetCDF raw data for external parameter + #data_dir=/store_new/mch/c2sm/extpar_raw_data/linked_data + data_dir=/scratch/mch/csteger/ExtPar/linked_data + +# unkown host +else + + # exit script in case of unknown host + echo ERROR: Unkown host: $hostname >> ${logfile} + exit 1 +fi + +#-------------------------------------------------------------------------------- +# define model (COSMO-1, COSMO-2, ICON-1, ICON-2) dependent variables +#-------------------------------------------------------------------------------- + +if [[ $model == "c1" ]]; then + + #output file names + netcdf_output_filename='external_parameter_mch_cosmo1.nc' + + # grid definition + rot_pol_lon=-170.0 + rot_pol_lat=43.0 + startlon_tot=-9.0 + startlat_tot=-9.0 + dlon=0.01 + dlat=0.01 + ie_tot=1801 + je_tot=1801 + + lsso_param=".TRUE." + lsubtract_mean_slope=".FALSE." + + # orography raw data + ntiles_column=2 + ntiles_row=4 + topo_files="'ASTER_orig_T006.nc' 'ASTER_orig_T007.nc' 'ASTER_orig_T018.nc' 'ASTER_orig_T019.nc' 'ASTER_orig_T030.nc' 'ASTER_orig_T031.nc' 'ASTER_orig_T042.nc' 'ASTER_orig_T043.nc'" + + #orography smoothing + lsmooth_oro=".TRUE." + ilow_pass_oro=1 + numfilt_oro=2 + eps_filter=1.7 + + # soil: tiles + itile_mode=0 + + # model grid + model_grid_type=2 + name_model_grid="INPUT_COSMO_GRID" + +elif [[ $model == "c2" ]]; then + + #output file names + netcdf_output_filename='external_parameter_mch_cosmo2.nc' + + # grid definition + rot_pol_lon=-170.0 + rot_pol_lat=43.0 + startlon_tot=-9.0 + startlat_tot=-9.0 + dlon=0.02 + dlat=0.02 + ie_tot=901 + je_tot=901 + + # orography raw data + lsso_param='.TRUE.' + ntiles_column=2 + ntiles_row=4 + topo_files="'ASTER_orig_T006.nc' 'ASTER_orig_T007.nc' 'ASTER_orig_T018.nc' 'ASTER_orig_T019.nc' 'ASTER_orig_T030.nc' 'ASTER_orig_T031.nc' 'ASTER_orig_T042.nc' 'ASTER_orig_T043.nc'" + + #orography smoothing + lsmooth_oro='.TRUE.' + ilow_pass_oro=1 + numfilt_oro=2 + eps_filter=1.7 + + # soil: tiles + itile_mode=0 + + # model grid + model_grid_type=2 + name_model_grid='INPUT_COSMO_GRID' + +elif [[ $model == "i1_dev" ]]; then + + #output file names + netcdf_output_filename="external_parameter_mch_ICON_1E_R19B08_DOM1.nc" + + # grid definition + grid_dir="/store/s83/tsm/ICON_INPUT/icon-1e_dev/" + grid_nc="ICON-1E_DOM01.nc" + + lsso_param='.TRUE.' + lsubtract_mean_slope='.FALSE.' + + # orography raw data + ntiles_column=2 + ntiles_row=4 + topo_files="'ASTER_orig_T006.nc' 'ASTER_orig_T007.nc' 'ASTER_orig_T018.nc' 'ASTER_orig_T019.nc' 'ASTER_orig_T030.nc' 'ASTER_orig_T031.nc' 'ASTER_orig_T042.nc' 'ASTER_orig_T043.nc'" + + #orography smoothing + lsmooth_oro='.FALSE.' + ilow_pass_oro=1 + numfilt_oro=2 + eps_filter=1.7 + + # soil: tiles + itile_mode=1 + + # model grid type + model_grid_type=1 + name_model_grid='INPUT_ICON_GRID' + +elif [[ $model == "i2_dev" ]]; then + + #output file names + netcdf_output_filename="external_parameter_mch_ICON_2E_R19B07_DOM1.nc" + + # grid definition + grid_dir="/store/s83/tsm/ICON_INPUT/icon-2e_dev/" + grid_nc="ICON-2E_DOM01.nc" + + lsso_param='.TRUE.' + lsubtract_mean_slope='.FALSE.' + + # orography raw data + ntiles_column=2 + ntiles_row=4 + topo_files="'ASTER_orig_T006.nc' 'ASTER_orig_T007.nc' 'ASTER_orig_T018.nc' 'ASTER_orig_T019.nc' 'ASTER_orig_T030.nc' 'ASTER_orig_T031.nc' 'ASTER_orig_T042.nc' 'ASTER_orig_T043.nc'" + + #orography smoothing + lsmooth_oro='.FALSE.' + ilow_pass_oro=1 + numfilt_oro=2 + eps_filter=1.7 + + # soil: tiles + itile_mode=1 + + # model grid type + model_grid_type=1 + name_model_grid='INPUT_ICON_GRID' + +elif [[ $model == "i05" ]]; then + + #output file names + netcdf_output_filename="extpar_icon_grid_00005_R19B09_DOM02.nc" + + # grid definition + grid_dir="/store_new/mch/msopr/glori/glori-ch500-nested/grid/" + grid_nc="icon_grid_00005_R19B09_DOM02.nc" + + lsso_param='.TRUE.' + lsubtract_mean_slope='.FALSE.' + + # orography raw data + ntiles_column=2 + ntiles_row=2 + topo_files="'Copernicus_DEM_N60-N50_W020-E000.nc' 'Copernicus_DEM_N60-N50_E000-E020.nc' 'Copernicus_DEM_N50-N40_W020-E000.nc' 'Copernicus_DEM_N50-N40_E000-E020.nc'" + + #orography smoothing + lsmooth_oro='.FALSE.' + ilow_pass_oro=1 + numfilt_oro=2 + eps_filter=1.7 + + # soil: tiles + itile_mode=1 + + # model grid type + model_grid_type=1 + name_model_grid='INPUT_ICON_GRID' + +elif [[ $model == "i1" ]]; then + + #output file names + netcdf_output_filename="external_parameter_icon_grid_0001_R19B08_mch.nc" + + # grid definition + grid_dir="/oprusers/osm/opr/data/grid_descriptions/" + grid_nc="icon_grid_0001_R19B08_mch.nc" + + lsso_param='.TRUE.' + lsubtract_mean_slope='.FALSE.' + + # orography raw data + ntiles_column=2 + ntiles_row=2 + topo_files="'Copernicus_DEM_N60-N50_W020-E000.nc' 'Copernicus_DEM_N60-N50_E000-E020.nc' 'Copernicus_DEM_N50-N40_W020-E000.nc' 'Copernicus_DEM_N50-N40_E000-E020.nc'" + + #orography smoothing + lsmooth_oro='.FALSE.' + ilow_pass_oro=1 + numfilt_oro=2 + eps_filter=1.7 + + # soil: tiles + itile_mode=1 + + # model grid type + model_grid_type=1 + name_model_grid='INPUT_ICON_GRID' + +elif [[ $model == "i2" ]]; then + + #output file names + netcdf_output_filename="external_parameter_icon_grid_0002_R19B07_mch.nc" + + # grid definition + grid_dir="/oprusers/osm/opr/data/grid_descriptions/" + grid_nc="icon_grid_0002_R19B07_mch.nc" + + lsso_param='.TRUE.' + lsubtract_mean_slope='.FALSE.' + + # orography raw data + ntiles_column=2 + ntiles_row=2 + topo_files="'Copernicus_DEM_N60-N50_W020-E000.nc' 'Copernicus_DEM_N60-N50_E000-E020.nc' 'Copernicus_DEM_N50-N40_W020-E000.nc' 'Copernicus_DEM_N50-N40_E000-E020.nc'" + + #orography smoothing + lsmooth_oro='.FALSE.' + ilow_pass_oro=1 + numfilt_oro=2 + eps_filter=1.7 + + # soil: tiles + itile_mode=1 + + # model grid type + model_grid_type=1 + name_model_grid='INPUT_ICON_GRID' + + +else + echo $model + echo " Please specify one of the following models: c1, c2, i1, i2, i1_dev, i2_dev" + exit + +fi + +#-------------------------------------------------------------------------------- +# define paths and variables independent from host or model +#-------------------------------------------------------------------------------- + +# Names of executables + +# python executables +binary_alb=extpar_alb_to_buffer.py +binary_ndvi=extpar_ndvi_to_buffer.py +binary_emiss=extpar_emiss_to_buffer.py +binary_era=extpar_era_to_buffer.py +binary_ahf=extpar_ahf_to_buffer.py +binary_isa=extpar_isa_to_buffer.py +binary_tclim=extpar_cru_to_buffer.py + +# fortran executables +binary_lu=extpar_landuse_to_buffer.exe +binary_topo=extpar_topo_to_buffer.exe +binary_aot=extpar_aot_to_buffer.py +binary_soil=extpar_soil_to_buffer.exe +binary_flake=extpar_flake_to_buffer.exe +binary_consistency_check=extpar_consistency_check.exe + +# Names of raw data for INPUT_* namelists +raw_data_alb='alb_new.nc' +raw_data_alnid='alnid_new.nc' +raw_data_aluvd='aluvd_new.nc' +buffer_alb='month_alb_buffer.nc' + +raw_data_aot='aot_GACP.nc' +buffer_aot='extpar_buffer_aot.nc' + +raw_data_tclim_coarse='CRU_T2M_SURF_clim.nc' +raw_data_tclim_fine='CRU_T_SOIL_clim.nc' +buffer_tclim='crutemp_clim_extpar_buffer.nc' + +raw_data_glc2000='GLC2000_byte.nc' +buffer_glc2000='extpar_landuse_buffer.nc' +raw_data_glcc='GLCC_usgs_class_byte.nc' +buffer_glcc='glcc_landuse_buffer.nc' + +raw_data_globcover_0='GLOBCOVER_0_16bit.nc' +raw_data_globcover_1='GLOBCOVER_1_16bit.nc' +raw_data_globcover_2='GLOBCOVER_2_16bit.nc' +raw_data_globcover_3='GLOBCOVER_3_16bit.nc' +raw_data_globcover_4='GLOBCOVER_4_16bit.nc' +raw_data_globcover_5='GLOBCOVER_5_16bit.nc' + +buffer_lu='extpar_landuse_buffer.nc' + +# lanczos filter is recommended when activating scale separation +raw_filt_globe_A10='GLOBE_A_filt_lanczos_window.nc' +raw_filt_globe_B10='GLOBE_B_filt_lanczos_window.nc' +raw_filt_globe_C10='GLOBE_C_filt_lanczos_window.nc' +raw_filt_globe_D10='GLOBE_D_filt_lanczos_window.nc' +raw_filt_globe_E10='GLOBE_E_filt_lanczos_window.nc' +raw_filt_globe_F10='GLOBE_F_filt_lanczos_window.nc' +raw_filt_globe_G10='GLOBE_G_filt_lanczos_window.nc' +raw_filt_globe_H10='GLOBE_H_filt_lanczos_window.nc' +raw_filt_globe_I10='GLOBE_I_filt_lanczos_window.nc' +raw_filt_globe_J10='GLOBE_J_filt_lanczos_window.nc' +raw_filt_globe_K10='GLOBE_K_filt_lanczos_window.nc' +raw_filt_globe_L10='GLOBE_L_filt_lanczos_window.nc' +raw_filt_globe_M10='GLOBE_M_filt_lanczos_window.nc' +raw_filt_globe_N10='GLOBE_N_filt_lanczos_window.nc' +raw_filt_globe_O10='GLOBE_O_filt_lanczos_window.nc' +raw_filt_globe_P10='GLOBE_P_filt_lanczos_window.nc' + + +buffer_topo='topography_buffer.nc' +output_topo="topography_${model}.nc" + +raw_data_ndvi='NDVI_1998_2003.nc' +buffer_ndvi='ndvi_buffer.nc' + +raw_data_soil_FAO='FAO_DSMW_double.nc' +raw_data_soil_HWSD='HWSD0_30_topsoil.nc' +raw_data_deep_soil='HWSD30_100_subsoil.nc' +buffer_soil='soil_buffer.nc' + +raw_lookup_table_HWSD='LU_TAB_HWSD_UF.data' +raw_HWSD_data='HWSD_DATA_ICON.data' +raw_HWSD_data_deep='HWSD_DATA_ICON_S.data' +raw_HWSD_data_extpar='HWSD_DATA_ICON_EXTPAR.asc' + +raw_data_flake='GLDB_lakedepth.nc' +buffer_flake='flake_buffer.nc' + +#-------------------------------------------------------------------------------- +# Prepare working directory and create namelists +#-------------------------------------------------------------------------------- + +if [[ ! -d ${sandboxdir} ]] ; then + mkdir -p ${sandboxdir} +fi + +# link raw data files to local workdir +ln -s -f ${data_dir}/*.nc ${sandboxdir} + +cp $scriptdir/* ${sandboxdir}/. +cp $exedir/* ${sandboxdir}/. +cd ${sandboxdir} +if [[ -e ${logfile} ]] ; then + rm -f ${logfile} +fi + +cd ${sandboxdir} + +echo "\n>>>> Data will be processed and produced in `pwd` <<<<\n" + +echo PYTHONPATH: ${PYTHONPATH} >> ${logfile} + +# create input namelists + +cat > namelist.py << EOF_namelist_python +input_alb = { + 'ialb_type': 1, + 'raw_data_alb_path': '', + 'raw_data_alb_filename': '${raw_data_alb}', + 'raw_data_alnid_filename': '${raw_data_alnid}', + 'raw_data_aluvd_filename': '${raw_data_aluvd}', + 'alb_buffer_file': '${buffer_alb}', + } + +input_aot = { + 'iaot_type': 1, + 'raw_data_aot_path': '', + 'raw_data_aot_filename': '${raw_data_aot}', + 'aot_buffer_file': '${buffer_aot}', + } + +input_tclim = { + 'raw_data_t_clim_path': '', + 'raw_data_tclim_coarse': '', + 'raw_data_tclim_fine': '${raw_data_tclim_fine}', + 't_clim_buffer_file': '${buffer_tclim}', + 'it_cl_type': 1 + } + +input_ndvi = { + 'raw_data_ndvi_path': '', + 'raw_data_ndvi_filename': '${raw_data_ndvi}', + 'ndvi_buffer_file': '${buffer_ndvi}', + } + +input_emiss = { + 'iemiss_type': 2, + 'raw_data_emiss_path': '', + 'raw_data_emiss_filename': 'CAMEL_bbe_lw_2010-2015.nc', + 'emiss_buffer_file': 'emiss_buffer.nc' + } + +input_era = { + 'iera_type': 1, + 'raw_data_era_path': '', + 'raw_data_era_ORO': 'ERA5_ORO_1990.nc', + 'raw_data_era_T2M': 'ERA5_T2M_1990_2019.nc', + 'raw_data_era_SST': 'ERA5_SST_1990_2019.nc', + 'raw_data_era_SD': 'ERA5_SD_1990_2019.nc', + 'era_buffer_file': 'era_buffer.nc', + } + +input_ahf = { + 'iahf_type': 1, + 'raw_data_ahf_path': '', + 'raw_data_ahf_filename': 'AHF_2006_CDO.nc', + 'ahf_buffer_file': 'ahf_buffer.nc', + } + +input_isa = { + 'raw_data_isa_path': '', + 'raw_data_isa_filename': 'NOAA_ISA_CDO.nc', + 'isa_buffer_file': 'isa_buffer.nc', + } +EOF_namelist_python + +# set target grid definition +cat > INPUT_grid_org << EOF_go +&GRID_DEF + igrid_type=${model_grid_type}, + domain_def_namelist='${name_model_grid}' +/ +EOF_go + +# COSMO grid +cat > INPUT_COSMO_GRID << EOF_grid_cosmo +&lmgrid + pollon=${rot_pol_lon}, + pollat=${rot_pol_lat}, + startlon_tot=${startlon_tot}, + startlat_tot=${startlat_tot}, + dlon=${dlon}, + dlat=${dlat}, + ie_tot=${ie_tot}, + je_tot=${je_tot}, +/ +EOF_grid_cosmo + +# ICON grid +cat > INPUT_ICON_GRID << EOF_grid_icon +&icon_grid_info + icon_grid_dir="${grid_dir}", + icon_grid_nc_file="${grid_nc}", +/ +EOF_grid_icon + +cat > INPUT_LU << EOF_lu +&lu_raw_data + raw_data_lu_path='', + raw_data_lu_filename='${raw_data_globcover_0}' '${raw_data_globcover_1}' '${raw_data_globcover_2}' '${raw_data_globcover_3}' '${raw_data_globcover_4}' '${raw_data_globcover_5}', + i_landuse_data=1, + ilookup_table_lu=1 +/ +&lu_io_extpar + lu_buffer_file='${buffer_lu}', +/ +&glcc_raw_data + raw_data_glcc_path='', + raw_data_glcc_filename='${raw_data_glcc}' +/ +&glcc_io_extpar + glcc_buffer_file='${buffer_glcc}', +/ +EOF_lu + +cat > INPUT_ORO << EOF_oro +&oro_runcontrol + lcompute_sgsl=.FALSE. , + / +&orography_io_extpar + orography_buffer_file='${buffer_topo}', + orography_output_file='${output_topo}' +/ +&orography_raw_data +itopo_type=4, +lsso_param=${lsso_param}, +raw_data_orography_path='./', +ntiles_column=${ntiles_column}, +ntiles_row=${ntiles_row}, +topo_files=${topo_files} +/ +EOF_oro + +cat > INPUT_OROSMOOTH << EOF_orosm +&orography_smoothing + lfilter_oro=${lsmooth_oro}, + ilow_pass_oro=${ilow_pass_oro}, + numfilt_oro=${numfilt_oro}, + ilow_pass_xso=5, + lxso_first=.FALSE., + numfilt_xso=1, + rxso_mask=750.0, + eps_filter=${eps_filter}, + rfill_valley=0.0, + ifill_valley=2 +/ +EOF_orosm + +cat > INPUT_RADTOPO << EOF_radtopo +&radtopo + lradtopo=.TRUE. + itype_scaling=2 + nhori=24 + radius=40000 + min_circ_cov=2 + max_missing=0.95 +/ +EOF_radtopo + +cat > INPUT_SCALE_SEP << EOF_scale_sep +&scale_separated_raw_data + lscale_separation = .FALSE., + raw_data_scale_sep_path = '', + scale_sep_files = '${raw_filt_globe_A10}' '${raw_filt_globe_B10}' '${raw_filt_globe_C10}' '${raw_filt_globe_D10}' '${raw_filt_globe_E10}' '${raw_filt_globe_F10}' '${raw_filt_globe_G10}' '${raw_filt_globe_H10}' '${raw_filt_globe_I10}' '${raw_filt_globe_J10}' '${raw_filt_globe_K10}' '${raw_filt_globe_L10}' '${raw_filt_globe_M10}' '${raw_filt_globe_N10}' '${raw_filt_globe_O10}' '${raw_filt_globe_P10}' +/ +EOF_scale_sep + +cat > INPUT_SOIL << EOF_soil +&soil_raw_data +isoil_data = 3, +ldeep_soil = .FALSE., +raw_data_soil_path='', +raw_data_soil_filename='${raw_data_soil_HWSD}', +/ +&soil_io_extpar + soil_buffer_file='${buffer_soil}', +/ +&HWSD_index_files +path_HWSD_index_files='', +lookup_table_HWSD='${raw_lookup_table_HWSD}', +HWSD_data='${raw_HWSD_data}', +HWSD_data_deep='${raw_HWSD_data_deep}', +HWSD_data_extpar='${raw_HWSD_data_extpar}' +/ +EOF_soil + +cat > INPUT_FLAKE << EOF_flake +&flake_raw_data + raw_data_flake_path='', + raw_data_flake_filename='${raw_data_flake}' +/ +&flake_io_extpar + flake_buffer_file='${buffer_flake}' +/ +EOF_flake + +# consistency check +cat > INPUT_CHECK << EOF_check +&extpar_consistency_check_io + netcdf_output_filename="${netcdf_output_filename}", + tile_mode=${itile_mode}, + lwrite_netcdf=.TRUE. + i_lsm_data=1, + land_sea_mask_file="", + number_special_points=0, +/ +EOF_check + +#-------------------------------------------------------------------------------- +# launch extpar executables +#-------------------------------------------------------------------------------- + +# 1. run topography first (is needed for CRU data processing) +#-------------------------------- +run_sequential ${binary_topo} # topography + +# 2. run rest which is needed for ICON +#-------------------------------- +run_sequential ${binary_alb} # albedo +run_sequential ${binary_aot} # aerosol optical thickness +run_sequential ${binary_emiss} # emissivity +run_sequential ${binary_era} # era climatology +run_sequential ${binary_flake} # fraction lake +run_sequential ${binary_lu} # land use +run_sequential ${binary_ndvi} # normalized difference vegetation index +run_sequential ${binary_tclim} # cru: temperature climatology +run_sequential ${binary_soil} # soil + + +# 3. binaries which are normally not needed for icon-nwp +#-------------------------------- +# run_sequential ${binary_ahf} # anthropogenic heat flux +# run_sequential ${binary_isa} # impervious surface area + + +# 4. finally: run consistency check +#-------------------------------- +run_sequential ${binary_consistency_check} + +#-------------------------------------------------------------------------------- +# clean-up +#-------------------------------------------------------------------------------- +rm -f exit_status_* +rm -f time_* + +echo ">>>> External parameters for COSMO/ICON model generated <<<<" diff --git a/src/extpar_topo_to_buffer.f90 b/src/extpar_topo_to_buffer.f90 index 257fb832..ade7ada2 100644 --- a/src/extpar_topo_to_buffer.f90 +++ b/src/extpar_topo_to_buffer.f90 @@ -6,7 +6,7 @@ ! V1_0 2010/12/21 Hermann Asensio ! Initial release ! V1_1 2011/01/20 Hermann Asensio -! small bug fixes accroding to Fortran compiler warnings +! small bug fixes according to Fortran compiler warnings ! V1_2 2011/03/25 Hermann Asensio ! update to support ICON refinement grids ! V1_4 2011/04/21 Anne Roches @@ -37,7 +37,7 @@ !> !! @par extpar_topo_to_buffer !! -!! This program reads in the GLOBE/ASTER/MERIT orography data set and aggregates the orographic height to the target grid +!! This program reads in the orography data set and aggregates the orographic height to the target grid !! and computes the subgrid-scale orography parameters (SSO) required by the SSO-parameterization. !! !> Purpose: read in GLOBE/ASTER orography data and aggregate to COSMO grid @@ -86,6 +86,7 @@ PROGRAM extpar_topo_to_buffer USE mo_topo_data, ONLY: topo_aster, & & topo_merit, & + & topo_copernicus, & & itopo_type, & & topo_tiles_grid, & & topo_grid, & @@ -106,6 +107,10 @@ PROGRAM extpar_topo_to_buffer & merit_lat_max, & & merit_lon_min, & & merit_lon_max, & + & copernicus_lat_min, & + & copernicus_lat_max, & + & copernicus_lon_min, & + & copernicus_lon_max, & & num_tiles, & & allocate_topo_data,& & fill_topo_data, & @@ -361,8 +366,20 @@ PROGRAM extpar_topo_to_buffer CALL logging%warning(message_text) CALL logging%error('The chosen latitude edges are not within the MERIT domain.',__FILE__,__LINE__) END IF - - + CASE(topo_copernicus) + WRITE(message_text,*)'Edges of domain: ', copernicus_lon_min,' ', copernicus_lon_max,' ', copernicus_lat_min,' ' & + & ,copernicus_lat_max + CALL logging%info(message_text) + IF (lon_geo (tg%ie,tg%je,tg%ke) > copernicus_lon_max .OR. lon_geo(1,1,1) < copernicus_lon_min) THEN + WRITE(message_text,*) 'COPERNICUS min lon is: ', copernicus_lon_min, ' and COPERNICUS max lon is: ', copernicus_lon_max + CALL logging%warning(message_text) + CALL logging%error('The chosen longitude edges are not within the COPERNICUS domain.',__FILE__,__LINE__) + END IF + IF (lat_geo(tg%ie,tg%je,tg%ke) > copernicus_lat_max .OR. lat_geo(1,1,1) < copernicus_lat_min) THEN + WRITE(message_text,*) 'COPERNICUS min lat is: ', copernicus_lat_min, ' and COPERNICUS max lat is: ', copernicus_lat_max + CALL logging%warning(message_text) + CALL logging%error('The chosen latitude edges are not within the COPERNICUS domain.',__FILE__,__LINE__) + END IF END SELECT CALL det_topo_tiles_grid(topo_tiles_grid) diff --git a/src/mo_agg_sgsl.f90 b/src/mo_agg_sgsl.f90 index 7456bdbd..e57fa0e3 100644 --- a/src/mo_agg_sgsl.f90 +++ b/src/mo_agg_sgsl.f90 @@ -57,6 +57,7 @@ MODULE mo_agg_sgsl & topo_gl, & & topo_aster, & & topo_merit, & + & topo_copernicus, & & get_fill_value_sgsl, & & nr_tot !< total number of rows in GLOBE/ASTER data @@ -388,18 +389,23 @@ SUBROUTINE agg_sgsl_data_to_target_grid(sgsl_tiles_grid, & SELECT CASE(itopo_type) CASE(topo_aster) IF (sl(i) /= default_sgsl) THEN - ndata(ie,je,ke) = ndata(ie,je,ke) + 1 - sgsl(ie,je,ke) = sgsl(ie,je,ke) + sl(i) + ndata(ie,je,ke) = ndata(ie,je,ke) + 1 + sgsl(ie,je,ke) = sgsl(ie,je,ke) + sl(i) ENDIF CASE(topo_gl) IF (sl(i) /= undef_sgsl) THEN - ndata(ie,je,ke) = ndata(ie,je,ke) + 1 - sgsl(ie,je,ke) = sgsl(ie,je,ke) + sl(i) + ndata(ie,je,ke) = ndata(ie,je,ke) + 1 + sgsl(ie,je,ke) = sgsl(ie,je,ke) + sl(i) ENDIF CASE(topo_merit) IF (sl(i) /= undef_sgsl) THEN - ndata(ie,je,ke) = ndata(ie,je,ke) + 1 - sgsl(ie,je,ke) = sgsl(ie,je,ke) + sl(i) + ndata(ie,je,ke) = ndata(ie,je,ke) + 1 + sgsl(ie,je,ke) = sgsl(ie,je,ke) + sl(i) + ENDIF + CASE(topo_copernicus) + IF (sl(i) /= undef_sgsl) THEN + ndata(ie,je,ke) = ndata(ie,je,ke) + 1 + sgsl(ie,je,ke) = sgsl(ie,je,ke) + sl(i) ENDIF END SELECT !$OMP END CRITICAL diff --git a/src/mo_agg_topo_icon.f90 b/src/mo_agg_topo_icon.f90 index 6bbe306e..4cd4b690 100644 --- a/src/mo_agg_topo_icon.f90 +++ b/src/mo_agg_topo_icon.f90 @@ -68,15 +68,16 @@ MODULE mo_agg_topo_icon USE mo_icon_domain, ONLY: icon_domain, grid_cells - USE mo_topo_data, ONLY: ntiles, & !< there are 16/240 GLOBE/ASTER/MERIT tiles + USE mo_topo_data, ONLY: ntiles, & & max_tiles, & - & nc_tot, & !< number of total GLOBE/ASTER columns un a latitude circle - & nr_tot, & !< total number of rows in GLOBE/ASTER/MERIT data + & nc_tot, & !< number of total topography dataset columns un a latitude circle + & nr_tot, & !< total number of rows in topography dataset & get_fill_value, & & itopo_type, & & topo_gl, & & topo_aster, & & topo_merit, & + & topo_copernicus, & & get_varname USE mo_topo_sso, ONLY: auxiliary_sso_parameter_icon,& @@ -272,6 +273,9 @@ SUBROUTINE agg_topo_data_to_target_grid_icon(topo_tiles_grid, & CASE(topo_merit) hh = default_topo hh_red = default_topo + CASE(topo_copernicus) + hh = default_topo + hh_red = default_topo END SELECT ! initialize some variables diff --git a/src/mo_extpar_output_nc.f90 b/src/mo_extpar_output_nc.f90 index d7e6c943..fe1b6d0b 100644 --- a/src/mo_extpar_output_nc.f90 +++ b/src/mo_extpar_output_nc.f90 @@ -76,7 +76,7 @@ MODULE mo_extpar_output_nc USE mo_soil_data, ONLY: HWSD_data - USE mo_topo_data, ONLY: itopo_type, topo_aster, topo_gl, topo_merit + USE mo_topo_data, ONLY: itopo_type, topo_aster, topo_gl, topo_merit, topo_copernicus USE mo_cosmo_grid, ONLY: lon_rot, lat_rot @@ -1718,6 +1718,9 @@ SUBROUTINE set_cdi_global_att_icon(global_attributes,itopo_type,name_lookup_tabl CASE(3) ! topo_merit global_attributes(3)%attributetext=TRIM(lu_dataset)//', FAO DSMW, MERIT, Lake Database' IF(isoil_data == HWSD_data) global_attributes(3)%attributetext=TRIM(lu_dataset)//', HWSD, MERIT, Lake Database' + CASE(4) ! topo_copernicus + global_attributes(3)%attributetext=TRIM(lu_dataset)//', FAO DSMW, COPERNICUS, Lake Database' + IF(isoil_data == HWSD_data) global_attributes(3)%attributetext=TRIM(lu_dataset)//', HWSD, COPERNICUS, Lake Database' END SELECT @@ -1788,6 +1791,12 @@ SUBROUTINE set_global_att_extpar(global_attributes,name_lookup_table_lu,lu_datas ELSE global_attributes(3)%attributetext=TRIM(lu_dataset)//', FAO DSMW, MERIT, Lake Database' ENDIF + CASE(topo_copernicus) + IF (isoil_data >= HWSD_data) THEN + global_attributes(3)%attributetext=TRIM(lu_dataset)//', HWSD, COPERNICUS, Lake Database' + ELSE + global_attributes(3)%attributetext=TRIM(lu_dataset)//', FAO DSMW, COPERNICUS, Lake Database' + ENDIF END SELECT global_attributes(4)%attname = 'note' global_attributes(4)%attributetext='Landuse data look-up table: '//TRIM(name_lookup_table_lu) diff --git a/src/mo_preproc_for_sgsl.f90 b/src/mo_preproc_for_sgsl.f90 index 7171b9ed..9d9a714e 100644 --- a/src/mo_preproc_for_sgsl.f90 +++ b/src/mo_preproc_for_sgsl.f90 @@ -26,7 +26,7 @@ MODULE mo_preproc_for_sgsl CONTAINS ! wrapper function for the preprocessing of raw orography data - ! and calls the right subroutine for itopo_type (GLOBE, ASTER or MERIT) + ! and calls the right subroutine for itopo_type (GLOBE, ASTER, MERIT or COPERNICUS) SUBROUTINE preproc_orography (raw_data_orography_path, & & topo_files, & & sgsl_files, & @@ -99,19 +99,21 @@ SUBROUTINE topo_subgrid_slope( infile, outfile, itopo_type) REAL(KIND=wp) :: dx, dy, oolenx, ooleny, grad(9), & & dx0, dx2, oolen0, oolen2, dlat, dlon - REAL(KIND=wp), PARAMETER :: r_earth = 6371.229E3, & ! radius of the earth - & pi = 4.0 * ATAN (1.0), & - & degrad = pi / 180.0, & - & dlat_globe = 30./3600. , &! resolution - & dlon_globe = 30./3600. , &! resolution - & dlat_aster = 1./3600. , &! resolution - & dlon_aster = 1./3600. , &! resolution - & dlat_merit = 3./3600. , &! resolution - & dlon_merit = 3./3600. , &! resolution - & eps = 1.E-9, & - & add_offset = 0., & - & scale_factor = 0.001, & - & r_scfct = 1. / scale_factor + REAL(KIND=wp), PARAMETER :: r_earth = 6371.229E3, & ! radius of the earth + & pi = 4.0 * ATAN (1.0), & + & degrad = pi / 180.0, & + & dlat_globe = 30./3600. , &! resolution + & dlon_globe = 30./3600. , &! resolution + & dlat_aster = 1./3600. , &! resolution + & dlon_aster = 1./3600. , &! resolution + & dlat_merit = 3./3600. , &! resolution + & dlon_merit = 3./3600. , &! resolution + & dlat_copernicus = 1./3600. , &! resolution + & dlon_copernicus = 1./3600. , &! resolution + & eps = 1.E-9, & + & add_offset = 0., & + & scale_factor = 0.001, & + & r_scfct = 1. / scale_factor WRITE(message_text,*) TRIM(infile), ' --> ', TRIM(outfile) @@ -160,8 +162,10 @@ SUBROUTINE topo_subgrid_slope( infile, outfile, itopo_type) status = nf90_inq_varid(ncid,"altitude", varid) ELSE IF (itopo_type == 2) THEN status = nf90_inq_varid(ncid,"Z", varid) - ELSE + ELSE IF (itopo_type == 3) THEN status = nf90_inq_varid(ncid,"Elevation", varid) + ELSE + status = nf90_inq_varid(ncid,"elevation", varid) ENDIF CALL check_err(status, __FILE__, __LINE__) @@ -186,9 +190,12 @@ SUBROUTINE topo_subgrid_slope( infile, outfile, itopo_type) ELSE IF (itopo_type == 2) THEN !ASTER dlat = dlat_aster dlon = dlon_aster - ELSE !MERIT + ELSE IF (itopo_type == 3) THEN !MERIT dlat = dlat_merit dlon = dlon_merit + ELSE !COPERNICUS + dlat = dlat_copernicus + dlon = dlon_copernicus ENDIF !itopo_type ! compute values for inner domain @@ -467,8 +474,10 @@ SUBROUTINE topo_subgrid_slope( infile, outfile, itopo_type) status = nf90_put_att(ncido, outid,'comment',trim(comment)) ELSE IF (itopo_type == 2) THEN status = nf90_put_att(ncido, outid,'comment','ASTER tile: '//infile) + ELSE IF (itopo_type == 3) THEN + status = nf90_put_att(ncido, outid,'comment','MERIT tile: '//infile) ELSE - status = nf90_put_att(ncido, outid,'comment','MERIT tile: '//infile) + status = nf90_put_att(ncido, outid,'comment','COPERNICUS tile: '//infile) ENDIF CALL check_err(status, __FILE__, __LINE__) diff --git a/src/mo_topo_data.f90 b/src/mo_topo_data.f90 index 562f4060..4c5844a2 100644 --- a/src/mo_topo_data.f90 +++ b/src/mo_topo_data.f90 @@ -87,6 +87,7 @@ MODULE mo_topo_data & topo_gl, & & topo_aster, & & topo_merit, & + & topo_copernicus, & & aster_lat_min, & & aster_lon_min, & & aster_lat_max, & @@ -95,6 +96,10 @@ MODULE mo_topo_data & merit_lon_min, & & merit_lat_max, & & merit_lon_max, & + & copernicus_lat_min, & + & copernicus_lon_min, & + & copernicus_lat_max, & + & copernicus_lon_max, & & ntiles_row, & & ntiles_column, & & lradtopo, & @@ -130,6 +135,7 @@ MODULE mo_topo_data INTEGER(KIND=i4), PARAMETER :: topo_gl = 1, & & topo_aster = 2, & & topo_merit = 3, & + & topo_copernicus = 4, & & max_tiles = 1000 REAL(KIND=wp), ALLOCATABLE :: tiles_lon_min(:), & @@ -151,6 +157,11 @@ MODULE mo_topo_data & merit_lon_min, & & merit_lon_max + REAL(KIND=wp):: copernicus_lat_min, & + & copernicus_lat_max, & + & copernicus_lon_min, & + & copernicus_lon_max + LOGICAL :: lradtopo CHARACTER(LEN=80) :: varname @@ -162,7 +173,7 @@ SUBROUTINE num_tiles(columns, rows, ntiles) ! it gives the value of the number o SAVE INTEGER(KIND=i4), INTENT(IN) :: columns, & & rows - INTEGER, INTENT(OUT) :: ntiles ! if the user chooses GLOBE, ASTER, MERIT + INTEGER, INTENT(OUT) :: ntiles ntiles_column = columns ntiles_row = rows @@ -175,7 +186,7 @@ SUBROUTINE allocate_topo_data(ntiles) IMPLICIT NONE - INTEGER, INTENT (IN) :: ntiles ! number of tiles: 36 for ASTER and 16 for GLOBE + INTEGER, INTENT (IN) :: ntiles INTEGER :: errorcode CALL logging%info('Enter routine: allocate_topo_data') @@ -214,6 +225,10 @@ SUBROUTINE allocate_topo_data(ntiles) merit_lon_min = 0.0 merit_lon_max = 0.0 + copernicus_lat_min = 0.0 + copernicus_lat_max = 0.0 + copernicus_lon_min = 0.0 + copernicus_lon_max = 0.0 END SUBROUTINE allocate_topo_data @@ -241,23 +256,23 @@ SUBROUTINE fill_topo_data(raw_data_orography_path,topo_files, & INTEGER(KIND=i4) :: i, errorcode, ncid, & & dimID_lat, dimID_lon, varID_lat, varID_lon - REAL(KIND=wp) :: half_gridp! distance of half a grid point as the grid point is centered on a GLOBE / ASTER pixel + REAL(KIND=wp) :: half_gridp ! distance of half a grid cell CALL logging%info('Enter routine: fill_topo_data') - SELECT CASE (itopo_type) ! Also topo could additionally be used for SELECT CASE (must first be read in) - CASE(topo_aster) ! ASTER topography, as it has 36 tiles at the moment. + SELECT CASE (itopo_type) + CASE(topo_aster) ! ASTER topography: 240 tiles CALL logging%info('ASTER is used as topography') - half_gridp = 1./(3600.*2.) ! the resolution of the ASTER data is 1./3600. degrees as it is half a grid point - ! it is additionally divided by 2 - CASE (topo_gl) ! GLOBE topography is composed of 16 tiles + half_gridp = 1./(3600.*2.) ! half a ASTER grid cell [degrees] + CASE (topo_gl) ! GLOBE topography: 16 tiles CALL logging%info('GLOBE is used as topography') - half_gridp = 1./(120.*2.) ! GLOBE resolution is 1./120. degrees (30 arc-seconds) - - CASE(topo_merit) ! ASTER topography, as it has 36 tiles at the moment. + half_gridp = 1./(120.*2.) ! half a GLOBE grid cell [degrees] + CASE(topo_merit) ! MERIT topography: 60 tiles CALL logging%info('MERIT is used as topography') - half_gridp = 1./(1200.*2.) ! the resolution of the MERIT data is 1./1200. degrees as it is half a grid point - ! it is additionally divided by 2 + half_gridp = 1./(1200.*2.) ! half a MERIT grid cell [degrees] + CASE(topo_copernicus) ! Copernicus topography: 324 tiles + CALL logging%info('COPERNICUS is used as topography') + half_gridp = 1./(3600.*2.) ! half a Copernicus grid cell [degrees] END SELECT DO i = 1,ntiles @@ -278,10 +293,19 @@ SUBROUTINE fill_topo_data(raw_data_orography_path,topo_files, & ! reads in the last latitude value of tile i CALL check_netcdf(nf90_close(ncid)) ! the netcdf file is closed again - tiles_lon_min(i) = REAL(NINT(tiles_lon_min(i) - half_gridp)) !< half of a grid point must be - tiles_lon_max(i) = REAL(NINT(tiles_lon_max(i) + half_gridp)) !< added, as the ASTER/GLOBE/MERIT data - tiles_lat_min(i) = REAL(NINT(tiles_lat_min(i) + half_gridp)) !< is located at the pixel center - tiles_lat_max(i) = REAL(NINT(tiles_lat_max(i) - half_gridp)) + SELECT CASE (itopo_type) + ! compute domain extent of tiles (coordinates refer to domain edges) from pixel center coordinates + CASE(topo_aster, topo_gl) ! edges of raw data tiles align with integer coordinates + tiles_lon_min(i) = REAL(NINT(tiles_lon_min(i) - half_gridp)) + tiles_lon_max(i) = REAL(NINT(tiles_lon_max(i) + half_gridp)) + tiles_lat_min(i) = REAL(NINT(tiles_lat_min(i) - half_gridp)) + tiles_lat_max(i) = REAL(NINT(tiles_lat_max(i) + half_gridp)) + CASE(topo_merit, topo_copernicus) ! edges of raw data tiles do not align with integer coordinates + tiles_lon_min(i) = REAL(tiles_lon_min(i) - half_gridp) + tiles_lon_max(i) = REAL(tiles_lon_max(i) + half_gridp) + tiles_lat_min(i) = REAL(tiles_lat_min(i) - half_gridp) + tiles_lat_max(i) = REAL(tiles_lat_max(i) + half_gridp) + END SELECT END DO SELECT CASE(itopo_type) @@ -296,6 +320,12 @@ SUBROUTINE fill_topo_data(raw_data_orography_path,topo_files, & merit_lat_max = MAXVAL(tiles_lat_max) merit_lon_min = MINVAL(tiles_lon_min) merit_lon_max = MAXVAL(tiles_lon_max) + + CASE(topo_copernicus) + copernicus_lat_min = MINVAL(tiles_lat_min) + copernicus_lat_max = MAXVAL(tiles_lat_max) + copernicus_lon_min = MINVAL(tiles_lon_min) + copernicus_lon_max = MAXVAL(tiles_lon_max) END SELECT nc_tot = 0 @@ -350,28 +380,26 @@ SUBROUTINE get_fill_value(topo_file_1, undef_topo) CHARACTER (len=*), INTENT(in) :: topo_file_1 INTEGER(KIND=i4), INTENT(out) :: undef_topo - INTEGER(KIND=i4) :: ncid, varid, status + INTEGER(KIND=i4) :: ncid, varid, status, i + CHARACTER(len=10) :: var_names(4) + LOGICAL :: vn_found = .FALSE. CALL check_netcdf(nf90_open(path = topo_file_1, mode = nf90_nowrite, ncid = ncid)) - status = nf90_inq_varid(ncid, "altitude", varid) - IF (status == NF90_ENOTVAR) THEN - status = nf90_inq_varid(ncid, "Z", varid) - IF (status == NF90_ENOTVAR) THEN - status = nf90_inq_varid(ncid, "Elevation", varid) - IF (status == NF90_ENOTVAR) THEN - WRITE(message_text,*)'Could not find "altitude (GLOBE)" or "Z (ASTER)" & - & or "Elevation (MERIT/REMA)" in topography file ' & - & //TRIM(topo_file_1) - CALL logging%error(message_text,__FILE__,__LINE__) - ELSE - CALL check_netcdf(status, __FILE__, __LINE__) - ENDIF - ELSE - CALL check_netcdf(status, __FILE__, __LINE__) - END IF - ELSE - CALL check_netcdf(status, __FILE__, __LINE__) - ENDIF + var_names = [character(len=10) :: "altitude", "Z", "Elevation", "elevation"] + DO i = 1, SIZE(var_names) + status = nf90_inq_varid(ncid, TRIM(var_names(i)), varid) + IF (status /= NF90_ENOTVAR) THEN + vn_found = .TRUE. + CALL check_netcdf(status, __FILE__, __LINE__) + EXIT + END IF + END DO + IF (.NOT. vn_found) THEN + WRITE(message_text,*)'Could not find "altitude" (GLOBE), "Z" (ASTER), & + & "Elevation" (MERIT/REMA) or "elevation" (COPERNICUS) in topography file ' & + & //TRIM(topo_file_1) + CALL logging%error(message_text,__FILE__,__LINE__) + END IF CALL check_netcdf(nf90_get_att(ncid, varid, "_FillValue", undef_topo), __FILE__, __LINE__) CALL check_netcdf(nf90_close(ncid)) @@ -397,6 +425,10 @@ SUBROUTINE get_varname(topo_file_1,varname) CALL check_netcdf(nf90_open(path = trim(topo_file_1), mode = nf90_nowrite, ncid = ncid)) CALL check_netcdf(nf90_inquire_variable(ncid,3,varname,type,ndims,dimids)) CALL check_netcdf(nf90_close(ncid)) + CASE(topo_copernicus) + CALL check_netcdf(nf90_open(path = trim(topo_file_1), mode = nf90_nowrite, ncid = ncid)) + CALL check_netcdf(nf90_inquire_variable(ncid,3,varname,type,ndims,dimids)) + CALL check_netcdf(nf90_close(ncid)) END SELECT END SUBROUTINE get_varname @@ -491,6 +523,13 @@ SUBROUTINE get_fill_value_sgsl(sgsl_file_1,undef_sgsl) undef_sgsl = fillval * scale_factor CALL check_netcdf(nf90_close(ncid)) + CASE(topo_copernicus) + CALL check_netcdf(nf90_open(path = sgsl_file_1, mode = nf90_nowrite, ncid = ncid)) + CALL check_netcdf(nf90_get_att(ncid, 3, "_FillValue", fillval)) + CALL check_netcdf(nf90_get_att(ncid, 3, "scale_factor", scale_factor)) + undef_sgsl = fillval * scale_factor + CALL check_netcdf(nf90_close(ncid)) + END SELECT END SUBROUTINE get_fill_value_sgsl @@ -526,6 +565,11 @@ SUBROUTINE get_varname_sgsl(sgsl_file_1,varname) CALL check_netcdf(nf90_inquire_variable(ncid,3,varname,type,ndims,dimids)) CALL check_netcdf(nf90_close(ncid)) + CASE(topo_copernicus) + CALL check_netcdf(nf90_open(path = sgsl_file_1, mode = nf90_nowrite, ncid = ncid)) + CALL check_netcdf(nf90_inquire_variable(ncid,3,varname,type,ndims,dimids)) + CALL check_netcdf(nf90_close(ncid)) + END SELECT END SUBROUTINE get_varname_sgsl diff --git a/src/mo_topo_output_nc.f90 b/src/mo_topo_output_nc.f90 index 01733f30..24fbb584 100644 --- a/src/mo_topo_output_nc.f90 +++ b/src/mo_topo_output_nc.f90 @@ -35,7 +35,7 @@ MODULE mo_topo_output_nc USE mo_cosmo_grid, ONLY: cosmo_grid, nborder, lon_rot, lat_rot USE mo_icon_grid_data, ONLY: icon_grid - USE mo_topo_data, ONLY: itopo_type, topo_gl, topo_aster, topo_merit + USE mo_topo_data, ONLY: itopo_type, topo_gl, topo_aster, topo_merit, topo_copernicus USE mo_io_utilities, ONLY: netcdf_attributes, dim_meta_info, & & netcdf_get_var, netcdf_put_var, & @@ -704,6 +704,8 @@ SUBROUTINE set_global_att_topo(global_attributes) global_attributes(1)%attributetext='GLOBE data ' CASE(topo_merit) global_attributes(1)%attributetext='MERIT data ' + CASE(topo_copernicus) + global_attributes(1)%attributetext='COPERNICUS data ' END SELECT global_attributes(2)%attname = 'institution' global_attributes(2)%attributetext='Deutscher Wetterdienst' @@ -717,6 +719,8 @@ SUBROUTINE set_global_att_topo(global_attributes) global_attributes(3)%attributetext='GLOBE, Global Land One-km Base Elevation' CASE(topo_merit) global_attributes(3)%attributetext='MERIT DEM: Multi-Error-Removed Improved-Terrain DEM ' + CASE(topo_copernicus) + global_attributes(3)%attributetext='Copernicus DEM GLO-30 ' END SELECT @@ -738,6 +742,9 @@ SUBROUTINE set_global_att_topo(global_attributes) global_attributes(5)%attributetext='http://www.ngdc.noaa.gov/mgg/topo/globe.html' CASE(topo_merit) global_attributes(5)%attributetext='http://hydro.iis.u-tokyo.ac.jp/~yamadai/MERIT_DEM/ ' + CASE(topo_copernicus) + global_attributes(5)%attributetext='https://dataspace.copernicus.eu/explore-data/data-collections/'// & + & 'copernicus-contributing-missions/collections-description/COP-DEM ' END SELECT global_attributes(6)%attname = 'comment' diff --git a/src/mo_topo_routines.f90 b/src/mo_topo_routines.f90 index b7bba29a..e7af3d79 100644 --- a/src/mo_topo_routines.f90 +++ b/src/mo_topo_routines.f90 @@ -51,6 +51,7 @@ MODULE mo_topo_routines & topo_aster, & & topo_gl, & & topo_merit, & + & topo_copernicus, & & aster_lat_min, & & aster_lat_max, & & aster_lon_min, & @@ -59,6 +60,10 @@ MODULE mo_topo_routines & merit_lat_max, & & merit_lon_min, & & merit_lon_max, & + & copernicus_lat_min, & + & copernicus_lat_max, & + & copernicus_lon_min, & + & copernicus_lon_max, & & get_varname, & & h_tile_row @@ -344,6 +349,23 @@ SUBROUTINE det_topo_grid(topo_grid) topo_grid%start_lat_reg = merit_lat_max + 0.5_wp * dlat ! latitude from north to south, note the negative increment! topo_grid%end_lat_reg = merit_lat_min - 0.5_wp * dlat ! latitude from north to south, note the negative increment! + CASE(topo_copernicus) + + dlon = (copernicus_lon_max - copernicus_lon_min) / REAL(nc_tot,wp) + + dlat = -1. * (copernicus_lat_max - copernicus_lat_min) / REAL(nr_tot,wp) + + WRITE(message_text,*)'Latitude increment for COPERNICUS data, dlat = ', dlat + CALL logging%info(message_text) + + ! latitude from north to south, negative increment + + topo_grid%start_lon_reg = copernicus_lon_min + 0.5_wp * dlon + topo_grid%end_lon_reg = copernicus_lon_max - 0.5_wp * dlon + + topo_grid%start_lat_reg = copernicus_lat_max + 0.5_wp * dlat ! latitude from north to south, note the negative increment! + topo_grid%end_lat_reg = copernicus_lat_min - 0.5_wp * dlat ! latitude from north to south, note the negative increment! + CASE(topo_gl) dlon = 360._wp / REAL(nc_tot,wp) @@ -510,6 +532,10 @@ SUBROUTINE get_topo_tile_block_indices(ta_grid, & m = 1 n = 1 o = ntiles + CASE(topo_copernicus) + m = 1 + n = 1 + o = ntiles CASE(topo_gl) m = 1 n = ntiles/4 diff --git a/src/mo_var_meta_data.f90 b/src/mo_var_meta_data.f90 index 17de4fb2..90787bbd 100644 --- a/src/mo_var_meta_data.f90 +++ b/src/mo_var_meta_data.f90 @@ -2941,6 +2941,7 @@ SUBROUTINE def_topo_meta(diminfo,itopo_type,igrid_type,coordinates,grid_mapping, INTEGER (KIND=i4), PARAMETER :: topo_aster = 2, & & topo_gl = 1, & & topo_merit = 3, & + & topo_copernicus = 4, & & igrid_icon = 1, & & igrid_cosmo = 2 @@ -2956,6 +2957,8 @@ SUBROUTINE def_topo_meta(diminfo,itopo_type,igrid_type,coordinates,grid_mapping, dataset = 'GLOBE' CASE(topo_merit) dataset = 'MERIT' + CASE(topo_copernicus) + dataset = 'COPERNICUS' END SELECT IF (PRESENT(grid_mapping)) gridmp = TRIM(grid_mapping) @@ -3079,6 +3082,8 @@ SUBROUTINE def_topo_meta(diminfo,itopo_type,igrid_type,coordinates,grid_mapping, fr_land_topo_meta%long_name = 'fraction land due to GLOBE data' CASE(topo_merit) fr_land_topo_meta%long_name = 'fraction land due to MERIT data' + CASE(topo_copernicus) + fr_land_topo_meta%long_name = 'fraction land due to COPERNICUS data' END SELECT fr_land_topo_meta%shortName = 'FR_LAND' fr_land_topo_meta%stepType = 'instant'