diff --git a/run_scripts/extpar_mch.aster.sh b/run_scripts/extpar_mch.aster.sh index fdec1c61..2045e8fe 100755 --- a/run_scripts/extpar_mch.aster.sh +++ b/run_scripts/extpar_mch.aster.sh @@ -3,20 +3,20 @@ #SBATCH --job-name="extpar" #SBATCH --nodes=1 #SBATCH --output="job.out" -#SBATCH --time=01:00:00 +#SBATCH --time=00:23:00 #SBATCH --partition=postproc -#SBATCH --cpus-per-task=12 +#SBATCH --cpus-per-task=64 #-------------------------------------------------------------------------------- # variables to define by user #-------------------------------------------------------------------------------- -# define model for which Extpar should run: c1, c2, i1, i2, i1_dev, i2_dev +# define model for which Extpar should run: c1, c2, i05, i1, i2, i1_dev, i2_dev model="i1" # Sandbox (make sure you have enough disk place at that location)! -sandboxdir=$SCRATCH/output_extpar/i1 +sandboxdir=$SCRATCH/ExtPar/output/${model} #-------------------------------------------------------------------------------- @@ -210,10 +210,40 @@ elif [[ $model == "i2_dev" ]]; then 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=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 == "i1" ]]; then #output file names - netcdf_output_filename="external_parameter_icon_grid_0001_R19B08_mch.nc" + netcdf_output_filename="extpar_icon_grid_0001_R19B08_mch.nc" # grid definition grid_dir="/oprusers/osm/opr/data/grid_descriptions/" @@ -292,11 +322,11 @@ 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 +binary_aot=extpar_aot_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.exe binary_soil=extpar_soil_to_buffer.exe binary_flake=extpar_flake_to_buffer.exe binary_consistency_check=extpar_consistency_check.exe @@ -402,6 +432,13 @@ input_alb = { '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': '', @@ -477,16 +514,6 @@ cat > INPUT_ICON_GRID << EOF_grid_icon / EOF_grid_icon -cat > INPUT_AOT << EOF_aot -&aerosol_raw_data - raw_data_aot_path='./', - raw_data_aot_filename='${raw_data_aot}' -/ -&aerosol_io_extpar - aot_buffer_file='${buffer_aot}', -/ -EOF_aot - cat > INPUT_LU << EOF_lu &lu_raw_data raw_data_lu_path='', @@ -564,7 +591,6 @@ isoil_data = 3, ldeep_soil = .FALSE., raw_data_soil_path='', raw_data_soil_filename='${raw_data_soil_HWSD}', -raw_data_deep_soil_filename='${raw_data_deep_soil}' / &soil_io_extpar soil_buffer_file='${buffer_soil}', diff --git a/src/extpar_topo_to_buffer.f90 b/src/extpar_topo_to_buffer.f90 index 00c0cbd1..8322aae9 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 @@ -29,6 +29,8 @@ ! Add namelist switch lfilter_topo ! V2_10 2018-02-19 Juergen Helmert ! lsubtract_mean_slope, ERA-I surface temp for land points +! V2_11 2024-12-10 Christian R. Steger +! introduction of the COPERNICUS topography raw data set for external parameters ! ! Code Description: ! Language: Fortran 2003. @@ -37,7 +39,8 @@ !> !! @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 GLOBE/ASTER/MERIT/COPERNICUS 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 +89,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 +110,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 +369,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..4538fa26 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 @@ -401,6 +402,11 @@ SUBROUTINE agg_sgsl_data_to_target_grid(sgsl_tiles_grid, & 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 ENDIF diff --git a/src/mo_agg_topo_cosmo.f90 b/src/mo_agg_topo_cosmo.f90 index f7069846..e56caa95 100644 --- a/src/mo_agg_topo_cosmo.f90 +++ b/src/mo_agg_topo_cosmo.f90 @@ -50,6 +50,7 @@ MODULE mo_agg_topo_cosmo & itopo_type, & & topo_gl, & & topo_merit, & + & topo_copernicus, & & topo_aster USE mo_topo_sso, ONLY: auxiliary_sso_parameter_cosmo, & @@ -271,6 +272,9 @@ SUBROUTINE agg_topo_data_to_target_grid_cosmo(topo_tiles_grid, & CASE(topo_merit) hh = undef_topo h_3rows = undef_topo + CASE(topo_copernicus) + hh = undef_topo + h_3rows = undef_topo END SELECT ! initialize some variables @@ -623,6 +627,25 @@ SUBROUTINE agg_topo_data_to_target_grid_cosmo(topo_tiles_grid, & h22(ie,je,ke) = h22(ie,je,ke) + dhdydy(i) ENDIF ENDIF + CASE(topo_copernicus) + IF (h_3rows(i,j_c) /= undef_topo) THEN + ndata(ie,je,ke) = ndata(ie,je,ke) + 1 + hh_target(ie,je,ke) = hh_target(ie,je,ke) + h_3rows(i,j_c) + hh2_target(ie,je,ke) = hh2_target(ie,je,ke) + (h_3rows(i,j_c) * h_3rows(i,j_c)) + IF (lscale_separation) THEN + hh_target_scale(ie,je,ke) = hh_target_scale(ie,je,ke) + & + & h_3rows_scale(i,j_c) + hh2_target_scale(ie,je,ke) = hh2_target_scale(ie,je,ke) + & + & (h_3rows_scale(i,j_c) * h_3rows_scale(i,j_c)) + hh_sqr_diff(ie,je,ke) = hh_sqr_diff(ie,je,ke)+ & + & (h_3rows(i,j_c) - h_3rows_scale(i,j_c))**2 + ENDIF + IF(lsso_param) THEN + h11(ie,je,ke) = h11(ie,je,ke) + dhdxdx(i) + h12(ie,je,ke) = h12(ie,je,ke) + dhdxdy(i) + h22(ie,je,ke) = h22(ie,je,ke) + dhdydy(i) + ENDIF + ENDIF END SELECT !$OMP END CRITICAL ENDIF diff --git a/src/mo_agg_topo_icon.f90 b/src/mo_agg_topo_icon.f90 index 6bbe306e..fa36c8ad 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, & !< there are 16/240 GLOBE/ASTER/MERIT/COPERNICUS tiles & 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 + & nr_tot, & !< total number of rows in GLOBE/ASTER/MERIT/COPERNICUS data & 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 0f832d89..54fab444 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 @@ -1717,6 +1717,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 @@ -1787,6 +1790,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..a396eb9e 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, & @@ -108,6 +108,8 @@ SUBROUTINE topo_subgrid_slope( infile, outfile, itopo_type) & 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, & @@ -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..e8a1bed5 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 ! if the user chooses GLOBE, ASTER, MERIT or COPERNICUS ntiles_column = columns ntiles_row = rows @@ -214,6 +225,11 @@ 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 @@ -246,18 +262,18 @@ SUBROUTINE fill_topo_data(raw_data_orography_path,topo_files, & 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. + 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 grid cell resolution (ASTER resolution: 1./3600. 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.) ! GLOBE resolution: 1./120. degrees (30 arc-seconds) + 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.) ! MERIT resolution: 1./1200. degrees + CASE(topo_copernicus) ! Copernicus topography: 324 tiles + CALL logging%info('COPERNICUS is used as topography') + half_gridp = 1./(3600.*2.) ! Copernicus resolution: 1./3600. degrees END SELECT DO i = 1,ntiles @@ -279,7 +295,7 @@ SUBROUTINE fill_topo_data(raw_data_orography_path,topo_files, & 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_lon_max(i) = REAL(NINT(tiles_lon_max(i) + half_gridp)) !< added, as the ASTER/GLOBE/MERIT/COPERNICUS 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)) END DO @@ -296,6 +312,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 +372,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 +417,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 +515,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 +557,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'