Correct grid shift in topography data (#426)

Co-authored-by: Christian Steger <csteger@balfrin-ln003.cscs.ch>
Co-authored-by: stelliom <67868694+stelliom@users.noreply.github.com>
Co-authored-by: Christian Steger <csteger@balfrin-ln002.cscs.ch>
Co-authored-by: Mikael Stellio <mikael.stellio@c2sm.ethz.ch>

Author: 4 people

src/mo_topo_data.f90

@@ -248,15 +248,15 @@ SUBROUTINE fill_topo_data(raw_data_orography_path,topo_files,  &
    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.
        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
+       half_gridp = 1._wp/7200._wp          ! 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
        CALL logging%info('GLOBE is used as topography')
-       half_gridp = 1./(120.*2.)                              ! GLOBE resolution is 1./120. degrees (30 arc-seconds)
+       half_gridp = 1._wp/240._wp                             ! GLOBE resolution is 1./120. degrees (30 arc-seconds)
 
      CASE(topo_merit)                                         ! ASTER topography, as it has 36 tiles at the moment.
        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
+       half_gridp = 1._wp/2400._wp           ! the resolution of the MERIT data is 1./1200. degrees as it is half a grid point
                                             ! it is additionally divided by 2
    END SELECT
 
@@ -278,10 +278,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), wp)
+         tiles_lon_max(i) = REAL(NINT(tiles_lon_max(i) + half_gridp), wp)
+         tiles_lat_min(i) = REAL(NINT(tiles_lat_min(i) - half_gridp), wp)
+         tiles_lat_max(i) = REAL(NINT(tiles_lat_max(i) + half_gridp), wp)
+       CASE(topo_merit) ! edges of raw data tiles do not align with integer coordinates
+         tiles_lon_min(i) = tiles_lon_min(i) - half_gridp
+         tiles_lon_max(i) = tiles_lon_max(i) + half_gridp
+         tiles_lat_min(i) = tiles_lat_min(i) - half_gridp
+         tiles_lat_max(i) = tiles_lat_max(i) + half_gridp
+     END SELECT
    END DO
 
    SELECT CASE(itopo_type)

src/mo_topo_routines.f90

@@ -285,13 +285,13 @@ SUBROUTINE det_topo_tiles_grid(topo_tiles_grid)
       dlat = -1.0_wp * (tiles_lat_max(k) - tiles_lat_min(k)) / REAL(tiles_nrows(k),wp)
 
       ! latitude from north to south, negative increment
-      topo_tiles_grid(k)%start_lon_reg  = tiles_lon_min(k) + 0.5 * dlon
-      topo_tiles_grid(k)%end_lon_reg    = tiles_lon_max(k) - 0.5 * dlon
+      topo_tiles_grid(k)%start_lon_reg  = tiles_lon_min(k) + 0.5_wp * dlon
+      topo_tiles_grid(k)%end_lon_reg    = tiles_lon_max(k) - 0.5_wp * dlon
 
       ! latitude from north to south, note the negative increment!
-      topo_tiles_grid(k)%start_lat_reg  = tiles_lat_max(k) + 0.5 * dlat
+      topo_tiles_grid(k)%start_lat_reg  = tiles_lat_max(k) + 0.5_wp * dlat
       ! latitude from north to south, note the negative increment!
-      topo_tiles_grid(k)%end_lat_reg    = tiles_lat_min(k) - 0.5 * dlat
+      topo_tiles_grid(k)%end_lat_reg    = tiles_lat_min(k) - 0.5_wp * dlat
 
       topo_tiles_grid(k)%dlon_reg = dlon
       topo_tiles_grid(k)%dlat_reg = dlat
@@ -318,7 +318,7 @@ SUBROUTINE det_topo_grid(topo_grid)
 
         dlon = (aster_lon_max - aster_lon_min) / REAL(nc_tot,wp)
 
-        dlat = -1. * (aster_lat_max - aster_lat_min) / REAL(nr_tot,wp)
+        dlat = -1._wp * (aster_lat_max - aster_lat_min) / REAL(nr_tot,wp)
         ! latitude from north to south, negative increment
 
         topo_grid%start_lon_reg  =  aster_lon_min + 0.5_wp * dlon
@@ -331,7 +331,7 @@ SUBROUTINE det_topo_grid(topo_grid)
 
         dlon = (merit_lon_max - merit_lon_min) / REAL(nc_tot,wp)
 
-        dlat = -1. * (merit_lat_max - merit_lat_min) / REAL(nr_tot,wp)
+        dlat = -1._wp * (merit_lat_max - merit_lat_min) / REAL(nr_tot,wp)
 
          WRITE(message_text,*)'Latitude increment for MERIT data, dlat = ', dlat
          CALL logging%info(message_text)
@@ -347,7 +347,7 @@ SUBROUTINE det_topo_grid(topo_grid)
       CASE(topo_gl)
 
         dlon = 360._wp / REAL(nc_tot,wp)
