Skip to content

Commit 29d2433

Browse files
author
Christian Steger
committed
New terrain horizon algorithm implemented
1 parent 496e175 commit 29d2433

File tree

2 files changed

+946
-47
lines changed

2 files changed

+946
-47
lines changed

‎src/extpar_topo_to_buffer.f90‎

Lines changed: 160 additions & 47 deletions
Original file line numberDiff line numberDiff line change
@@ -6,7 +6,7 @@
66
! V1_0 2010/12/21 Hermann Asensio
77
! Initial release
88
! V1_1 2011/01/20 Hermann Asensio
9-
! small bug fixes accroding to Fortran compiler warnings
9+
! small bug fixes according to Fortran compiler warnings
1010
! V1_2 2011/03/25 Hermann Asensio
1111
! update to support ICON refinement grids
1212
! V1_4 2011/04/21 Anne Roches
@@ -44,6 +44,8 @@
4444
!> \author Hermann Asensio
4545
PROGRAM extpar_topo_to_buffer
4646

47+
USE, INTRINSIC :: iso_c_binding
48+
4749
USE mo_logging
4850
USE info_extpar, ONLY: info_print
4951
USE mo_kind, ONLY: wp, i4
@@ -59,7 +61,8 @@ PROGRAM extpar_topo_to_buffer
5961

6062
USE mo_cosmo_grid, ONLY: COSMO_grid
6163

62-
USE mo_icon_grid_data, ONLY: ICON_grid
64+
USE mo_icon_grid_data, ONLY: ICON_grid, &
65+
& icon_grid_region
6366

6467
USE mo_io_units, ONLY: filename_max
6568

@@ -136,51 +139,92 @@ PROGRAM extpar_topo_to_buffer
136139

137140
IMPLICIT NONE
138141

139-
CHARACTER (len=filename_max) :: netcdf_filename, &
140-
& namelist_grid_def, &
141-
& namelist_topo_data_input, & !< file with input namelist with GLOBE data information
142-
& namelist_scale_sep_data_input, &!< file with input namelist with scale separated data information
143-
& namelist_oro_smooth, & !< file with orography smoothing information (switches)
144-
& namelist_lrad, & !< file with opo information (switches)
145-
& topo_files(1:max_tiles), & !< filenames globe raw data
146-
& sgsl_files(1:max_tiles), & !< filenames subgrid-slope
147-
& orography_buffer_file, & !< name for orography buffer file
148-
& orography_output_file, & !< name for orography output file
149-
& sgsl_output_file, & !< name for sgsl output file
150-
& raw_data_orography_path, & !< path to raw data
151-
& raw_data_scale_sep_orography_path, & !< path to raw data
152-
& scale_sep_files(1:max_tiles) !< filenames globe raw data
142+
CHARACTER (len=filename_max) :: netcdf_filename, &
143+
& namelist_grid_def, &
144+
& namelist_topo_data_input, & !< file with input namelist with GLOBE data information
145+
& namelist_scale_sep_data_input, &!< file with input namelist with scale separated data information
146+
& namelist_oro_smooth, & !< file with orography smoothing information (switches)
147+
& namelist_lrad, & !< file with opo information (switches)
148+
& topo_files(1:max_tiles), & !< filenames globe raw data
149+
& sgsl_files(1:max_tiles), & !< filenames subgrid-slope
150+
& orography_buffer_file, & !< name for orography buffer file
151+
& orography_output_file, & !< name for orography output file
152+
& sgsl_output_file, & !< name for sgsl output file
153+
& raw_data_orography_path, & !< path to raw data
154+
& raw_data_scale_sep_orography_path, & !< path to raw data
155+
& scale_sep_files(1:max_tiles) !< filenames globe raw data
153156

154-
REAL(KIND=wp) :: undefined !< value to indicate undefined grid elements
157+
REAL(KIND=wp) :: undefined !< value to indicate undefined grid elements
155158

