Skip to content

Enable Copernicus DEM #427

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
ChristianSteger wants to merge 5 commits into
base: master
Choose a base branch
from
Open

Enable Copernicus DEM #427

Show file tree
Hide file tree
Changes from 4 commits
Commits
Show all changes
5 commits
Select commit Hold shift + click to select a range
2f6a865
Copernicus DEM implemented analogous to MERIT
Jun 30, 2025
0bbb128
Typo corrected, unnecessary/incorrect comments removed
Jun 30, 2025
2ac0aac
Run script added to test Copernicus DEM input
Jul 1, 2025
b43947e
Merge branch 'master' into copernicus_dem
stelliom Jul 22, 2025
1887034
Pull request review of 'stelliom' adressed (except test for Copernicu…
Sep 8, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
666 changes: 666 additions & 0 deletions run_scripts/extpar_mch.copernicus.sh
View file
Open in desktop

Large diffs are not rendered by default.

25 changes: 21 additions & 4 deletions src/extpar_topo_to_buffer.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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, &
Expand All @@ -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, &
Expand Down Expand Up @@ -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)
Expand Down
6 changes: 6 additions & 0 deletions src/mo_agg_sgsl.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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

Expand Down Expand Up @@ -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)
Copy link
Contributor

@stelliom stelliom Jul 23, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
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)

Also, could you do me a favor and apply the same changes to the cases above? 🙂

Copy link
Collaborator Author

@ChristianSteger ChristianSteger Sep 8, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Done!

ENDIF
END SELECT
!$OMP END CRITICAL
ENDIF
Expand Down
23 changes: 23 additions & 0 deletions src/mo_agg_topo_cosmo.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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, &
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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)
Copy link
Contributor

@stelliom stelliom Jul 23, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Unless you have a reason to add this code here I think we could just skip it, since COSMO is not really used anymore and I will probably start to remove COSMO related code at some point.

Copy link
Collaborator Author

@ChristianSteger ChristianSteger Sep 8, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Done - I removed the implementation of Copernicus DEM for COSMO.

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
Expand Down
10 changes: 7 additions & 3 deletions src/mo_agg_topo_icon.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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,&
Expand Down Expand Up @@ -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
Expand Down
11 changes: 10 additions & 1 deletion src/mo_extpar_output_nc.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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

Expand Down Expand Up @@ -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


Expand Down Expand Up @@ -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)
Expand Down
43 changes: 26 additions & 17 deletions src/mo_preproc_for_sgsl.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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, &
Expand Down Expand Up @@ -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)
Expand Down Expand Up @@ -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__)

Expand All @@ -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
Expand Down Expand Up @@ -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__)

Expand Down
Loading
Loading