-        dlat = -1. * 180._wp / REAL(nr_tot,wp)
+        dlat = -1._wp * 180._wp / REAL(nr_tot,wp)
         ! latitude from north to south, negative increment
 
         topo_grid%start_lon_reg  = -180._wp + 0.5_wp * dlon

test/testsuite/data/get_data.sh

@@ -14,38 +14,38 @@ fi
 # mch
 test -d mch || exit 1
 cd mch/c7_globe
-wget --quiet 'ftp://iacftp.ethz.ch/pub_read/stelliom/external_parameter_mch_c7_PR417.nc'
+wget --quiet 'ftp://iacftp.ethz.ch/pub_read/stelliom/external_parameter_mch_c7_PR426.nc'
 cd -
 
 test -d mch || exit 1
 cd mch/c1_aster
-wget --quiet 'ftp://iacftp.ethz.ch/pub_read/stelliom/external_parameter_mch_c1_PR417.nc'
+wget --quiet 'ftp://iacftp.ethz.ch/pub_read/stelliom/external_parameter_mch_c1_PR426.nc'
 cd -
 
 # clm
 test -d clm || exit 1
 cd clm/12km_globe
-wget --quiet 'ftp://iacftp.ethz.ch/pub_read/stelliom/external_parameter_12km_globe_PR417.nc'
+wget --quiet 'ftp://iacftp.ethz.ch/pub_read/stelliom/external_parameter_12km_globe_PR426.nc'
 cd -
 
 cd clm/ecoclimap_sg
 wget --quiet 'ftp://iacftp.ethz.ch/pub_read/stelliom/icon_grid_bolivia.nc'
-wget --quiet 'ftp://iacftp.ethz.ch/pub_read/stelliom/external_parameter_icon_eco_PR417.nc'
+wget --quiet 'ftp://iacftp.ethz.ch/pub_read/stelliom/external_parameter_icon_eco_PR426.nc'
 cd -
 
 # dwd
 test -d dwd || exit 1
 
 cd dwd/icon_d2
 wget --quiet 'ftp://iacftp.ethz.ch/pub_read/stelliom/icon_grid_DOM01.nc'
-wget --quiet 'ftp://iacftp.ethz.ch/pub_read/stelliom/external_parameter_icon_d2_PR417.nc'
+wget --quiet 'ftp://iacftp.ethz.ch/pub_read/stelliom/external_parameter_icon_d2_PR426.nc'
 cd -
 
 cd dwd/icon_ecci
 wget --quiet 'ftp://iacftp.ethz.ch/pub_read/stelliom/icon_grid_bolivia.nc'
 wget --quiet 'ftp://iacftp.ethz.ch/pub_read/stelliom/clim_t2m_icon_ecci.nc'
 wget --quiet 'ftp://iacftp.ethz.ch/pub_read/stelliom/clim_tsea_icon_ecci.nc'
-wget --quiet 'ftp://iacftp.ethz.ch/pub_read/stelliom/external_parameter_icon_ecci_PR417.nc'
+wget --quiet 'ftp://iacftp.ethz.ch/pub_read/stelliom/external_parameter_icon_ecci_PR426.nc'
 cd -
 
 # mpim
@@ -54,14 +54,14 @@ cd mpim/icon_r2b4
 wget --quiet 'http://icon-downloads.mpimet.mpg.de/grids/public/mpim/0013/icon_grid_0013_R02B04_G.nc'
 wget --quiet 'ftp://iacftp.ethz.ch/pub_read/stelliom/clim_t2m_icon_r2b4.nc'
 wget --quiet 'ftp://iacftp.ethz.ch/pub_read/stelliom/clim_tsea_icon_r2b4.nc'
-wget --quiet 'ftp://iacftp.ethz.ch/pub_read/stelliom/external_parameter_icon_mpim_PR417.nc'
+wget --quiet 'ftp://iacftp.ethz.ch/pub_read/stelliom/external_parameter_icon_mpim_PR426.nc'
 cd -
 
 # ecmwf
 test -d ecmwf || exit 1
 cd ecmwf/corine_icon
 wget --quiet 'ftp://iacftp.ethz.ch/pub_read/stelliom/corine/icon_grid_0099_R19B10.nc'
-wget --quiet 'ftp://iacftp.ethz.ch/pub_read/stelliom/corine/external_parameter_icon_corine_PR417.nc'
+wget --quiet 'ftp://iacftp.ethz.ch/pub_read/stelliom/corine/external_parameter_icon_corine_PR426.nc'
 wget --quiet 'ftp://iacftp.ethz.ch/pub_read/stelliom/corine/clim_tsea_corine.nc'
 wget --quiet 'ftp://iacftp.ethz.ch/pub_read/stelliom/corine/clim_t2m_corine.nc'
 cd -