156-
INTEGER (KIND=i4) :: &
157-
& k,ie,je,ke, &
158-
& igrid_type, & !< target grid type, 1 for ICON, 2 for COSMO, 3 for GME grid
159-
& ntiles_column, & !< number of tile columns in total domain
160-
& ntiles_row, & !< number of tile rows in total domain
161-
& ilow_pass_oro, &
162-
& numfilt_oro, &
163-
& ifill_valley, &
164-
& ilow_pass_xso, &
165-
& numfilt_xso
166-
167-
INTEGER (KIND=i4), ALLOCATABLE :: topo_startrow(:), & !< startrow indeces for each GLOBE tile
168-
& topo_endrow(:), & !< endrow indeces for each GLOBE tile
169-
& topo_startcolumn(:), & !< starcolumn indeces for each GLOBE tile
170-
& topo_endcolumn(:) !< endcolumn indeces for each GLOBE tile
171-
172-
REAL (KIND=wp) :: eps_filter, &
173-
& rfill_valley, &
174-
& rxso_mask
175-
176-
LOGICAL :: lsso_param, &
177-
& lcompute_sgsl=.FALSE., & !compute subgrid slope
178-
& lpreproc_oro = .FALSE., & !TRUE: preproc raw oro data FALSE: read directly from NetCDF
179-
& lscale_separation=.FALSE., &
180-
& lscale_file= .FALSE., &
181-
& lsubtract_mean_slope = .FALSE., &
182-
& lfilter_oro, &
183-
& lxso_first
159+
INTEGER (KIND=i4) :: &
160+
& k,ie,je,ke, &
161+
& igrid_type, & !< target grid type, 1 for ICON, 2 for COSMO, 3 for GME grid
162+
& ntiles_column, & !< number of tile columns in total domain
163+
& ntiles_row, & !< number of tile rows in total domain
164+
& ilow_pass_oro, &
165+
& numfilt_oro, &
166+
& ifill_valley, &
167+
& ilow_pass_xso, &
168+
& numfilt_xso
169+
170+
INTEGER (KIND=i4), ALLOCATABLE :: topo_startrow(:), & !< startrow indeces for each GLOBE tile
171+
& topo_endrow(:), & !< endrow indeces for each GLOBE tile
172+
& topo_startcolumn(:), & !< starcolumn indeces for each GLOBE tile
173+
& topo_endcolumn(:) !< endcolumn indeces for each GLOBE tile
174+
175+
INTEGER(c_int) :: num_cell_c, num_vertex_c, num_hori_c
176+
INTEGER(c_int) :: grid_type_c
177+
REAL (c_double) :: radius_c, ray_org_elev_c
178+
INTEGER(c_int) :: refine_factor_c, itype_scaling_c
179+
CHARACTER(KIND=c_char, len=2000) :: buffer_c
180+
INTEGER(c_int) :: buffer_len_c
181+
REAL(c_double), ALLOCATABLE :: clon_c(:), &
182+
& clat_c(:), &
183+
& hsurf_c(:), &
184+
& vlon_c(:), &
185+
& vlat_c(:), &
186+
& horizon_topo_c(:, :), &
187+
& skyview_topo_c(:)
188+
INTEGER(c_int), ALLOCATABLE :: cells_of_vertex_c(:, :)
189+
190+
REAL (KIND=wp) :: eps_filter, &
191+
& rfill_valley, &
192+
& rxso_mask
193+
194+
LOGICAL :: lsso_param, &
195+
& lcompute_sgsl=.FALSE., & !compute subgrid slope
196+
& lpreproc_oro = .FALSE., & !TRUE: preproc raw oro data FALSE: read directly from NetCDF
197+
& lscale_separation=.FALSE., &
198+
& lscale_file= .FALSE., &
199+
& lsubtract_mean_slope, &
200+
& lfilter_oro, &
201+
& lxso_first
202+
203+
INTERFACE
204+
SUBROUTINE horizon_svf_comp(clon_c, clat_c, hsurf_c, &
205+
& vlon_c, vlat_c, &
206+
& cells_of_vertex_c, &
207+
& horizon_topo_c, skyview_topo_c, &
208+
& num_cell_c, num_vertex_c, num_hori_c, &
209+
& grid_type_c, radius_c, &
210+
& ray_org_elev_c, refine_factor_c, &
211+
& itype_scaling_c, buffer_c, buffer_len_c) bind(C, name="horizon_svf_comp")
212+
USE iso_c_binding
213+
IMPLICIT NONE
214+
REAL(c_double), DIMENSION(*), INTENT(in) :: clon_c, clat_c
215+
REAL(c_double), DIMENSION(*), INTENT(in) :: hsurf_c
216+
REAL(c_double), DIMENSION(*), INTENT(in) :: vlon_c, vlat_c
217+
INTEGER(c_int), DIMENSION(*), INTENT(in) :: cells_of_vertex_c
218+
REAL(c_double), DIMENSION(*), INTENT(inout) :: horizon_topo_c
219+
REAL(c_double), DIMENSION(*), INTENT(inout) :: skyview_topo_c
220+
INTEGER(c_int), value :: num_cell_c, num_vertex_c, num_hori_c
221+
INTEGER(c_int), value :: grid_type_c
222+
REAL(c_double), value :: radius_c, ray_org_elev_c
223+
INTEGER(c_int), value :: refine_factor_c, itype_scaling_c
224+
CHARACTER(KIND=c_char), dimension(*) :: buffer_c
225+
INTEGER(c_int) :: buffer_len_c
226+
END SUBROUTINE horizon_svf_comp
227+
END INTERFACE
184228

