Improved algorithm to compute terrain horizon and Sky View Factor

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}',

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

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

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

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)