185229
namelist_grid_def = 'INPUT_grid_org'
186230
namelist_scale_sep_data_input = 'INPUT_SCALE_SEP'
@@ -518,8 +562,77 @@ PROGRAM extpar_topo_to_buffer
518562
CALL compute_lradtopo(nhori,tg,hh_topo,slope_asp_topo,slope_ang_topo, &
519563
& horizon_topo,skyview_topo)
520564
ELSEIF ( igrid_type == igrid_icon ) THEN
521-
CALL lradtopo_icon(nhori, radius, min_circ_cov,tg, hh_topo, horizon_topo, &
522-
& skyview_topo, max_missing, itype_scaling)
565+
566+
! -----------------------------------------------------------------------
567+
! Old Fortran implementation
568+
! -----------------------------------------------------------------------
569+
570+
! CALL lradtopo_icon(nhori, radius, min_circ_cov,tg, hh_topo, horizon_topo, &
571+
! & skyview_topo, max_missing, itype_scaling)
572+
573+
! -----------------------------------------------------------------------
574+
! New C++ ray-tracing based implementation
575+
! -----------------------------------------------------------------------
576+
577+
! Cast non-contiguous arrays to C types
578+
ALLOCATE(clon_c(icon_grid_region%ncells))
579+
ALLOCATE(clat_c(icon_grid_region%ncells))
580+
DO k = 1, icon_grid_region%ncells
581+
clon_c(k) = REAL(icon_grid_region%cells%center(k)%lon, KIND=c_double)
582+
clat_c(k) = REAL(icon_grid_region%cells%center(k)%lat, KIND=c_double)
583+
END DO
584+
ALLOCATE(vlon_c(icon_grid_region%nverts))
585+
ALLOCATE(vlat_c(icon_grid_region%nverts))
586+
DO k = 1, icon_grid_region%nverts
587+
vlon_c(k) = REAL(icon_grid_region%verts%vertex(k)%lon, KIND=c_double)
588+
vlat_c(k) = REAL(icon_grid_region%verts%vertex(k)%lat, KIND=c_double)
589+
END DO
590+
591+
! Cast contiguous array to C type
592+
ALLOCATE(hsurf_c(icon_grid_region%ncells))
593+
hsurf_c = REAL(hh_topo(:,1,1), KIND=c_double)
594+
595+
! Adjust two-dimensional index arary: type and starting index
596+
ALLOCATE(cells_of_vertex_c(icon_grid_region%nverts, 6))
597+
cells_of_vertex_c = INT(icon_grid_region%verts%cell_index - 1, &
598+
& KIND=c_int)
599+
600+
! Cast scalar values to C types
601+
num_cell_c = INT(icon_grid_region%ncells, KIND=c_int)
602+
num_vertex_c = INT(icon_grid_region%nverts, KIND=c_int)
603+
num_hori_c = INT(nhori, KIND=c_int)
604+
radius_c = REAL(radius, KIND=c_double)
605+
itype_scaling_c = INT(itype_scaling, KIND=c_int)
606+
607+
! Constant settings
608+
grid_type_c = 1 ! triangle mesh buidling from ICON grid (0 or 1)
609+
ray_org_elev_c = 0.2 ! elevation of ray origin above ground level [m]
610+
refine_factor_c = 10 ! number of sub-sampling within azimuth sector
611+
612+
! Allocate output arrays
613+
ALLOCATE(horizon_topo_c(icon_grid_region%ncells, nhori))
614+
ALLOCATE(skyview_topo_c(icon_grid_region%ncells))
615+
616+
buffer_len_c = len(buffer_c)
617+
CALL horizon_svf_comp(clon_c, clat_c, hsurf_c, &
618+
& vlon_c, vlat_c, &
619+
& cells_of_vertex_c, &
620+
& horizon_topo_c, skyview_topo_c, &
621+
& num_cell_c, num_vertex_c, num_hori_c, &
622+
& grid_type_c, radius_c, &
623+
& ray_org_elev_c, refine_factor_c, &
624+
& itype_scaling_c, buffer_c, buffer_len_c)
625+
WRITE(message_text,*) buffer_c(3:buffer_len_c)
626+
CALL logging%info(message_text)
627+
628+
! Cast output to Fortran types
629+
horizon_topo(:,1,1,:) = REAL(horizon_topo_c)
630+
skyview_topo(:,1,1) = REAL(skyview_topo_c)
631+
632+
DEALLOCATE(clon_c, clat_c, hsurf_c, vlon_c, vlat_c)
633+
DEALLOCATE(cells_of_vertex_c)
634+
DEALLOCATE(horizon_topo_c, skyview_topo_c)
635+
523636
ENDIF
524637
ENDIF
525638

0 commit comments

Comments
 (0)