source: palm/trunk/UTIL/inifor/src/inifor_grid.f90 @ 3866

Last change on this file since 3866 was 3866, checked in by eckhard, 2 years ago

inifor: Use PALM's working precision; improved error handling, coding style, and comments

  • Property svn:keywords set to Id
File size: 205.2 KB
Line 
1!> @file src/inifor_grid.f90
2!------------------------------------------------------------------------------!
3! This file is part of the PALM model system.
4!
5! PALM is free software: you can redistribute it and/or modify it under the
6! terms of the GNU General Public License as published by the Free Software
7! Foundation, either version 3 of the License, or (at your option) any later
8! version.
9!
10! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
11! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
12! A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
13!
14! You should have received a copy of the GNU General Public License along with
15! PALM. If not, see <http://www.gnu.org/licenses/>.
16!
17! Copyright 2017-2019 Leibniz Universitaet Hannover
18! Copyright 2017-2019 Deutscher Wetterdienst Offenbach
19!------------------------------------------------------------------------------!
20!
21! Current revisions:
22! -----------------
23!
24!
25! Former revisions:
26! -----------------
27! $Id: inifor_grid.f90 3866 2019-04-05 14:25:01Z eckhard $
28! Use PALM's working precision
29! Catch errors while reading namelists
30! Improved coding style
31!
32!
33! 3802 2019-03-17 13:33:42Z raasch
34! unused variable removed
35!
36! 3801 2019-03-15 17:14:25Z eckhard
37! Read COSMO rotated pole from HHL file
38! Check for PALM and averaging domains extending outside COSMO domain
39!
40! 3785 2019-03-06 10:41:14Z eckhard
41! Assigned names to averaging grids
42! Improved variable naming and minor clean-up
43!
44!
45! 3765 2019-02-26 13:45:46Z eckhard
46! Removed dependency on radiation input files
47!
48!
49! 3716 2019-02-05 17:02:38Z eckhard
50! Removed dependency on soilmoisture input files
51!
52!
53! 3680 2019-01-18 14:54:12Z knoop
54! bugfix: Avoid empty averaging regions for small PALM domains
55! change: moved get_input_file_list() to io module
56!
57!
58! 3618 2018-12-10 13:25:22Z eckhard
59! Prefixed all INIFOR modules with inifor_, removed unused variables
60!
61!
62! 3615 2018-12-10 07:21:03Z raasch
63! bugfix: abort replaced by inifor_abort
64!
65! 3614 2018-12-10 07:05:46Z raasch
66! unused variables removed
67!
68! 3613 2018-12-07 18:20:37Z eckhard
69! Average initial profiles only over PALM domain region, not the
70!     geostrophic-wind averaging region
71! Fix unintended modification of COSMO heights
72! Fix commenting out of setup_volumentric setting
73! Moved version output to setup_parameters()
74!
75!
76! 3557 2018-11-22 16:01:22Z eckhard
77! Updated documentation
78!
79!
80! 3537 2018-11-20 10:53:14Z eckhard
81! Read COSMO domain extents and soil depths from input files
82! Report averaging mode and debugging mode in log
83!
84!
85! 3456 2018-10-30 14:29:54Z eckhard
86! Remove surface forcing from netCDF output (radiation + precipitation)
87! NetCDf output of internal arrays only with --debug option
88!
89!
90! 3447 2018-10-29 15:52:54Z eckhard
91! Renamed source files for compatibilty with PALM build system
92!
93!
94! 3395 2018-10-22 17:32:49Z eckhard
95! Added computation of geostrophic winds form COSMO pressure fields
96! Introduced averaging grids and internal 'output' variables for computation of
97!     geostrophic and large-scale forcings
98! Merged qv, uv and vg output variables and 'temperature' IO group into new
99!     'thermodynamics' IP group
100! Cleaned up terminal output, show some messages only with --debug option
101! Improved variable naming and readability
102!
103!
104! 3183 2018-07-27 14:25:55Z suehring
105! Introduced new PALM grid stretching
106! Updated variable names and metadata for PIDS v1.9 compatibility
107! Better compatibility with older Intel compilers:
108! - avoiding implicit array allocation with new get_netcdf_variable()
109!   subroutine instead of function
110! Improved command line interface:
111! - Produce forcing variables when '-mode profile' is being used
112! - Renamend initial-condition mode variable 'mode' to 'ic_mode'
113! - Improved handling of the start date string
114! Removed unnecessary variables and routines
115!
116!
117! 3182 2018-07-27 13:36:03Z suehring
118! Initial revision
119!
120!
121!
122! Authors:
123! --------
124!> @author Eckhard Kadasch (Deutscher Wetterdienst, Offenbach)
125!
126!------------------------------------------------------------------------------!
127! Description:
128! ------------
129!> The grid module contains all variables and routines to specify and work with
130!> the numerical grids in INIFOR. By convention, all angles are stored in
131!> radians.
132!------------------------------------------------------------------------------!
133
134#if defined ( __netcdf )
135 MODULE inifor_grid
136
137    USE inifor_control
138    USE inifor_defs,                                                           &
139        ONLY:  DATE, EARTH_RADIUS, TO_RADIANS, TO_DEGREES, PI,                 &
140               SNAME, LNAME, PATH, FORCING_STEP, WATER_ID, FILL_ITERATIONS,    &
141               BETA, P_SL, T_SL, BETA, RD, RV, G, P_REF, RD_PALM, CP_PALM,     &
142               RHO_L, OMEGA, HECTO, wp, iwp
143    USE inifor_io,                                                             &
144        ONLY:  get_cosmo_grid, get_netcdf_attribute, get_netcdf_dim_vector,    &
145               get_netcdf_variable, parse_command_line_arguments,              &
146               get_input_file_list, validate_config
147    USE inifor_transform,                                                      &
148        ONLY:  average_2d, rotate_to_cosmo, find_horizontal_neighbours,        &
149               compute_horizontal_interp_weights,                              &
150               find_vertical_neighbours_and_weights_interp,                    &
151               find_vertical_neighbours_and_weights_average, interpolate_2d,   &
152               gamma_from_hemisphere, phic_to_phin, lamc_to_lamn,              &
153               project, centre_velocities, phi2phirot, rla2rlarot, uv2uvrot,   &
154               phirot2phi, rlarot2rla
155    USE inifor_types
156    USE inifor_util
157    USE netcdf,                                                                &
158        ONLY:  NF90_MAX_NAME, NF90_MAX_VAR_DIMS
159   
160    IMPLICIT NONE
161   
162    SAVE
163   
164    REAL(wp) ::  averaging_angle   = 0.0_wp       !< latitudal and longitudal width of averaging regions [rad]
165    REAL(wp) ::  averaging_width_ns = 0.0_wp       !< longitudal width of averaging regions [m]
166    REAL(wp) ::  averaging_width_ew = 0.0_wp       !< latitudal width of averaging regions [m]
167    REAL(wp) ::  phi_equat         = 0.0_wp       !< latitude of rotated equator of COSMO-DE grid [rad]
168    REAL(wp) ::  phi_n             = 0.0_wp       !< latitude of rotated pole of COSMO-DE grid [rad]
169    REAL(wp) ::  lambda_n          = 0.0_wp       !< longitude of rotaded pole of COSMO-DE grid [rad]
170    REAL(wp) ::  phi_c             = 0.0_wp       !< rotated-grid latitude of the center of the PALM domain [rad]
171    REAL(wp) ::  lambda_c          = 0.0_wp       !< rotated-grid longitude of the centre of the PALM domain [rad]
172    REAL(wp) ::  phi_cn            = 0.0_wp       !< latitude of the rotated pole relative to the COSMO-DE grid [rad]
173    REAL(wp) ::  lambda_cn         = 0.0_wp       !< longitude of the rotated pole relative to the COSMO-DE grid [rad]
174    REAL(wp) ::  lam_centre        = 0.0_wp       !< longitude of the PLAM domain centre in the source (COSMO rotated-pole) system [rad]
175    REAL(wp) ::  phi_centre        = 0.0_wp       !< latitude of the PLAM domain centre in the source (COSMO rotated-pole) system [rad]
176    REAL(wp) ::  lam_east          = 0.0_wp       !< longitude of the east central-averaging-domain boundary in the source (COSMO rotated-pole) system [rad]
177    REAL(wp) ::  lam_west          = 0.0_wp       !< longitude of the west central-averaging-domain boundary in the source (COSMO rotated-pole) system [rad]
178    REAL(wp) ::  phi_north         = 0.0_wp       !< latitude of the north central-averaging-domain boundary in the source (COSMO rotated-pole) system [rad]
179    REAL(wp) ::  phi_south         = 0.0_wp       !< latitude of the south central-averaging-domain boundary in the source (COSMO rotated-pole) system [rad]
180    REAL(wp) ::  gam               = 0.0_wp       !< angle for working around phirot2phi/rlarot2rla bug
181    REAL(wp) ::  dx                = 0.0_wp       !< PALM-4U grid spacing in x direction [m]
182    REAL(wp) ::  dy                = 0.0_wp       !< PALM-4U grid spacing in y direction [m]
183    REAL(wp) ::  dz(10)            = -1.0_wp      !< PALM-4U grid spacing in z direction [m]
184    REAL(wp) ::  dz_max            = 1000.0_wp    !< maximum vertical grid spacing [m]
185    REAL(wp) ::  dz_stretch_factor = 1.08_wp      !< factor for vertical grid stretching [m]
186    REAL(wp) ::  dz_stretch_level  = -9999999.9_wp!< height above which the vertical grid will be stretched [m]
187    REAL(wp) ::  dz_stretch_level_start(9) = -9999999.9_wp !< namelist parameter
188    REAL(wp) ::  dz_stretch_level_end(9) = 9999999.9_wp !< namelist parameter
189    REAL(wp) ::  dz_stretch_factor_array(9) = 1.08_wp !< namelist parameter
190    REAL(wp) ::  dxi               = 0.0_wp       !< inverse PALM-4U grid spacing in x direction [m^-1]
191    REAL(wp) ::  dyi               = 0.0_wp       !< inverse PALM-4U grid spacing in y direction [m^-1]
192    REAL(wp) ::  dzi               = 0.0_wp       !< inverse PALM-4U grid spacing in z direction [m^-1]
193    REAL(wp) ::  f3                = 0.0_wp       !< Coriolis parameter
194    REAL(wp) ::  lx                = 0.0_wp       !< PALM-4U domain size in x direction [m]
195    REAL(wp) ::  ly                = 0.0_wp       !< PALM-4U domain size in y direction [m]
196    REAL(wp) ::  p0                = 0.0_wp       !< PALM-4U surface pressure, at z0 [Pa]
197    REAL(wp) ::  x0                = 0.0_wp       !< x coordinate of PALM-4U Earth tangent [m]
198    REAL(wp) ::  y0                = 0.0_wp       !< y coordinate of PALM-4U Earth tangent [m]
199    REAL(wp) ::  z0                = 0.0_wp       !< Elevation of the PALM-4U domain above sea level [m]
200    REAL(wp) ::  z_top             = 0.0_wp       !< height of the scalar top boundary [m]
201    REAL(wp) ::  zw_top            = 0.0_wp       !< height of the vertical velocity top boundary [m]
202    REAL(wp) ::  lonmin_cosmo      = 0.0_wp       !< Minimunm longitude of COSMO-DE's rotated-pole grid [COSMO rotated-pole rad]
203    REAL(wp) ::  lonmax_cosmo      = 0.0_wp       !< Maximum longitude of COSMO-DE's rotated-pole grid [COSMO rotated-pole rad]
204    REAL(wp) ::  latmin_cosmo      = 0.0_wp       !< Minimunm latitude of COSMO-DE's rotated-pole grid [COSMO rotated-pole rad]
205    REAL(wp) ::  latmax_cosmo      = 0.0_wp       !< Maximum latitude of COSMO-DE's rotated-pole grid [COSMO rotated-pole rad]
206    REAL(wp) ::  lonmin_palm       = 0.0_wp       !< Minimunm longitude of PALM grid [COSMO rotated-pole rad]
207    REAL(wp) ::  lonmax_palm       = 0.0_wp       !< Maximum longitude of PALM grid [COSMO rotated-pole rad]
208    REAL(wp) ::  latmin_palm       = 0.0_wp       !< Minimunm latitude of PALM grid [COSMO rotated-pole rad]
209    REAL(wp) ::  latmax_palm       = 0.0_wp       !< Maximum latitude of PALM grid [COSMO rotated-pole rad]
210    REAL(wp) ::  lonmin_tot        = 0.0_wp       !< Minimunm longitude of required COSMO data [COSMO rotated-pole rad]
211    REAL(wp) ::  lonmax_tot        = 0.0_wp       !< Maximum longitude of required COSMO data [COSMO rotated-pole rad]
212    REAL(wp) ::  latmin_tot        = 0.0_wp       !< Minimunm latitude of required COSMO data [COSMO rotated-pole rad]
213    REAL(wp) ::  latmax_tot        = 0.0_wp       !< Maximum latitude of required COSMO data [COSMO rotated-pole rad]
214    REAL(wp) ::  latitude          = 0.0_wp       !< geographical latitude of the PALM-4U origin, from inipar namelist [deg]
215    REAL(wp) ::  longitude         = 0.0_wp       !< geographical longitude of the PALM-4U origin, from inipar namelist [deg]
216    REAL(wp) ::  origin_lat        = 0.0_wp       !< geographical latitude of the PALM-4U origin, from static driver netCDF file [deg]
217    REAL(wp) ::  origin_lon        = 0.0_wp       !< geographical longitude of the PALM-4U origin, from static driver netCDF file [deg]
218    REAL(wp) ::  rotation_angle    = 0.0_wp       !< clockwise angle the PALM-4U north is rotated away from geographical north [deg]
219    REAL(wp) ::  end_time          = 0.0_wp       !< PALM-4U simulation time [s]
220
221    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  hhl             !< heights of half layers (cell faces) above sea level in COSMO-DE, read in from external file
222    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  hfl             !< heights of full layers (cell centres) above sea level in COSMO-DE, computed as arithmetic average of hhl
223    REAL(wp), DIMENSION(:), ALLOCATABLE, TARGET     ::  depths          !< COSMO-DE's TERRA-ML soil layer depths
224    REAL(wp), DIMENSION(:), ALLOCATABLE, TARGET     ::  d_depth         !< COSMO-DE's TERRA-ML soil layer thicknesses
225    REAL(wp), DIMENSION(:), ALLOCATABLE, TARGET     ::  d_depth_rho_inv !< inverted soil water mass
226    REAL(wp), DIMENSION(:), ALLOCATABLE, TARGET     ::  rlon            !< longitudes of COSMO-DE's rotated-pole grid
227    REAL(wp), DIMENSION(:), ALLOCATABLE, TARGET     ::  rlat            !< latitudes of COSMO-DE's rotated-pole grid
228    REAL(wp), DIMENSION(:), ALLOCATABLE, TARGET     ::  time            !< output times
229    REAL(wp), DIMENSION(:), ALLOCATABLE, TARGET     ::  x               !< base palm grid x coordinate vector pointed to by grid_definitions
230    REAL(wp), DIMENSION(:), ALLOCATABLE, TARGET     ::  xu              !< base palm grid xu coordinate vector pointed to by grid_definitions
231    REAL(wp), DIMENSION(:), ALLOCATABLE, TARGET     ::  y               !< base palm grid y coordinate vector pointed to by grid_definitions
232    REAL(wp), DIMENSION(:), ALLOCATABLE, TARGET     ::  yv              !< base palm grid yv coordinate vector pointed to by grid_definitions
233    REAL(wp), DIMENSION(:), ALLOCATABLE, TARGET     ::  z_column        !< base palm grid z coordinate vector including the top boundary coordinate (entire column)
234    REAL(wp), DIMENSION(:), ALLOCATABLE, TARGET     ::  zw_column       !< base palm grid zw coordinate vector including the top boundary coordinate (entire column)
235    REAL(wp), DIMENSION(:), ALLOCATABLE, TARGET     ::  z               !< base palm grid z coordinate vector pointed to by grid_definitions
236    REAL(wp), DIMENSION(:), ALLOCATABLE, TARGET     ::  zw              !< base palm grid zw coordinate vector pointed to by grid_definitions
237
238    INTEGER(iwp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  soiltyp     !< COSMO-DE soil type map
239    INTEGER ::  dz_stretch_level_end_index(9)               !< vertical grid level index until which the vertical grid spacing is stretched
240    INTEGER ::  dz_stretch_level_start_index(9)             !< vertical grid level index above which the vertical grid spacing is stretched
241    INTEGER ::  iostat !< return status of READ statement
242    INTEGER ::  nt    !< number of output time steps
243    INTEGER ::  nx    !< number of PALM-4U grid points in x direction
244    INTEGER ::  ny    !< number of PALM-4U grid points in y direction
245    INTEGER ::  nz    !< number of PALM-4U grid points in z direction
246    INTEGER ::  nlon  !< number of longitudal points in target grid (COSMO-DE)
247    INTEGER ::  nlat  !< number of latitudal points in target grid (COSMO-DE)
248    INTEGER ::  nlev  !< number of levels in target grid (COSMO-DE)
249    INTEGER ::  ndepths !< number of COSMO-DE soil layers
250    INTEGER ::  start_hour_flow         !< start of flow forcing in number of hours relative to start_date
251    INTEGER ::  start_hour_soil         !< start of soil forcing in number of hours relative to start_date, typically equals start_hour_flow
252    INTEGER ::  start_hour_radiation    !< start of radiation forcing in number of hours relative to start_date, 0 to 2 hours before start_hour_flow to reconstruct hourly averages from one- to three hourly averages of the input data
253    INTEGER ::  start_hour_soilmoisture !< start of forcing for the soil moisture spin-up in number of hours relative to start_date, typically -672 (-4 weeks)
254    INTEGER ::  end_hour  !< simulation time in hours
255    INTEGER ::  end_hour_soilmoisture  !< end of soil moisture spin-up in hours relative to start_hour_flow
256    INTEGER ::  step_hour !< number of hours between forcing time steps
257
258    LOGICAL ::  init_variables_required       !< flag controlling whether init variables are to be processed
259    LOGICAL ::  boundary_variables_required   !< flag controlling whether boundary grids are to be allocated and boundary variables are to be computed
260    LOGICAL ::  ls_forcing_variables_required !< flag controlling whether large-scale forcing variables are to be computed
261    LOGICAL ::  surface_forcing_required      !< flag controlling whether surface forcing variables are to be computed
262    LOGICAL ::  palm_domain_outside_cosmo     !< indicates whether COSMO grid covers the PALM domain and the geostrophic averaging domains
263
264    TYPE(nc_var), ALLOCATABLE, TARGET ::  input_var_table(:)  !< table of input variables
265    TYPE(nc_var), ALLOCATABLE, TARGET ::  output_var_table(:) !< table of input variables
266    TYPE(nc_var) ::  cosmo_var                                !< COSMO dummy variable, used for reading HHL, rlon, rlat
267
268    TYPE(grid_definition), TARGET ::  palm_grid                       !< PALM-4U grid in the target system (COSMO-DE rotated-pole)
269    TYPE(grid_definition), TARGET ::  palm_intermediate               !< PALM-4U grid with coarse vertical grid wiht levels interpolated from COSMO-DE grid
270    TYPE(grid_definition), TARGET ::  cosmo_grid                      !< target system (COSMO-DE rotated-pole)
271    TYPE(grid_definition), TARGET ::  scalars_east_grid               !< grid for eastern scalar boundary condition
272    TYPE(grid_definition), TARGET ::  scalars_west_grid               !< grid for western scalar boundary condition
273    TYPE(grid_definition), TARGET ::  scalars_north_grid              !< grid for northern scalar boundary condition
274    TYPE(grid_definition), TARGET ::  scalars_south_grid              !< grid for southern scalar boundary condition
275    TYPE(grid_definition), TARGET ::  scalars_top_grid                !< grid for top scalar boundary condition
276    TYPE(grid_definition), TARGET ::  scalars_east_intermediate       !< intermediate grid for eastern scalar boundary condition
277    TYPE(grid_definition), TARGET ::  scalars_west_intermediate       !< intermediate grid for western scalar boundary condition
278    TYPE(grid_definition), TARGET ::  scalars_north_intermediate      !< intermediate grid for northern scalar boundary condition
279    TYPE(grid_definition), TARGET ::  scalars_south_intermediate      !< intermediate grid for southern scalar boundary condition
280    TYPE(grid_definition), TARGET ::  scalars_top_intermediate        !< intermediate grid for top scalar boundary condition
281    TYPE(grid_definition), TARGET ::  u_initial_grid                  !< grid for u initial condition
282    TYPE(grid_definition), TARGET ::  u_east_grid                     !< grid for eastern u boundary condition
283    TYPE(grid_definition), TARGET ::  u_west_grid                     !< grid for western u boundary condition
284    TYPE(grid_definition), TARGET ::  u_north_grid                    !< grid for northern u boundary condition
285    TYPE(grid_definition), TARGET ::  u_south_grid                    !< grid for southern u boundary condition
286    TYPE(grid_definition), TARGET ::  u_top_grid                      !< grid for top u boundary condition
287    TYPE(grid_definition), TARGET ::  u_initial_intermediate          !< intermediate grid for u initial condition
288    TYPE(grid_definition), TARGET ::  u_east_intermediate             !< intermediate grid for eastern u boundary condition
289    TYPE(grid_definition), TARGET ::  u_west_intermediate             !< intermediate grid for western u boundary condition
290    TYPE(grid_definition), TARGET ::  u_north_intermediate            !< intermediate grid for northern u boundary condition
291    TYPE(grid_definition), TARGET ::  u_south_intermediate            !< intermediate grid for southern u boundary condition
292    TYPE(grid_definition), TARGET ::  u_top_intermediate              !< intermediate grid for top u boundary condition
293    TYPE(grid_definition), TARGET ::  v_initial_grid                  !< grid for v initial condition
294    TYPE(grid_definition), TARGET ::  v_east_grid                     !< grid for eastern v boundary condition
295    TYPE(grid_definition), TARGET ::  v_west_grid                     !< grid for western v boundary condition
296    TYPE(grid_definition), TARGET ::  v_north_grid                    !< grid for northern v boundary condition
297    TYPE(grid_definition), TARGET ::  v_south_grid                    !< grid for southern v boundary condition
298    TYPE(grid_definition), TARGET ::  v_top_grid                      !< grid for top v boundary condition
299    TYPE(grid_definition), TARGET ::  v_initial_intermediate          !< intermediate grid for v initial condition
300    TYPE(grid_definition), TARGET ::  v_east_intermediate             !< intermediate grid for eastern v boundary condition
301    TYPE(grid_definition), TARGET ::  v_west_intermediate             !< intermediate grid for western v boundary condition
302    TYPE(grid_definition), TARGET ::  v_north_intermediate            !< intermediate grid for northern v boundary condition
303    TYPE(grid_definition), TARGET ::  v_south_intermediate            !< intermediate grid for southern v boundary condition
304    TYPE(grid_definition), TARGET ::  v_top_intermediate              !< intermediate grid for top v boundary condition
305    TYPE(grid_definition), TARGET ::  w_initial_grid                  !< grid for w initial condition
306    TYPE(grid_definition), TARGET ::  w_east_grid                     !< grid for eastern w boundary condition
307    TYPE(grid_definition), TARGET ::  w_west_grid                     !< grid for western w boundary condition
308    TYPE(grid_definition), TARGET ::  w_north_grid                    !< grid for northern w boundary condition
309    TYPE(grid_definition), TARGET ::  w_south_grid                    !< grid for southern w boundary condition
310    TYPE(grid_definition), TARGET ::  w_top_grid                      !< grid for top w boundary condition
311    TYPE(grid_definition), TARGET ::  w_initial_intermediate          !< intermediate grid for w initial condition
312    TYPE(grid_definition), TARGET ::  w_east_intermediate             !< intermediate grid for eastern w boundary condition
313    TYPE(grid_definition), TARGET ::  w_west_intermediate             !< intermediate grid for western w boundary condition
314    TYPE(grid_definition), TARGET ::  w_north_intermediate            !< intermediate grid for northern w boundary condition
315    TYPE(grid_definition), TARGET ::  w_south_intermediate            !< intermediate grid for southern w boundary condition
316    TYPE(grid_definition), TARGET ::  w_top_intermediate              !< intermediate grid for top w boundary condition
317    TYPE(grid_definition), TARGET ::  north_averaged_scalar_profile   !< grid of the northern geostrophic scalar averaging region
318    TYPE(grid_definition), TARGET ::  south_averaged_scalar_profile   !< grid of the southern geostrophic scalar averaging region
319    TYPE(grid_definition), TARGET ::  west_averaged_scalar_profile    !< grid of the western geostrophic scalar averaging region
320    TYPE(grid_definition), TARGET ::  east_averaged_scalar_profile    !< grid of the eastern geostrophic scalar averaging region
321    TYPE(grid_definition), TARGET ::  averaged_scalar_profile         !< grid of the central geostrophic scalar averaging region
322    TYPE(grid_definition), TARGET ::  averaged_w_profile              !< grid of the central geostrophic w-velocity averaging region
323    TYPE(grid_definition), TARGET ::  averaged_initial_scalar_profile !< averaging grid for initial scalar profiles
324    TYPE(grid_definition), TARGET ::  averaged_initial_w_profile      !< averaging grid for the initial w profile
325
326    TYPE(io_group), ALLOCATABLE, TARGET ::  io_group_list(:)  !< List of I/O groups, which group together output variables that share the same input variable
327 
328    NAMELIST /inipar/ nx, ny, nz, dx, dy, dz, longitude, latitude,             &
329                      dz_max, dz_stretch_factor, dz_stretch_level,             &
330                      dz_stretch_level_start, dz_stretch_level_end
331    NAMELIST /d3par/  end_time
332   
333    CHARACTER(LEN=LNAME) ::  nc_source_text = ''  !< Text describing the source of the output data, e.g. 'COSMO-DE analysis from ...'
334
335    CHARACTER(LEN=PATH), ALLOCATABLE, DIMENSION(:) ::  flow_files          !< list of atmospheric input files (<prefix>YYYYMMDDHH-flow.nc)
336    CHARACTER(LEN=PATH), ALLOCATABLE, DIMENSION(:) ::  soil_moisture_files !< list of precipitation input files (<prefix>YYYYMMDDHH-soilmoisture.nc)
337    CHARACTER(LEN=PATH), ALLOCATABLE, DIMENSION(:) ::  soil_files          !< list of soil input files (temperature, moisture, <prefix>YYYYMMDDHH-soil.nc)
338    CHARACTER(LEN=PATH), ALLOCATABLE, DIMENSION(:) ::  radiation_files     !< list of radiation input files (<prefix>YYYYMMDDHH-rad.nc)
339
340    CHARACTER(LEN=SNAME) ::  input_prefix         !< prefix of input files, e.g. 'laf' for COSMO-DE analyses
341    CHARACTER(LEN=SNAME) ::  flow_prefix          !< prefix of flow input files, e.g. 'laf' for COSMO-DE analyses
342    CHARACTER(LEN=SNAME) ::  soil_prefix          !< prefix of soil input files, e.g. 'laf' for COSMO-DE analyses
343    CHARACTER(LEN=SNAME) ::  radiation_prefix     !< prefix of radiation input files, e.g. 'laf' for COSMO-DE analyses
344    CHARACTER(LEN=SNAME) ::  soilmoisture_prefix  !< prefix of input files for soil moisture spin-up, e.g. 'laf' for COSMO-DE analyses
345    CHARACTER(LEN=SNAME) ::  flow_suffix          !< suffix of flow input files, e.g. 'flow'
346    CHARACTER(LEN=SNAME) ::  soil_suffix          !< suffix of soil input files, e.g. 'soil'
347    CHARACTER(LEN=SNAME) ::  radiation_suffix     !< suffix of radiation input files, e.g. 'radiation'
348    CHARACTER(LEN=SNAME) ::  soilmoisture_suffix  !< suffix of input files for soil moisture spin-up, e.g. 'soilmoisture'
349                         
350    TYPE(nc_file) ::  output_file !< metadata of the dynamic driver
351
352    TYPE(inifor_config) ::  cfg !< container of the INIFOR command-line configuration
353
354 CONTAINS
355   
356!------------------------------------------------------------------------------!
357! Description:
358! ------------
359!> This routine initializes INIFOR. This includes parsing command-line options,
360!> setting the names of the input and output files, reading INIFOR's namelist
361!> file as well as reading in and setting grid parameters for defining
362!> interplation grids later in setup_grids().
363!------------------------------------------------------------------------------!
364 SUBROUTINE setup_parameters()
365
366!
367!------------------------------------------------------------------------------
368! Section 1: Define default parameters
369!------------------------------------------------------------------------------
370    cfg%start_date = '2013072100'
371    end_hour = 2
372    start_hour_soil = -2
373    start_hour_soilmoisture = - (4 * 7 * 24) - 2
374
375!
376!-- Defaultmain centre (_c) of the PALM-4U grid in the geographical system (_g)
377    origin_lat = 52.325079_wp * TO_RADIANS ! south-west of Berlin, origin used for the Dec 2017 showcase simulation
378    origin_lon = 13.082744_wp * TO_RADIANS
379    cfg%z0 = 35.0_wp
380
381!
382!-- Default atmospheric parameters
383    cfg%ug = 0.0_wp
384    cfg%vg = 0.0_wp
385    cfg%p0 = P_SL
386
387!
388!-- Parameters for file names
389    start_hour_flow = 0
390    start_hour_soil = 0
391    start_hour_radiation = 0
392    start_hour_soilmoisture = start_hour_flow - 2
393    end_hour_soilmoisture = start_hour_flow
394    step_hour = FORCING_STEP
395
396    input_prefix = 'laf'
397    cfg%flow_prefix = input_prefix
398    cfg%input_prefix = input_prefix
399    cfg%soil_prefix = input_prefix
400    cfg%radiation_prefix = input_prefix
401    cfg%soilmoisture_prefix  = input_prefix
402
403    flow_suffix = '-flow'
404    soil_suffix = '-soil'
405    radiation_suffix = '-rad'
406    soilmoisture_suffix = '-soilmoisture'
407
408    cfg%debug = .FALSE.
409    cfg%averaging_angle = 2.0_wp
410!
411!------------------------------------------------------------------------------
412! Section 2: Read command-line arguments, namelist, and grid parameters
413!------------------------------------------------------------------------------
414
415!
416!-- Set default paths and modes
417    cfg%input_path         = './'
418    cfg%hhl_file           = ''
419    cfg%soiltyp_file       = ''
420    cfg%namelist_file      = './namelist'
421    cfg%static_driver_file = ''
422    cfg%output_file = './palm-4u-input.nc'
423    cfg%ic_mode = 'volume'
424    cfg%bc_mode = 'real'
425    cfg%averaging_mode = 'level'
426
427!
428!-- Overwrite defaults with user configuration
429    CALL parse_command_line_arguments( cfg )
430    CALL report('main_loop', 'Running INIFOR version ' // VERSION)
431
432    flow_prefix = TRIM(cfg%input_prefix)
433    radiation_prefix = TRIM(cfg%input_prefix)
434    soil_prefix = TRIM(cfg%input_prefix)
435    soilmoisture_prefix = TRIM(cfg%input_prefix)
436    IF (cfg%flow_prefix_is_set)  flow_prefix = TRIM(cfg%flow_prefix)
437    IF (cfg%radiation_prefix_is_set)  radiation_prefix = TRIM(cfg%radiation_prefix)
438    IF (cfg%soil_prefix_is_set)  soil_prefix = TRIM(cfg%soil_prefix)
439    IF (cfg%soilmoisture_prefix_is_set)  soilmoisture_prefix = TRIM(cfg%soilmoisture_prefix)
440
441    output_file%name = cfg%output_file
442    z0 = cfg%z0
443    p0 = cfg%p0
444
445    init_variables_required = .TRUE.
446    boundary_variables_required = TRIM( cfg%bc_mode ) == 'real'
447    ls_forcing_variables_required = TRIM( cfg%bc_mode ) == 'ideal'
448    surface_forcing_required = .FALSE.
449
450    IF ( ls_forcing_variables_required )  THEN
451       message = "Averaging of large-scale forcing profiles " //            &
452                 "has not been implemented, yet."
453       CALL inifor_abort('setup_parameters', message)
454    ENDIF
455
456!
457!-- Set default file paths, if not specified by user.
458    CALL normalize_path(cfg%input_path)
459    IF (TRIM(cfg%hhl_file) == '')  cfg%hhl_file = TRIM(cfg%input_path) // 'hhl.nc'
460    IF (TRIM(cfg%soiltyp_file) == '')  cfg%soiltyp_file = TRIM(cfg%input_path) // 'soil.nc'
461
462    CALL validate_config( cfg ) 
463
464    CALL report('setup_parameters', "initialization mode: " // TRIM(cfg%ic_mode))
465    CALL report('setup_parameters', "       forcing mode: " // TRIM(cfg%bc_mode))
466    CALL report('setup_parameters', "     averaging mode: " // TRIM(cfg%averaging_mode))
467    CALL report('setup_parameters', "    averaging angle: " // real_to_str(cfg%averaging_angle))
468    CALL report('setup_parameters', "          data path: " // TRIM(cfg%input_path))
469    CALL report('setup_parameters', "           hhl file: " // TRIM(cfg%hhl_file))
470    CALL report('setup_parameters', "       soiltyp file: " // TRIM(cfg%soiltyp_file))
471    CALL report('setup_parameters', "      namelist file: " // TRIM(cfg%namelist_file))
472    CALL report('setup_parameters', "   output data file: " // TRIM(output_file%name))
473    IF (cfg%debug )  CALL report('setup_parameters', "     debugging mode: enabled")
474
475    CALL log_runtime('time', 'init')
476!
477!-- Read in namelist parameters
478    OPEN(10, FILE=cfg%namelist_file, STATUS='old')
479    READ(10, NML=inipar, IOSTAT=iostat) ! nx, ny, nz, dx, dy, dz   
480    IF ( iostat > 0 )  THEN     
481       message = "Failed to read namelist 'inipar' from file '" //             &
482                 TRIM( cfg%namelist_file ) // "'. "
483       CALL inifor_abort( 'setup_parameters', message )
484       CLOSE(10)
485    ENDIF
486
487    READ(10, NML=d3par, IOSTAT=iostat)  ! end_time
488    IF ( iostat > 0 )  THEN
489       message = "Failed to read namelist 'd3par' from file '" //              &
490                 TRIM( cfg%namelist_file ) // "'. "
491       CALL inifor_abort( 'setup_parameters', message )
492       CLOSE(10)
493    ENDIF
494    CLOSE(10)
495   
496    CALL log_runtime('time', 'read')
497
498    end_hour = CEILING( end_time / 3600.0 * step_hour )
499
500!
501!-- Generate input file lists
502    CALL get_input_file_list(                                               &
503       cfg%start_date, start_hour_flow, end_hour, step_hour,              &
504       cfg%input_path, flow_prefix, flow_suffix, flow_files)
505    CALL get_input_file_list(                                               &
506       cfg%start_date, start_hour_soil, end_hour, step_hour,              &
507       cfg%input_path, soil_prefix, soil_suffix, soil_files)
508    CALL get_input_file_list(                                               &
509       cfg%start_date, start_hour_radiation, end_hour, step_hour,         &
510       cfg%input_path, radiation_prefix, radiation_suffix, radiation_files, nocheck=.TRUE.)
511    CALL get_input_file_list(                                               &
512       cfg%start_date, start_hour_soilmoisture, end_hour_soilmoisture, step_hour, &
513       cfg%input_path, soilmoisture_prefix, soilmoisture_suffix, soil_moisture_files, nocheck=.TRUE.)
514
515!
516!------------------------------------------------------------------------------
517! Section 3: Check for consistency
518!------------------------------------------------------------------------------
519
520!
521!------------------------------------------------------------------------------
522! Section 4: Compute additional parameters
523!------------------------------------------------------------------------------
524!------------------------------------------------------------------------------
525! Section 4.1: COSMO-DE parameters
526!------------------------------------------------------------------------------
527
528
529    CALL log_runtime('time', 'init')
530!
531!-- Read COSMO soil type map
532    cosmo_var%name = 'SOILTYP'
533    CALL get_netcdf_variable(cfg%soiltyp_file, cosmo_var, soiltyp)
534
535    message = 'Reading PALM-4U origin from'
536    IF (TRIM(cfg%static_driver_file) .NE. '')  THEN
537
538       origin_lon = get_netcdf_attribute(cfg%static_driver_file, 'origin_lon')
539       origin_lat = get_netcdf_attribute(cfg%static_driver_file, 'origin_lat')
540
541       message = TRIM(message) // " static driver file '"                   &
542                               // TRIM(cfg%static_driver_file) // "'"
543
544
545    ELSE
546
547       origin_lon = longitude
548       origin_lat = latitude
549
550       message = TRIM(message) // " namlist file '"                         &
551                               // TRIM(cfg%namelist_file) // "'"
552
553    ENDIF
554    origin_lon = origin_lon * TO_RADIANS
555    origin_lat = origin_lat * TO_RADIANS
556
557    CALL report('setup_parameters', message)
558
559    CALL log_runtime('time', 'read')
560
561    CALL get_cosmo_grid( cfg, soil_files(1), rlon, rlat, hhl, hfl, depths,  &
562                         d_depth, d_depth_rho_inv, phi_n, lambda_n,         &
563                         phi_equat,                                         &
564                         lonmin_cosmo, lonmax_cosmo,                        &
565                         latmin_cosmo, latmax_cosmo,                        &
566                         nlon, nlat, nlev, ndepths )
567
568
569!------------------------------------------------------------------------------
570! Section 4.2: PALM-4U parameters
571!------------------------------------------------------------------------------
572!
573!-- PALM-4U domain extents
574    lx = (nx+1) * dx
575    ly = (ny+1) * dy
576   
577!
578!-- PALM-4U point of Earth tangency
579    x0 = 0.0_wp
580    y0 = 0.0_wp
581
582!
583!-- time vector
584    nt = CEILING(end_time / (step_hour * 3600.0_wp)) + 1
585    ALLOCATE( time(nt) )
586    CALL linspace(0.0_wp, 3600.0_wp * (nt-1), time)
587    output_file%time => time
588    CALL log_runtime('time', 'init')
589
590!
591!-- Convert the PALM-4U origin coordinates to COSMO's rotated-pole grid
592    phi_c    = TO_RADIANS *                                                 &
593               phi2phirot( origin_lat * TO_DEGREES, origin_lon * TO_DEGREES,&
594                           phi_n * TO_DEGREES, lambda_n * TO_DEGREES )
595    lambda_c = TO_RADIANS *                                                 &
596               rla2rlarot( origin_lat * TO_DEGREES, origin_lon * TO_DEGREES,&
597                           phi_n * TO_DEGREES, lambda_n * TO_DEGREES,     &
598                           0.0_wp )
599
600!
601!-- Set gamma according to whether PALM domain is in the northern or southern
602!-- hemisphere of the COSMO rotated-pole system. Gamma assumes either the
603!-- value 0 or PI and is needed to work around around a bug in the
604!-- rotated-pole coordinate transformations.
605    gam = gamma_from_hemisphere(origin_lat, phi_equat)
606
607!
608!-- Compute the north pole of the rotated-pole grid centred at the PALM-4U
609!-- domain centre. The resulting (phi_cn, lambda_cn) are coordinates in
610!-- COSMO-DE's rotated-pole grid.
611    phi_cn    = phic_to_phin(phi_c) 
612    lambda_cn = lamc_to_lamn(phi_c, lambda_c) 
613
614    message =   "PALM-4U origin:" // NEW_LINE('') // &
615       "           lon (lambda) = " // &
616       TRIM(real_to_str_f(origin_lon * TO_DEGREES)) // " deg"// NEW_LINE(' ') //&
617       "           lat (phi   ) = " // &
618       TRIM(real_to_str_f(origin_lat * TO_DEGREES)) // " deg (geographical)" // NEW_LINE(' ') //&
619       "           lon (lambda) = " // &
620       TRIM(real_to_str_f(lambda_c * TO_DEGREES)) // " deg" // NEW_LINE(' ') // &
621       "           lat (phi   ) = " // &
622       TRIM(real_to_str_f(phi_c * TO_DEGREES)) // " deg (COSMO-DE rotated-pole)"
623    CALL report ('setup_parameters', message)
624
625    message = "North pole of the rotated COSMO-DE system:" // NEW_LINE(' ') // &
626       "           lon (lambda) = " // &
627       TRIM(real_to_str_f(lambda_n * TO_DEGREES)) // " deg" // NEW_LINE(' ') //&
628       "           lat (phi   ) = " // &
629       TRIM(real_to_str_f(phi_n * TO_DEGREES)) // " deg (geographical)"
630    CALL report ('setup_parameters', message)
631       
632    message = "North pole of the rotated palm system:" // NEW_LINE(' ') // &
633       "           lon (lambda) = " // &
634       TRIM(real_to_str_f(lambda_cn * TO_DEGREES)) // " deg" // NEW_LINE(' ') // &
635       "           lat (phi   ) = " // &
636       TRIM(real_to_str_f(phi_cn * TO_DEGREES)) // " deg (COSMO-DE rotated-pole)"
637    CALL report ('setup_parameters', message)
638
639    CALL log_runtime('time', 'comp')
640
641!------------------------------------------------------------------------------
642! Section 4.3: INIFOR averaging domains
643!------------------------------------------------------------------------------
644
645!
646!-- Compute coordiantes of the PALM centre in the source (COSMO) system
647    phi_centre = phirot2phi(                                                   &
648       phirot = project(0.5_wp*ly, y0, EARTH_RADIUS) * TO_DEGREES,             &
649       rlarot = project(0.5_wp*lx, x0, EARTH_RADIUS) * TO_DEGREES,             &
650       polphi = phi_cn * TO_DEGREES,                                           &
651       polgam = gam * TO_DEGREES                                               &
652    ) * TO_RADIANS
653
654    lam_centre = rlarot2rla(                                                   &
655       phirot = project(0.5_wp*ly, y0, EARTH_RADIUS) * TO_DEGREES,             &
656       rlarot = project(0.5_wp*lx, x0, EARTH_RADIUS) * TO_DEGREES,             &
657       polphi = phi_cn * TO_DEGREES, pollam = lambda_cn * TO_DEGREES,          &
658       polgam = gam * TO_DEGREES                                               &
659    ) * TO_RADIANS
660
661    message = "PALM-4U centre:" // NEW_LINE('') // &
662       "           lon (lambda) = " // &
663       TRIM(real_to_str_f(lam_centre * TO_DEGREES)) // " deg" // NEW_LINE(' ') // &
664       "           lat (phi   ) = " // &
665       TRIM(real_to_str_f(phi_centre * TO_DEGREES)) // " deg (COSMO-DE rotated-pole)"
666    CALL report( 'setup_parameters', message )
667
668!
669!-- Compute boundaries of the central averaging box
670    averaging_angle = cfg%averaging_angle * TO_RADIANS
671    lam_east = lam_centre + 0.5_wp * averaging_angle
672    lam_west = lam_centre - 0.5_wp * averaging_angle
673    phi_north = phi_centre + 0.5_wp * averaging_angle
674    phi_south = phi_centre - 0.5_wp * averaging_angle
675    averaging_width_ew = averaging_angle * COS(phi_centre) * EARTH_RADIUS
676    averaging_width_ns = averaging_angle * EARTH_RADIUS
677
678    lonmin_tot = MIN(lam_centre - averaging_angle, lonmin_palm)
679    lonmax_tot = MAX(lam_centre + averaging_angle, lonmax_palm)
680    latmin_tot = MIN(phi_centre - averaging_angle, latmin_palm)
681    latmax_tot = MAX(phi_centre + averaging_angle, latmax_palm)
682
683    palm_domain_outside_cosmo = ANY(                                           &
684       (/ lonmin_tot,   -lonmax_tot,   latmin_tot,   -latmax_tot/) .LT.        &
685       (/ lonmin_cosmo, -lonmax_cosmo, latmin_cosmo, -latmax_cosmo/)           &
686    )
687
688    IF ( palm_domain_outside_cosmo )  THEN
689       message = 'PALM domain or geostrophic averaging domains extend ' //     &
690                 'outside COSMO domain.'
691       CALL inifor_abort( 'setup_parameters', message )
692    ENDIF
693
694
695!
696!-- Coriolis parameter
697    f3 = 2.0_wp * OMEGA * SIN(                                                 &
698       TO_RADIANS*phirot2phi( phi_centre * TO_DEGREES, lam_centre * TO_DEGREES,&
699                              phi_n * TO_DEGREES,                              &
700                              gam * TO_DEGREES )                               &
701    )
702
703 END SUBROUTINE setup_parameters
704
705
706!------------------------------------------------------------------------------!
707! Description:
708! ------------
709!> Defines the COSMO, PALM-4U, PALM-4U boundary grids, in particular their
710!> coordinates and interpolation weights
711!------------------------------------------------------------------------------!
712 SUBROUTINE setup_grids()
713    CHARACTER ::  interp_mode
714
715!------------------------------------------------------------------------------
716! Section 0: Define base PALM-4U grid coordinate vectors
717!------------------------------------------------------------------------------
718!
719!-- palm x y z, we allocate the column to nz+1 in order to include the top
720!-- scalar boundary. The interpolation grids will be associated with
721!-- a shorter column that omits the top element.
722    ALLOCATE( x(0:nx), y(0:ny), z(1:nz), z_column(1:nz+1) )
723    CALL linspace(0.5_wp * dx, lx - 0.5_wp * dx, x)
724    CALL linspace(0.5_wp * dy, ly - 0.5_wp * dy, y)
725    CALL stretched_z(z_column, dz, dz_max=dz_max,                           &
726                     dz_stretch_factor=dz_stretch_factor,                   &
727                     dz_stretch_level=dz_stretch_level,                     &
728                     dz_stretch_level_start=dz_stretch_level_start,         &
729                     dz_stretch_level_end=dz_stretch_level_end,             &
730                     dz_stretch_factor_array=dz_stretch_factor_array)
731    z(1:nz) = z_column(1:nz)
732    z_top = z_column(nz+1)
733
734!
735!-- palm xu yv zw, compared to the scalar grid, velocity coordinates
736!-- contain one element less.
737    ALLOCATE( xu(1:nx),  yv(1:ny), zw(1:nz-1), zw_column(1:nz))
738    CALL linspace(dx, lx - dx, xu)
739    CALL linspace(dy, ly - dy, yv)
740    CALL midpoints(z_column, zw_column)
741    zw(1:nz-1) = zw_column(1:nz-1)
742    zw_top     = zw_column(nz)
743
744
745!------------------------------------------------------------------------------
746! Section 1: Define initialization and boundary grids
747!------------------------------------------------------------------------------
748    CALL init_grid_definition('palm', grid=palm_grid,                       &
749            xmin=0.0_wp, xmax=lx,                                           &
750            ymin=0.0_wp, ymax=ly,                                           &
751            x0=x0, y0=y0, z0=z0,                                            &
752            nx=nx, ny=ny, nz=nz, z=z, zw=zw, ic_mode=cfg%ic_mode)
753
754!
755!-- Subtracting 1 because arrays will be allocated with nlon + 1 elements.
756    CALL init_grid_definition('cosmo-de', grid=cosmo_grid,                  &
757            xmin=lonmin_cosmo, xmax=lonmax_cosmo,                           &
758            ymin=latmin_cosmo, ymax=latmax_cosmo,                           &
759            x0=x0, y0=y0, z0=0.0_wp,                                        &
760            nx=nlon-1, ny=nlat-1, nz=nlev-1)
761
762!
763!-- Define intermediate grid. This is the same as palm_grid except with a
764!-- much coarser vertical grid. The vertical levels are interpolated in each
765!-- PALM column from COSMO's secondary levels. The main levels are then
766!-- computed as the averages of the bounding secondary levels.
767    CALL init_grid_definition('palm intermediate', grid=palm_intermediate,  &
768            xmin=0.0_wp, xmax=lx,                                           &
769            ymin=0.0_wp, ymax=ly,                                           &
770            x0=x0, y0=y0, z0=z0,                                            &
771            nx=nx, ny=ny, nz=nlev-2)
772
773    CALL init_grid_definition('boundary', grid=u_initial_grid,              &
774            xmin = dx, xmax = lx - dx,                                      &
775            ymin = 0.5_wp * dy, ymax = ly - 0.5_wp * dy,                    &
776            x0=x0, y0=y0, z0 = z0,                                          &
777            nx = nx-1, ny = ny, nz = nz,                                    &
778            z=z, ic_mode=cfg%ic_mode)
779
780    CALL init_grid_definition('boundary', grid=v_initial_grid,              &
781            xmin = 0.5_wp * dx, xmax = lx - 0.5_wp * dx,                    &
782            ymin = dy, ymax = ly - dy,                                      &
783            x0=x0, y0=y0, z0 = z0,                                          &
784            nx = nx, ny = ny-1, nz = nz,                                    &
785            z=z, ic_mode=cfg%ic_mode)
786
787    CALL init_grid_definition('boundary', grid=w_initial_grid,              &
788            xmin = 0.5_wp * dx, xmax = lx - 0.5_wp * dx,                    &
789            ymin = 0.5_wp * dy, ymax = ly - 0.5_wp * dy,                    &
790            x0=x0, y0=y0, z0 = z0,                                          &
791            nx = nx, ny = ny, nz = nz-1,                                    &
792            z=zw, ic_mode=cfg%ic_mode)
793
794    CALL init_grid_definition('boundary intermediate', grid=u_initial_intermediate,      &
795            xmin = dx, xmax = lx - dx,                                      &
796            ymin = 0.5_wp * dy, ymax = ly - 0.5_wp * dy,                    &
797            x0=x0, y0=y0, z0 = z0,                                          &
798            nx = nx-1, ny = ny, nz = nlev - 2)
799
800    CALL init_grid_definition('boundary intermediate', grid=v_initial_intermediate,      &
801            xmin = 0.5_wp * dx, xmax = lx - 0.5_wp * dx,                    &
802            ymin = dy, ymax = ly - dy,                                      &
803            x0=x0, y0=y0, z0 = z0,                                          &
804            nx = nx, ny = ny-1, nz = nlev - 2)
805
806    CALL init_grid_definition('boundary intermediate', grid=w_initial_intermediate,      &
807            xmin = 0.5_wp * dx, xmax = lx - 0.5_wp * dx,                    &
808            ymin = 0.5_wp * dy, ymax = ly - 0.5_wp * dy,                    &
809            x0=x0, y0=y0, z0 = z0,                                          &
810            nx = nx, ny = ny, nz = nlev - 1)
811
812    IF (boundary_variables_required)  THEN
813!
814!------------------------------------------------------------------------------
815! Section 2: Define PALM-4U boundary grids
816!------------------------------------------------------------------------------
817       CALL init_grid_definition('boundary', grid=scalars_east_grid,           &
818               xmin = lx + 0.5_wp * dx, xmax = lx + 0.5_wp * dx,               &
819               ymin =  0.5_wp * dy, ymax = ly - 0.5_wp * dy,                   &
820               x0=x0, y0=y0, z0 = z0,                                          &
821               nx = 0, ny = ny, nz = nz, z=z)
822
823       CALL init_grid_definition('boundary', grid=scalars_west_grid,           &
824               xmin = -0.5_wp * dx, xmax = -0.5_wp * dx,                       &
825               ymin =  0.5_wp * dy, ymax = ly - 0.5_wp * dy,                   &
826               x0=x0, y0=y0, z0 = z0,                                          &
827               nx = 0, ny = ny, nz = nz, z=z)
828
829       CALL init_grid_definition('boundary', grid=scalars_north_grid,          &
830               xmin = 0.5_wp * dx, xmax = lx - 0.5_wp * dx,                    &
831               ymin = ly + 0.5_wp * dy, ymax = ly + 0.5_wp * dy,               &
832               x0=x0, y0=y0, z0 = z0,                                          &
833               nx = nx, ny = 0, nz = nz, z=z)
834
835       CALL init_grid_definition('boundary', grid=scalars_south_grid,          &
836               xmin =  0.5_wp * dx, xmax = lx - 0.5_wp * dx,                   &
837               ymin = -0.5_wp * dy, ymax = -0.5_wp * dy,                       &
838               x0=x0, y0=y0, z0 = z0,                                          &
839               nx = nx, ny = 0, nz = nz, z=z)
840
841       CALL init_grid_definition('boundary', grid=scalars_top_grid,            &
842               xmin =  0.5_wp * dx, xmax = lx - 0.5_wp * dx,                   &
843               ymin =  0.5_wp * dy, ymax = ly - 0.5_wp * dy,                   &
844               x0=x0, y0=y0, z0 = z0,                                          &
845               nx = nx, ny = ny, nz = 1, z=(/z_top/))
846
847       CALL init_grid_definition('boundary', grid=u_east_grid,                 &
848               xmin = lx, xmax = lx,                                           &
849               ymin = 0.5_wp * dy, ymax = ly - 0.5_wp * dy,                    &
850               x0=x0, y0=y0, z0 = z0,                                          &
851               nx = 0, ny = ny, nz = nz, z=z)
852
853       CALL init_grid_definition('boundary', grid=u_west_grid,                 &
854               xmin = 0.0_wp, xmax = 0.0_wp,                                   &
855               ymin = 0.5_wp * dy, ymax = ly - 0.5_wp * dy,                    &
856               x0=x0, y0=y0, z0 = z0,                                          &
857               nx = 0, ny = ny, nz = nz, z=z)
858
859       CALL init_grid_definition('boundary', grid=u_north_grid,                &
860               xmin = dx, xmax = lx - dx,                                      &
861               ymin = ly + 0.5_wp * dy, ymax = ly + 0.5_wp * dy,               &
862               x0=x0, y0=y0, z0 = z0,                                          &
863               nx = nx-1, ny = 0, nz = nz, z=z)
864   
865       CALL init_grid_definition('boundary', grid=u_south_grid,                &
866               xmin = dx, xmax = lx - dx,                                      &
867               ymin = -0.5_wp * dy, ymax = -0.5_wp * dy,                       &
868               x0=x0, y0=y0, z0 = z0,                                          &
869               nx = nx-1, ny = 0, nz = nz, z=z)
870
871       CALL init_grid_definition('boundary', grid=u_top_grid,                  &
872               xmin = dx, xmax = lx - dx,                                      &
873               ymin = 0.5_wp * dy, ymax = ly - 0.5_wp * dy,                    &
874               x0=x0, y0=y0, z0 = z0,                                          &
875               nx = nx-1, ny = ny, nz = 1, z=(/z_top/))
876
877       CALL init_grid_definition('boundary', grid=v_east_grid,                 &
878               xmin = lx + 0.5_wp * dx, xmax = lx + 0.5_wp * dx,               &
879               ymin = dy, ymax = ly - dy,                                      &
880               x0=x0, y0=y0, z0 = z0,                                          &
881               nx = 0, ny = ny-1, nz = nz, z=z)
882
883       CALL init_grid_definition('boundary', grid=v_west_grid,                 &
884               xmin = -0.5_wp * dx, xmax = -0.5_wp * dx,                       &
885               ymin = dy, ymax = ly - dy,                                      &
886               x0=x0, y0=y0, z0 = z0,                                          &
887               nx = 0, ny = ny-1, nz = nz, z=z)
888
889       CALL init_grid_definition('boundary', grid=v_north_grid,                &
890               xmin = 0.5_wp * dx, xmax = lx - 0.5_wp * dx,                    &
891               ymin = ly, ymax = ly,                                           &
892               x0=x0, y0=y0, z0 = z0,                                          &
893               nx = nx, ny = 0, nz = nz, z=z)
894
895       CALL init_grid_definition('boundary', grid=v_south_grid,                &
896               xmin = 0.5_wp * dx, xmax = lx - 0.5_wp * dx,                    &
897               ymin = 0.0_wp, ymax = 0.0_wp,                                   &
898               x0=x0, y0=y0, z0 = z0,                                          &
899               nx = nx, ny = 0, nz = nz, z=z)
900
901       CALL init_grid_definition('boundary', grid=v_top_grid,                  &
902               xmin = 0.5_wp * dx, xmax = lx - 0.5_wp * dx,                    &
903               ymin = dy, ymax = ly - dy,                                      &
904               x0=x0, y0=y0, z0 = z0,                                          &
905               nx = nx, ny = ny-1, nz = 1, z=(/z_top/))
906
907       CALL init_grid_definition('boundary', grid=w_east_grid,                 &
908               xmin = lx + 0.5_wp * dx, xmax = lx + 0.5_wp * dx,               &
909               ymin =  0.5_wp * dy, ymax = ly - 0.5_wp * dy,                   &
910               x0=x0, y0=y0, z0 = z0,                                          &
911               nx = 0, ny = ny, nz = nz - 1, z=zw)
912
913       CALL init_grid_definition('boundary', grid=w_west_grid,                 &
914               xmin = -0.5_wp * dx, xmax = -0.5_wp * dx,                       &
915               ymin = 0.5_wp * dy, ymax = ly - 0.5_wp * dy,                    &
916               x0=x0, y0=y0, z0 = z0,                                          &
917               nx = 0, ny = ny, nz = nz - 1, z=zw)
918
919       CALL init_grid_definition('boundary', grid=w_north_grid,                &
920               xmin = 0.5_wp * dx, xmax = lx - 0.5_wp * dx,                    &
921               ymin = ly + 0.5_wp * dy, ymax = ly + 0.5_wp * dy,               &
922               x0=x0, y0=y0, z0 = z0,                                          &
923               nx = nx, ny = 0, nz = nz - 1, z=zw)
924
925       CALL init_grid_definition('boundary', grid=w_south_grid,                &
926               xmin =  0.5_wp * dx, xmax = lx - 0.5_wp * dx,                   &
927               ymin = -0.5_wp * dy, ymax = -0.5_wp * dy,                       &
928               x0=x0, y0=y0, z0 = z0,                                          &
929               nx = nx, ny = 0, nz = nz - 1, z=zw)
930
931       CALL init_grid_definition('boundary', grid=w_top_grid,                  &
932               xmin = 0.5_wp * dx, xmax = lx - 0.5_wp * dx,                    &
933               ymin = 0.5_wp * dy, ymax = ly - 0.5_wp * dy,                    &
934               x0=x0, y0=y0, z0 = z0,                                          &
935               nx = nx, ny = ny, nz = 1, z=(/zw_top/))
936
937       CALL init_grid_definition('boundary intermediate', grid=scalars_east_intermediate,   &
938               xmin = lx + 0.5_wp * dx, xmax = lx + 0.5_wp * dx,               &
939               ymin =  0.5_wp * dy, ymax = ly - 0.5_wp * dy,                   &
940               x0=x0, y0=y0, z0 = z0,                                          &
941               nx = 0, ny = ny, nz = nlev - 2)
942
943       CALL init_grid_definition('boundary intermediate', grid=scalars_west_intermediate,   &
944               xmin = -0.5_wp * dx, xmax = -0.5_wp * dx,                       &
945               ymin =  0.5_wp * dy, ymax = ly - 0.5_wp * dy,                   &
946               x0=x0, y0=y0, z0 = z0,                                          &
947               nx = 0, ny = ny, nz = nlev - 2)
948
949       CALL init_grid_definition('boundary intermediate', grid=scalars_north_intermediate,  &
950               xmin = 0.5_wp * dx, xmax = lx - 0.5_wp * dx,                    &
951               ymin = ly + 0.5_wp * dy, ymax = ly + 0.5_wp * dy,               &
952               x0=x0, y0=y0, z0 = z0,                                          &
953               nx = nx, ny = 0, nz = nlev - 2)
954
955       CALL init_grid_definition('boundary intermediate', grid=scalars_south_intermediate,  &
956               xmin =  0.5_wp * dx, xmax = lx - 0.5_wp * dx,                   &
957               ymin = -0.5_wp * dy, ymax = -0.5_wp * dy,                       &
958               x0=x0, y0=y0, z0 = z0,                                          &
959               nx = nx, ny = 0, nz = nlev - 2)
960
961       CALL init_grid_definition('boundary intermediate', grid=scalars_top_intermediate,    &
962               xmin =  0.5_wp * dx, xmax = lx - 0.5_wp * dx,                   &
963               ymin =  0.5_wp * dy, ymax = ly - 0.5_wp * dy,                   &
964               x0=x0, y0=y0, z0 = z0,                                          &
965               nx = nx, ny = ny, nz = nlev - 2)
966
967       CALL init_grid_definition('boundary intermediate', grid=u_east_intermediate,         &
968               xmin = lx, xmax = lx,                                           &
969               ymin = 0.5_wp * dy, ymax = ly - 0.5_wp * dy,                    &
970               x0=x0, y0=y0, z0 = z0,                                          &
971               nx = 0, ny = ny, nz = nlev - 2)
972
973       CALL init_grid_definition('boundary intermediate', grid=u_west_intermediate,         &
974               xmin = 0.0_wp, xmax = 0.0_wp,                                   &
975               ymin = 0.5_wp * dy, ymax = ly - 0.5_wp * dy,                    &
976               x0=x0, y0=y0, z0 = z0,                                          &
977               nx = 0, ny = ny, nz = nlev - 2)
978
979       CALL init_grid_definition('boundary intermediate', grid=u_north_intermediate,        &
980               xmin = dx, xmax = lx - dx,                                      &
981               ymin = ly + 0.5_wp * dy, ymax = ly + 0.5_wp * dy,               &
982               x0=x0, y0=y0, z0 = z0,                                          &
983               nx = nx-1, ny = 0, nz = nlev - 2)
984
985       CALL init_grid_definition('boundary intermediate', grid=u_south_intermediate,        &
986               xmin = dx, xmax = lx - dx,                                      &
987               ymin = -0.5_wp * dy, ymax = -0.5_wp * dy,                       &
988               x0=x0, y0=y0, z0 = z0,                                          &
989               nx = nx-1, ny = 0, nz = nlev - 2)
990
991       CALL init_grid_definition('boundary intermediate', grid=u_top_intermediate,          &
992               xmin = dx, xmax = lx - dx,                                      &
993               ymin = 0.5_wp * dy, ymax = ly - 0.5_wp * dy,                    &
994               x0=x0, y0=y0, z0 = z0,                                          &
995               nx = nx-1, ny = ny, nz = nlev - 2)
996
997       CALL init_grid_definition('boundary intermediate', grid=v_east_intermediate,         &
998               xmin = lx + 0.5_wp * dx, xmax = lx + 0.5_wp * dx,               &
999               ymin = dy, ymax = ly - dy,                                      &
1000               x0=x0, y0=y0, z0 = z0,                                          &
1001               nx = 0, ny = ny-1, nz = nlev - 2)
1002
1003       CALL init_grid_definition('boundary intermediate', grid=v_west_intermediate,         &
1004               xmin = -0.5_wp * dx, xmax = -0.5_wp * dx,                       &
1005               ymin = dy, ymax = ly - dy,                                      &
1006               x0=x0, y0=y0, z0 = z0,                                          &
1007               nx = 0, ny = ny-1, nz = nlev - 2)
1008
1009       CALL init_grid_definition('boundary intermediate', grid=v_north_intermediate,        &
1010               xmin = 0.5_wp * dx, xmax = lx - 0.5_wp * dx,                    &
1011               ymin = ly, ymax = ly,                                           &
1012               x0=x0, y0=y0, z0 = z0,                                          &
1013               nx = nx, ny = 0, nz = nlev - 2)
1014
1015       CALL init_grid_definition('boundary intermediate', grid=v_south_intermediate,        &
1016               xmin = 0.5_wp * dx, xmax = lx - 0.5_wp * dx,                    &
1017               ymin = 0.0_wp, ymax = 0.0_wp,                                   &
1018               x0=x0, y0=y0, z0 = z0,                                          &
1019               nx = nx, ny = 0, nz = nlev - 2)
1020
1021       CALL init_grid_definition('boundary intermediate', grid=v_top_intermediate,          &
1022               xmin = 0.5_wp * dx, xmax = lx - 0.5_wp * dx,                    &
1023               ymin = dy, ymax = ly - dy,                                      &
1024               x0=x0, y0=y0, z0 = z0,                                          &
1025               nx = nx, ny = ny-1, nz = nlev - 2)
1026
1027       CALL init_grid_definition('boundary intermediate', grid=w_east_intermediate,         &
1028               xmin = lx + 0.5_wp * dx, xmax = lx + 0.5_wp * dx,               &
1029               ymin =  0.5_wp * dy, ymax = ly - 0.5_wp * dy,                   &
1030               x0=x0, y0=y0, z0 = z0,                                          &
1031               nx = 0, ny = ny, nz = nlev - 1)
1032
1033       CALL init_grid_definition('boundary intermediate', grid=w_west_intermediate,         &
1034               xmin = -0.5_wp * dx, xmax = -0.5_wp * dx,                       &
1035               ymin =  0.5_wp * dy, ymax = ly - 0.5_wp * dy,                   &
1036               x0=x0, y0=y0, z0 = z0,                                          &
1037               nx = 0, ny = ny, nz = nlev - 1)
1038
1039       CALL init_grid_definition('boundary intermediate', grid=w_north_intermediate,        &
1040               xmin = 0.5_wp * dx, xmax = lx - 0.5_wp * dx,                    &
1041               ymin = ly + 0.5_wp * dy, ymax = ly + 0.5_wp * dy,               &
1042               x0=x0, y0=y0, z0 = z0,                                          &
1043               nx = nx, ny = 0, nz = nlev - 1)
1044
1045       CALL init_grid_definition('boundary intermediate', grid=w_south_intermediate,        &
1046               xmin =  0.5_wp * dx, xmax = lx - 0.5_wp * dx,                   &
1047               ymin = -0.5_wp * dy, ymax = -0.5_wp * dy,                       &
1048               x0=x0, y0=y0, z0 = z0,                                          &
1049               nx = nx, ny = 0, nz = nlev - 1)
1050
1051       CALL init_grid_definition('boundary intermediate', grid=w_top_intermediate,          &
1052               xmin =  0.5_wp * dx, xmax = lx - 0.5_wp * dx,                   &
1053               ymin =  0.5_wp * dy, ymax = ly - 0.5_wp * dy,                   &
1054               x0=x0, y0=y0, z0 = z0,                                          &
1055               nx = nx, ny = ny, nz = nlev - 1)
1056    ENDIF
1057
1058!                                                                             
1059!------------------------------------------------------------------------------
1060! Section 3: Define profile grids
1061!------------------------------------------------------------------------------
1062
1063    lonmin_palm = MINVAL(palm_intermediate%clon)
1064    lonmax_palm = MAXVAL(palm_intermediate%clon)
1065    latmin_palm = MINVAL(palm_intermediate%clat)
1066    latmax_palm = MAXVAL(palm_intermediate%clat)
1067
1068    CALL init_averaging_grid(averaged_initial_scalar_profile, cosmo_grid,   &
1069            x = 0.5_wp * lx, y = 0.5_wp * ly, z = z, z0 = z0,               &
1070            lonmin = lonmin_palm, lonmax = lonmax_palm,                     &
1071            latmin = latmin_palm, latmax = latmax_palm,                     &
1072            kind='scalar', name='averaged initial scalar')
1073
1074    CALL init_averaging_grid(averaged_initial_w_profile, cosmo_grid,        &
1075            x = 0.5_wp * lx, y = 0.5_wp * ly, z = zw, z0 = z0,              &
1076            lonmin = lonmin_palm, lonmax = lonmax_palm,                     &
1077            latmin = latmin_palm, latmax = latmax_palm,                     &
1078            kind='w', name='averaged initial w')
1079
1080    CALL init_averaging_grid(averaged_scalar_profile, cosmo_grid,           &
1081            x = 0.5_wp * lx, y = 0.5_wp * ly, z = z, z0 = z0,               &
1082            lonmin = lam_west, lonmax = lam_east,                           &
1083            latmin = phi_south, latmax = phi_north,                         &
1084            kind='scalar', name='centre geostrophic scalar')
1085
1086    CALL init_averaging_grid(averaged_w_profile, cosmo_grid,                &
1087            x = 0.5_wp * lx, y = 0.5_wp * ly, z = zw, z0 = z0,              &
1088            lonmin = lam_west, lonmax = lam_east,                           &
1089            latmin = phi_south, latmax = phi_north,                         &
1090            kind='w', name='centre geostrophic w')
1091
1092    CALL init_averaging_grid(south_averaged_scalar_profile, cosmo_grid,     &
1093            x = 0.5_wp * lx, y = 0.5_wp * ly, z = z, z0 = z0,               &
1094            lonmin = lam_west, lonmax = lam_east,                           &
1095            latmin = phi_centre - averaging_angle,                          &
1096            latmax = phi_centre,                                            &
1097            kind='scalar', name='south geostrophic scalar')
1098
1099    CALL init_averaging_grid(north_averaged_scalar_profile, cosmo_grid,     &
1100            x = 0.5_wp * lx, y = 0.5_wp * ly, z = z, z0 = z0,               &
1101            lonmin = lam_west, lonmax = lam_east,                           &
1102            latmin = phi_centre,                                            &
1103            latmax = phi_centre + averaging_angle,                          &
1104            kind='scalar', name='north geostrophic scalar')
1105
1106    CALL init_averaging_grid(west_averaged_scalar_profile, cosmo_grid,      &
1107            x = 0.5_wp * lx, y = 0.5_wp * ly, z = z, z0 = z0,               &
1108            lonmin = lam_centre - averaging_angle,                          &
1109            lonmax = lam_centre,                                            &
1110            latmin = phi_south, latmax = phi_north,                         &
1111            kind='scalar', name='west geostrophic scalar')
1112
1113    CALL init_averaging_grid(east_averaged_scalar_profile, cosmo_grid,      &
1114            x = 0.5_wp * lx, y = 0.5_wp * ly, z = z, z0 = z0,               &
1115            lonmin = lam_centre,                                            &
1116            lonmax = lam_centre + averaging_angle,                          &
1117            latmin = phi_south, latmax = phi_north,                         &
1118            kind='scalar', name='east geostrophic scalar')
1119
1120
1121!                                                                             
1122!------------------------------------------------------------------------------
1123! Section 4: Precompute neighbours and weights for interpolation             
1124!------------------------------------------------------------------------------
1125    interp_mode = 's'
1126    CALL setup_interpolation(cosmo_grid, palm_grid, palm_intermediate, interp_mode, ic_mode=cfg%ic_mode)
1127    IF (boundary_variables_required)  THEN
1128       CALL setup_interpolation(cosmo_grid, scalars_east_grid, scalars_east_intermediate, interp_mode)
1129       CALL setup_interpolation(cosmo_grid, scalars_west_grid, scalars_west_intermediate, interp_mode)
1130       CALL setup_interpolation(cosmo_grid, scalars_north_grid, scalars_north_intermediate, interp_mode)
1131       CALL setup_interpolation(cosmo_grid, scalars_south_grid, scalars_south_intermediate, interp_mode)
1132       CALL setup_interpolation(cosmo_grid, scalars_top_grid, scalars_top_intermediate, interp_mode)
1133    ENDIF
1134
1135    interp_mode = 'u'
1136    CALL setup_interpolation(cosmo_grid, u_initial_grid, u_initial_intermediate, interp_mode, ic_mode=cfg%ic_mode)
1137    IF (boundary_variables_required)  THEN
1138       CALL setup_interpolation(cosmo_grid, u_east_grid, u_east_intermediate, interp_mode)
1139       CALL setup_interpolation(cosmo_grid, u_west_grid, u_west_intermediate, interp_mode)
1140       CALL setup_interpolation(cosmo_grid, u_north_grid, u_north_intermediate, interp_mode)
1141       CALL setup_interpolation(cosmo_grid, u_south_grid, u_south_intermediate, interp_mode)
1142       CALL setup_interpolation(cosmo_grid, u_top_grid, u_top_intermediate, interp_mode)
1143    ENDIF
1144
1145    interp_mode = 'v'
1146    CALL setup_interpolation(cosmo_grid, v_initial_grid, v_initial_intermediate, interp_mode, ic_mode=cfg%ic_mode)
1147    IF (boundary_variables_required)  THEN
1148       CALL setup_interpolation(cosmo_grid, v_east_grid, v_east_intermediate, interp_mode)
1149       CALL setup_interpolation(cosmo_grid, v_west_grid, v_west_intermediate, interp_mode)
1150       CALL setup_interpolation(cosmo_grid, v_north_grid, v_north_intermediate, interp_mode)
1151       CALL setup_interpolation(cosmo_grid, v_south_grid, v_south_intermediate, interp_mode)
1152       CALL setup_interpolation(cosmo_grid, v_top_grid, v_top_intermediate, interp_mode)
1153    ENDIF
1154
1155    interp_mode = 'w'
1156    CALL setup_interpolation(cosmo_grid, w_initial_grid, w_initial_intermediate, interp_mode, ic_mode=cfg%ic_mode)
1157    IF (boundary_variables_required)  THEN
1158       CALL setup_interpolation(cosmo_grid, w_east_grid, w_east_intermediate, interp_mode)
1159       CALL setup_interpolation(cosmo_grid, w_west_grid, w_west_intermediate, interp_mode)
1160       CALL setup_interpolation(cosmo_grid, w_north_grid, w_north_intermediate, interp_mode)
1161       CALL setup_interpolation(cosmo_grid, w_south_grid, w_south_intermediate, interp_mode)
1162       CALL setup_interpolation(cosmo_grid, w_top_grid, w_top_intermediate, interp_mode)
1163    ENDIF
1164
1165    IF (TRIM(cfg%ic_mode) == 'profile')  THEN
1166        !TODO: remove this conditional if not needed.
1167    ENDIF
1168
1169 END SUBROUTINE setup_grids
1170
1171
1172!------------------------------------------------------------------------------!
1173! Description:
1174! ------------
1175!> Driver for computing neighbour indices and weights for horizontal and
1176!> vertical interpolation.
1177!------------------------------------------------------------------------------!
1178 SUBROUTINE setup_interpolation(cosmo_grid, grid, intermediate_grid, kind, ic_mode)
1179
1180    TYPE(grid_definition), INTENT(IN), TARGET    ::  cosmo_grid
1181    TYPE(grid_definition), INTENT(INOUT), TARGET ::  grid, intermediate_grid
1182    CHARACTER, INTENT(IN)                        ::  kind
1183    CHARACTER(LEN=*), INTENT(IN), OPTIONAL       ::  ic_mode
1184
1185    REAL(wp), DIMENSION(:), POINTER     ::  cosmo_lat, cosmo_lon
1186    REAL(wp), DIMENSION(:,:,:), POINTER ::  cosmo_h
1187
1188    LOGICAL :: setup_volumetric
1189
1190!------------------------------------------------------------------------------
1191! Section 1: Horizontal interpolation                                       
1192!------------------------------------------------------------------------------
1193!
1194!-- Select horizontal coordinates according to kind of points (s/w, u, v)
1195    SELECT CASE(kind)
1196
1197!
1198!-- scalars
1199    CASE('s')
1200
1201       cosmo_lat => cosmo_grid%lat
1202       cosmo_lon => cosmo_grid%lon
1203       cosmo_h   => cosmo_grid%hfl
1204!
1205!-- vertical velocity
1206    CASE('w')
1207
1208       cosmo_lat => cosmo_grid%lat
1209       cosmo_lon => cosmo_grid%lon
1210       cosmo_h   => cosmo_grid%hhl
1211!
1212!-- x velocity
1213    CASE('u')
1214
1215       cosmo_lat => cosmo_grid%lat
1216       cosmo_lon => cosmo_grid%lonu
1217       cosmo_h   => cosmo_grid%hfl
1218
1219!
1220!-- y velocity
1221    CASE('v')
1222
1223       cosmo_lat => cosmo_grid%latv
1224       cosmo_lon => cosmo_grid%lon
1225       cosmo_h   => cosmo_grid%hfl
1226
1227    CASE DEFAULT
1228
1229       message = "Interpolation quantity '" // kind // "' is not supported."
1230       CALL inifor_abort('setup_interpolation', message)
1231
1232    END SELECT
1233
1234    CALL find_horizontal_neighbours(cosmo_lat, cosmo_lon,                   &
1235       intermediate_grid%clat, intermediate_grid%clon,                  &
1236       intermediate_grid%ii, intermediate_grid%jj)
1237
1238    CALL compute_horizontal_interp_weights(cosmo_lat, cosmo_lon,            &
1239       intermediate_grid%clat, intermediate_grid%clon,                  &
1240       intermediate_grid%ii, intermediate_grid%jj,                      &
1241       intermediate_grid%w_horiz)
1242
1243!------------------------------------------------------------------------------
1244! Section 2: Vertical interpolation
1245!------------------------------------------------------------------------------
1246
1247!
1248!-- If profile initialization is chosen, we--somewhat counterintuitively--
1249!-- don't need to compute vertical interpolation weights. At least, we
1250!-- don't need them on the intermediate grid, which fills the entire PALM
1251!-- domain volume. Instead we need vertical weights for the intermediate
1252!-- profile grids, which get computed in setup_averaging().
1253    setup_volumetric = .TRUE.
1254    IF (PRESENT(ic_mode))  THEN
1255       IF (TRIM(ic_mode) == 'profile')  setup_volumetric = .FALSE.
1256    ENDIF
1257
1258    IF (setup_volumetric)  THEN
1259       ALLOCATE( intermediate_grid%h(0:intermediate_grid%nx,            &
1260                                     0:intermediate_grid%ny,            &
1261                                     0:intermediate_grid%nz) ) 
1262       intermediate_grid%h(:,:,:) = - EARTH_RADIUS
1263
1264!
1265!--    For w points, use hhl, for scalars use hfl
1266!--    compute the full heights for the intermediate grids
1267       CALL interpolate_2d(cosmo_h, intermediate_grid%h, intermediate_grid)
1268       CALL find_vertical_neighbours_and_weights_interp(grid, intermediate_grid)
1269    ENDIF
1270       
1271 END SUBROUTINE setup_interpolation
1272
1273!------------------------------------------------------------------------------!
1274! Description:
1275! ------------
1276!> Initializes grid_definition-type variables.
1277!>
1278!> Input parameters:
1279!> -----------------
1280!> kind : Grid kind, distinguishes between PALM-4U and COSMO-DE grids
1281!>    as well as grids covering the boundary surfaces. Valid kinds are:
1282!>       - 'palm'
1283!>       - 'cosmo-de'
1284!>       - 'eastwest-scalar'
1285!>
1286!> <xyx>min, <xyz>max : Domain minima and maxima in x, y, and z direction. Note
1287!>    that these values do not necessarily translate to the outmost coordinates
1288!>    of the generated grid but rather refer to the extent of the underlying
1289!>    PALM-4U computational domain (i.e. the outer cell faces). The coordinates
1290!>    of the generated grid will be inferred from this information taking into
1291!>    account the initialization mode ic_mode. For example, the coordinates of a
1292!>    boundary grid initialized using mode 'eastwest-scalar' will be located in
1293!>    planes one half grid point outwards of xmin and xmax.
1294!>
1295!> z0 : Elevation of the PALM-4U domain above sea level [m]
1296!>
1297!> n<xyz> : Number of grod points in x, y, and z direction
1298!>
1299!> Output parameters:
1300!> ------------------
1301!> grid : Grid variable to be initialized.
1302!------------------------------------------------------------------------------!
1303 SUBROUTINE init_grid_definition(kind, xmin, xmax, ymin, ymax,              &
1304                                 x0, y0, z0, nx, ny, nz, z, zw, grid, ic_mode)
1305    CHARACTER(LEN=*), INTENT(IN)           ::  kind
1306    CHARACTER(LEN=*), INTENT(IN), OPTIONAL ::  ic_mode
1307    INTEGER, INTENT(IN)                    ::  nx, ny, nz
1308    REAL(wp), INTENT(IN)                   ::  xmin, xmax, ymin, ymax
1309    REAL(wp), INTENT(IN)                   ::  x0, y0, z0
1310    REAL(wp), INTENT(IN), TARGET, OPTIONAL ::  z(:)
1311    REAL(wp), INTENT(IN), TARGET, OPTIONAL ::  zw(:)
1312    TYPE(grid_definition), INTENT(INOUT)   ::  grid
1313
1314    grid%nx = nx
1315    grid%ny = ny
1316    grid%nz = nz
1317
1318    grid%lx = xmax - xmin
1319    grid%ly = ymax - ymin
1320
1321    grid%x0 = x0
1322    grid%y0 = y0
1323    grid%z0 = z0
1324
1325    SELECT CASE( TRIM(kind) )
1326
1327       CASE('boundary')
1328       
1329          IF (.NOT.PRESENT(z))  THEN
1330             message = "z has not been passed but is required for 'boundary' grids"
1331             CALL inifor_abort('init_grid_definition', message)
1332          ENDIF
1333       
1334          ALLOCATE( grid%x(0:nx) )
1335          CALL linspace(xmin, xmax, grid%x)
1336       
1337          ALLOCATE( grid%y(0:ny) )
1338          CALL linspace(ymin, ymax, grid%y)
1339       
1340          grid%z => z
1341       
1342!     
1343!--       Allocate neighbour indices and weights
1344          IF (TRIM(ic_mode) .NE. 'profile')  THEN
1345             ALLOCATE( grid%kk(0:nx, 0:ny, 1:nz, 2) )
1346             grid%kk(:,:,:,:) = -1
1347       
1348             ALLOCATE( grid%w_verti(0:nx, 0:ny, 1:nz, 2) )
1349             grid%w_verti(:,:,:,:) = 0.0_wp
1350          ENDIF
1351       
1352       CASE('boundary intermediate')
1353       
1354          ALLOCATE( grid%x(0:nx) )
1355          CALL linspace(xmin, xmax, grid%x)
1356       
1357          ALLOCATE( grid%y(0:ny) )
1358          CALL linspace(ymin, ymax, grid%y)
1359       
1360          ALLOCATE( grid%clon(0:nx, 0:ny), grid%clat(0:nx, 0:ny)  )
1361       
1362          CALL rotate_to_cosmo(                                               &
1363             phir = project( grid%y, y0, EARTH_RADIUS ) ,                   & ! = plate-carree latitude
1364             lamr = project( grid%x, x0, EARTH_RADIUS ) ,                   & ! = plate-carree longitude
1365             phip = phi_cn, lamp = lambda_cn,                                 &
1366             phi  = grid%clat,                                              &
1367             lam  = grid%clon,                                              &
1368             gam  = gam                                                       &
1369          )
1370       
1371!     
1372!--       Allocate neighbour indices and weights
1373          ALLOCATE( grid%ii(0:nx, 0:ny, 4),                                 &
1374                    grid%jj(0:nx, 0:ny, 4) )
1375          grid%ii(:,:,:)   = -1
1376          grid%jj(:,:,:)   = -1
1377       
1378          ALLOCATE( grid%w_horiz(0:nx, 0:ny, 4) )
1379          grid%w_horiz(:,:,:)   = 0.0_wp
1380       
1381!     
1382!--    This mode initializes a Cartesian PALM-4U grid and adds the
1383!--    corresponding latitudes and longitudes of the rotated pole grid.
1384       CASE('palm')
1385       
1386          IF (.NOT.PRESENT(z))  THEN
1387             message = "z has not been passed but is required for 'palm' grids"
1388             CALL inifor_abort('init_grid_definition', message)
1389          ENDIF
1390       
1391          IF (.NOT.PRESENT(zw))  THEN
1392             message = "zw has not been passed but is required for 'palm' grids"
1393             CALL inifor_abort('init_grid_definition', message)
1394          ENDIF
1395       
1396          grid%name(1) = 'x and lon'
1397          grid%name(2) = 'y and lat'
1398          grid%name(3) = 'z'
1399       
1400!     
1401!--       TODO: Remove use of global dx, dy, dz variables. Consider
1402!--       TODO: associating global x,y, and z arrays.
1403          ALLOCATE( grid%x(0:nx),   grid%y(0:ny) )
1404          ALLOCATE( grid%xu(1:nx),  grid%yv(1:ny) )
1405          CALL linspace(xmin + 0.5_wp* dx, xmax - 0.5_wp* dx, grid%x)
1406          CALL linspace(ymin + 0.5_wp* dy, ymax - 0.5_wp* dy, grid%y)
1407          grid%z => z
1408          CALL linspace(xmin +  dx, xmax -  dx, grid%xu)
1409          CALL linspace(ymin +  dy, ymax -  dy, grid%yv)
1410          grid%zw => zw
1411       
1412          grid%depths => depths
1413       
1414!     
1415!--       Allocate neighbour indices and weights
1416          IF (TRIM(ic_mode) .NE. 'profile')  THEN
1417             ALLOCATE( grid%kk(0:nx, 0:ny, 1:nz, 2) )
1418             grid%kk(:,:,:,:) = -1
1419       
1420             ALLOCATE( grid%w_verti(0:nx, 0:ny, 1:nz, 2) )
1421             grid%w_verti(:,:,:,:) = 0.0_wp
1422          ENDIF
1423       
1424       CASE('palm intermediate')
1425       
1426          grid%name(1) = 'x and lon'
1427          grid%name(2) = 'y and lat'
1428          grid%name(3) = 'interpolated hhl or hfl'
1429       
1430!     
1431!--       TODO: Remove use of global dx, dy, dz variables. Consider
1432!--       TODO: associating global x,y, and z arrays.
1433          ALLOCATE( grid%x(0:nx),   grid%y(0:ny) )
1434          ALLOCATE( grid%xu(1:nx),  grid%yv(1:ny) )
1435          CALL linspace(xmin + 0.5_wp*dx, xmax - 0.5_wp*dx, grid%x)
1436          CALL linspace(ymin + 0.5_wp*dy, ymax - 0.5_wp*dy, grid%y)
1437          CALL linspace(xmin + dx, xmax - dx, grid%xu)
1438          CALL linspace(ymin + dy, ymax - dy, grid%yv)
1439       
1440          grid%depths => depths
1441       
1442!     
1443!--       Allocate rotated-pole coordinates, clon is for (c)osmo-de (lon)gitude
1444          ALLOCATE( grid%clon(0:nx, 0:ny),   grid%clat(0:nx, 0:ny)  )
1445          ALLOCATE( grid%clonu(1:nx, 0:ny),  grid%clatu(1:nx, 0:ny) )
1446          ALLOCATE( grid%clonv(0:nx, 1:ny),  grid%clatv(0:nx, 1:ny) )
1447       
1448!     
1449!--       Compute rotated-pole coordinates of...
1450!--       ... PALM-4U centres
1451          CALL rotate_to_cosmo(                                               &
1452             phir = project( grid%y, y0, EARTH_RADIUS ) , & ! = plate-carree latitude
1453             lamr = project( grid%x, x0, EARTH_RADIUS ) , & ! = plate-carree longitude
1454             phip = phi_cn, lamp = lambda_cn,                                 &
1455             phi  = grid%clat,                                              &
1456             lam  = grid%clon,                                              &
1457             gam  = gam                                                       &
1458          )
1459       
1460!     
1461!--       ... PALM-4U u winds
1462          CALL rotate_to_cosmo(                                               &
1463             phir = project( grid%y,  y0, EARTH_RADIUS ), & ! = plate-carree latitude
1464             lamr = project( grid%xu, x0, EARTH_RADIUS ), & ! = plate-carree longitude
1465             phip = phi_cn, lamp = lambda_cn,                                 &
1466             phi  = grid%clatu,                                             &
1467             lam  = grid%clonu,                                             &
1468             gam  = gam                                                       &
1469          )
1470       
1471!     
1472!--       ... PALM-4U v winds
1473          CALL rotate_to_cosmo(                                               &
1474             phir = project( grid%yv, y0, EARTH_RADIUS ), & ! = plate-carree latitude
1475             lamr = project( grid%x,  x0, EARTH_RADIUS ), & ! = plate-carree longitude
1476             phip = phi_cn, lamp = lambda_cn,                                 &
1477             phi  = grid%clatv,                                             &
1478             lam  = grid%clonv,                                             &
1479             gam  = gam                                                       &
1480          )
1481       
1482!     
1483!--       Allocate neighbour indices and weights
1484          ALLOCATE( grid%ii(0:nx, 0:ny, 4),                                 &
1485                    grid%jj(0:nx, 0:ny, 4) )
1486          grid%ii(:,:,:)   = -1
1487          grid%jj(:,:,:)   = -1
1488       
1489          ALLOCATE( grid%w_horiz(0:nx, 0:ny, 4) )
1490          grid%w_horiz(:,:,:)   = 0.0_wp
1491       
1492       CASE('cosmo-de')
1493          grid%name(1) = 'rlon'         ! of COMSO-DE cell centres (scalars)
1494          grid%name(2) = 'rlat'         ! of COMSO-DE cell centres (scalars)
1495          grid%name(3) = 'height'
1496       
1497          ALLOCATE( grid%lon(0:nx),   grid%lat(0:ny)  )
1498          ALLOCATE( grid%lonu(0:nx),  grid%latv(0:ny) )
1499       
1500          CALL linspace(xmin, xmax, grid%lon)
1501          CALL linspace(ymin, ymax, grid%lat)
1502          grid%lonu(:) = grid%lon + 0.5_wp * (grid%lx / grid%nx)
1503          grid%latv(:) = grid%lat + 0.5_wp * (grid%ly / grid%ny)
1504       
1505!     
1506!--       Point to heights of half levels (hhl) and compute heights of full
1507!--       levels (hfl) as arithmetic averages
1508          grid%hhl => hhl
1509          grid%hfl => hfl
1510          grid%depths => depths
1511       
1512       CASE DEFAULT
1513           message = "Grid kind '" // TRIM(kind) // "' is not recognized."
1514           CALL inifor_abort('init_grid_definition', message)
1515
1516    END SELECT
1517
1518 END SUBROUTINE init_grid_definition
1519
1520
1521!------------------------------------------------------------------------------!
1522! Description:
1523! ------------
1524!> Initializes averaging grids
1525!>
1526!> Input parameters:
1527!> -----------------
1528!>
1529!> cosmo_grid : grid_definition-type varieble initialized as 'cosmo-de' grid
1530!>    providing longitudes and latitudes of the parent grid
1531!>
1532!> x, y : location of the profile to be averaged in the PALM system [m]
1533!>
1534!> z : vertical grid coordinates of the profile in the PALM system [m]
1535!>
1536!> <lat/lon><min/max>: boundaries of the averaging region in the parent system,
1537!>    i.e. the COSMO-DE rotated-pole system. [rad]
1538!>
1539!> kind : kind of quantity to be averaged using this averaging grid.
1540!>    Destinguishes COSMO-DE scalar and w-velocity levels. Note that finding the
1541!>    parent/COSMO columns for the region in get_cosmo_averaging_region() is
1542!>    independent of 'kind' b/ecause INIFOR uses column-centred u and v velocity
1543!>    components, which are computed in the preprocessing step.
1544!>
1545!> Output parameters:
1546!> ------------------
1547!> avg_grid : averagin grid to be initialized
1548!------------------------------------------------------------------------------!
1549 SUBROUTINE init_averaging_grid(avg_grid, cosmo_grid, x, y, z, z0,          &
1550    lonmin, lonmax, latmin, latmax, kind, name)
1551
1552    TYPE(grid_definition), INTENT(INOUT) ::  avg_grid
1553    TYPE(grid_definition), INTENT(IN)    ::  cosmo_grid
1554    REAL(wp), INTENT(IN)                 ::  x, y, z0
1555    REAL(wp), INTENT(IN), TARGET         ::  z(:)
1556    REAL(wp), INTENT(IN)                 ::  lonmin !< lower longitude bound of the averaging grid region [COSMO rotated-pole rad]
1557    REAL(wp), INTENT(IN)                 ::  lonmax !< upper longitude bound of the averaging grid region [COSMO rotated-pole rad]
1558    REAL(wp), INTENT(IN)                 ::  latmin !< lower latitude bound of the averaging grid region [COSMO rotated-pole rad]
1559    REAL(wp), INTENT(IN)                 ::  latmax !< lower latitude bound of the averaging grid region [COSMO rotated-pole rad]
1560
1561    CHARACTER(LEN=*), INTENT(IN)         ::  kind
1562    CHARACTER(LEN=*), INTENT(IN)         ::  name
1563
1564    LOGICAL                              ::  level_based_averaging
1565
1566    ALLOCATE( avg_grid%x(1) )
1567    ALLOCATE( avg_grid%y(1) )
1568    avg_grid%x(1) = x
1569    avg_grid%y(1) = y
1570    avg_grid%z => z
1571    avg_grid%z0 = z0
1572
1573    avg_grid%nz = SIZE(z, 1)
1574
1575    ALLOCATE( avg_grid%lon(2) )
1576    ALLOCATE( avg_grid%lat(2) )
1577    avg_grid%lon(1:2) = (/lonmin, lonmax/)
1578    avg_grid%lat(1:2) = (/latmin, latmax/)
1579
1580    avg_grid%kind = TRIM(kind)
1581    avg_grid%name(1) = TRIM(name)
1582
1583!
1584!-- Find and store COSMO columns that fall into the coordinate range
1585!-- given by avg_grid%clon, %clat
1586    CALL get_cosmo_averaging_region(avg_grid, cosmo_grid)
1587
1588    ALLOCATE (avg_grid%kkk(avg_grid%n_columns, avg_grid%nz, 2) )
1589    ALLOCATE (avg_grid%w(avg_grid%n_columns, avg_grid%nz, 2) )
1590!
1591!-- Compute average COSMO levels in the averaging region
1592    SELECT CASE(avg_grid%kind)
1593
1594       CASE('scalar', 'u', 'v')
1595          avg_grid%cosmo_h => cosmo_grid%hfl
1596
1597       CASE('w')
1598          avg_grid%cosmo_h => cosmo_grid%hhl
1599
1600       CASE DEFAULT
1601          message = "Averaging grid kind '" // TRIM(avg_grid%kind) // &
1602                    "' is not supported. Use 'scalar', 'u', or 'v'."
1603          CALL inifor_abort('get_cosmo_averaging_region', message)
1604
1605    END SELECT
1606
1607!
1608!-- For level-besed averaging, compute average heights
1609    level_based_averaging = ( TRIM(cfg%averaging_mode) == 'level' )
1610    IF (level_based_averaging)  THEN
1611       ALLOCATE(avg_grid%h(1,1,SIZE(avg_grid%cosmo_h, 3)) )
1612 
1613       CALL average_2d(avg_grid%cosmo_h, avg_grid%h(1,1,:),             &
1614                       avg_grid%iii, avg_grid%jjj)
1615
1616    ENDIF
1617
1618!
1619!-- Compute vertical weights and neighbours
1620    CALL find_vertical_neighbours_and_weights_average(                      &
1621       avg_grid, level_based_averaging                                      &
1622    )
1623
1624 END SUBROUTINE init_averaging_grid
1625
1626
1627 SUBROUTINE get_cosmo_averaging_region(avg_grid, cosmo_grid)
1628    TYPE(grid_definition), INTENT(INOUT)         ::  avg_grid
1629    TYPE(grid_definition), TARGET, INTENT(IN)    ::  cosmo_grid
1630
1631    REAL(wp), DIMENSION(:), POINTER              ::  cosmo_lon, cosmo_lat
1632    REAL(wp)                                     ::  dlon, dlat
1633
1634    INTEGER ::  i, j, imin, imax, jmin, jmax, l, nx, ny
1635
1636
1637    SELECT CASE( TRIM(avg_grid%kind) )
1638
1639       CASE('scalar', 'w')
1640          cosmo_lon => cosmo_grid%lon
1641          cosmo_lat => cosmo_grid%lat
1642
1643       CASE('u')
1644          cosmo_lon => cosmo_grid%lonu
1645          cosmo_lat => cosmo_grid%lat
1646
1647       CASE('v')
1648          cosmo_lon => cosmo_grid%lon
1649          cosmo_lat => cosmo_grid%latv
1650
1651       CASE DEFAULT
1652          message = "Averaging grid kind '" // TRIM(avg_grid%kind) // &
1653                    "' is not supported. Use 'scalar', 'u', or 'v'."
1654          CALL inifor_abort('get_cosmo_averaging_region', message)
1655
1656    END SELECT
1657
1658!
1659!-- FIXME: make dlon, dlat parameters of the grid_defintion type
1660    dlon = cosmo_lon(1) - cosmo_lon(0)
1661    dlat = cosmo_lat(1) - cosmo_lat(0)
1662
1663    imin = FLOOR  ( (avg_grid%lon(1) - cosmo_lon(0)) / dlon )
1664    imax = CEILING( (avg_grid%lon(2) - cosmo_lon(0)) / dlon )
1665
1666    jmin = FLOOR  ( (avg_grid%lat(1) - cosmo_lat(0)) / dlat )
1667    jmax = CEILING( (avg_grid%lat(2) - cosmo_lat(0)) / dlat )
1668   
1669    message = "Grid " // TRIM(avg_grid%name(1)) // " averages over " // &
1670              TRIM(str(imin)) // " <= i <= " // TRIM(str(imax)) //          &
1671              " and " //                                                    &
1672              TRIM(str(jmin)) // " <= j <= " // TRIM(str(jmax))
1673    CALL report( 'get_cosmo_averaging_region', message )
1674
1675    nx = imax - imin + 1
1676    ny = jmax - jmin + 1
1677    avg_grid%n_columns = nx * ny
1678
1679    ALLOCATE( avg_grid%iii(avg_grid%n_columns),                         &
1680              avg_grid%jjj(avg_grid%n_columns) )
1681
1682    l = 0
1683    DO  j = jmin, jmax
1684    DO  i = imin, imax
1685       l = l + 1
1686       avg_grid%iii(l) = i
1687       avg_grid%jjj(l) = j
1688    ENDDO
1689    ENDDO
1690
1691 END SUBROUTINE get_cosmo_averaging_region
1692
1693
1694!------------------------------------------------------------------------------!
1695! Description:
1696! ------------
1697!> PALM's stretched vertical grid generator. Forked from PALM revision 3139, see
1698!> https://palm.muk.uni-hannover.de/trac/browser/palm/trunk/SOURCE/init_grid.f90?rev=3139
1699!>
1700!> This routine computes the levels of scalar points. The levels of the velocity
1701!> points are then obtained as the midpoints inbetween using the INIFOR routine
1702!> 'modpoints'.
1703!------------------------------------------------------------------------------!
1704 SUBROUTINE stretched_z(z, dz, dz_max, dz_stretch_factor, dz_stretch_level, &
1705                        dz_stretch_level_start, dz_stretch_level_end,       &
1706                        dz_stretch_factor_array)
1707
1708    REAL(wp), DIMENSION(:), INTENT(INOUT) ::  z, dz, dz_stretch_factor_array
1709    REAL(wp), DIMENSION(:), INTENT(INOUT) ::  dz_stretch_level_start, dz_stretch_level_end
1710    REAL(wp), INTENT(IN) ::  dz_max, dz_stretch_factor, dz_stretch_level
1711
1712    INTEGER ::  number_stretch_level_start        !< number of user-specified start levels for stretching
1713    INTEGER ::  number_stretch_level_end          !< number of user-specified end levels for stretching
1714
1715    REAL(wp), DIMENSION(:), ALLOCATABLE ::  min_dz_stretch_level_end
1716    REAL(wp) ::  dz_level_end, dz_stretched
1717
1718    INTEGER ::  dz_stretch_level_end_index(9)      !< vertical grid level index until which the vertical grid spacing is stretched
1719    INTEGER ::  dz_stretch_level_start_index(9)    !< vertical grid level index above which the vertical grid spacing is stretched
1720    INTEGER ::  dz_stretch_level_index = 0
1721    INTEGER ::  k, n, number_dz
1722
1723!
1724!-- Compute height of u-levels from constant grid length and dz stretch factors
1725    IF ( dz(1) == -1.0_wp )  THEN
1726       message = 'missing dz'
1727       CALL inifor_abort( 'stretched_z', message)
1728    ELSEIF ( dz(1) <= 0.0_wp )  THEN
1729       WRITE( message, * ) 'dz=', dz(1),' <= 0.0'
1730       CALL inifor_abort( 'stretched_z', message)
1731    ENDIF
1732
1733!
1734!-- Initialize dz_stretch_level_start with the value of dz_stretch_level
1735!-- if it was set by the user
1736    IF ( dz_stretch_level /= -9999999.9_wp )  THEN
1737       dz_stretch_level_start(1) = dz_stretch_level
1738    ENDIF
1739       
1740!
1741!-- Determine number of dz values and stretching levels specified by the
1742!-- user to allow right controlling of the stretching mechanism and to
1743!-- perform error checks. The additional requirement that dz /= dz_max
1744!-- for counting number of user-specified dz values is necessary. Otherwise
1745!-- restarts would abort if the old stretching mechanism with dz_stretch_level
1746!-- is used (Attention: The user is not allowed to specify a dz value equal
1747!-- to the default of dz_max = 999.0).
1748    number_dz = COUNT( dz /= -1.0_wp .AND. dz /= dz_max )
1749    number_stretch_level_start = COUNT( dz_stretch_level_start /=           &
1750                                        -9999999.9_wp )
1751    number_stretch_level_end = COUNT( dz_stretch_level_end /=               &
1752                                      9999999.9_wp )
1753
1754!
1755!-- The number of specified end levels +1 has to be the same than the number
1756!-- of specified dz values
1757    IF ( number_dz /= number_stretch_level_end + 1 )  THEN
1758       WRITE( message, * ) 'The number of values for dz = ',                &
1759                           number_dz, 'has to be the same than ',           &
1760                           'the number of values for ',                     &
1761                           'dz_stretch_level_end + 1 = ',                   &
1762                           number_stretch_level_end+1
1763       CALL inifor_abort( 'stretched_z', message)
1764    ENDIF
1765   
1766!
1767!-- The number of specified start levels has to be the same or one less than
1768!-- the number of specified dz values
1769    IF ( number_dz /= number_stretch_level_start + 1 .AND.                  &
1770         number_dz /= number_stretch_level_start )  THEN
1771       WRITE( message, * ) 'The number of values for dz = ',         &
1772                           number_dz, 'has to be the same or one ', &
1773                           'more than& the number of values for ',  &
1774                           'dz_stretch_level_start = ',             &
1775                           number_stretch_level_start
1776       CALL inifor_abort( 'stretched_z', message)
1777    ENDIF
1778   
1779!-- The number of specified start levels has to be the same or one more than
1780!-- the number of specified end levels
1781    IF ( number_stretch_level_start /= number_stretch_level_end + 1 .AND.   &
1782         number_stretch_level_start /= number_stretch_level_end )  THEN
1783       WRITE( message, * ) 'The number of values for ',              &
1784                           'dz_stretch_level_start = ',              &
1785                           dz_stretch_level_start, 'has to be the ',&
1786                           'same or one more than& the number of ', &
1787                           'values for dz_stretch_level_end = ',    &
1788                           number_stretch_level_end
1789       CALL inifor_abort( 'stretched_z', message)
1790    ENDIF
1791
1792!
1793!-- Initialize dz for the free atmosphere with the value of dz_max
1794    IF ( dz(number_stretch_level_start+1) == -1.0_wp .AND.                     &
1795         number_stretch_level_start /= 0 )  THEN
1796       dz(number_stretch_level_start+1) = dz_max
1797    ENDIF
1798       
1799!
1800!-- Initialize the stretching factor if (infinitely) stretching in the free
1801!-- atmosphere is desired (dz_stretch_level_end was not specified for the
1802!-- free atmosphere)
1803    IF ( number_stretch_level_start == number_stretch_level_end + 1 )  THEN
1804       dz_stretch_factor_array(number_stretch_level_start) =                   &
1805       dz_stretch_factor
1806    ENDIF
1807
1808!-- Allocation of arrays for stretching
1809    ALLOCATE( min_dz_stretch_level_end(number_stretch_level_start) )
1810
1811!
1812!-- The stretching region has to be large enough to allow for a smooth
1813!-- transition between two different grid spacings
1814    DO  n = 1, number_stretch_level_start
1815       min_dz_stretch_level_end(n) = dz_stretch_level_start(n) +            &
1816                                     4 * MAX( dz(n),dz(n+1) )
1817    ENDDO
1818
1819    IF ( ANY( min_dz_stretch_level_end(1:number_stretch_level_start) >      &
1820              dz_stretch_level_end(1:number_stretch_level_start) ) )  THEN
1821    !IF ( ANY( min_dz_stretch_level_end >      &
1822    !          dz_stretch_level_end ) ) THEN
1823          message = 'Each dz_stretch_level_end has to be larger '  // &
1824                    'than its corresponding value for ' //            &
1825                    'dz_stretch_level_start + 4*MAX(dz(n),dz(n+1)) '//&
1826                    'to allow for smooth grid stretching'
1827          CALL inifor_abort('stretched_z', message)
1828    ENDIF
1829   
1830!
1831!-- Stretching must not be applied within the prandtl_layer
1832!-- (first two grid points). For the default case dz_stretch_level_start
1833!-- is negative. Therefore the absolut value is checked here.
1834    IF ( ANY( ABS( dz_stretch_level_start ) < dz(1) * 1.5_wp ) )  THEN
1835       WRITE( message, * ) 'Eeach dz_stretch_level_start has to be ',&
1836                           'larger than ', dz(1) * 1.5
1837          CALL inifor_abort( 'stretched_z', message)
1838    ENDIF
1839
1840!
1841!-- The stretching has to start and end on a grid level. Therefore
1842!-- user-specified values have to ''interpolate'' to the next lowest level
1843    IF ( number_stretch_level_start /= 0 )  THEN
1844       dz_stretch_level_start(1) = INT( (dz_stretch_level_start(1) -        &
1845                                         dz(1)/2.0) / dz(1) )               &
1846                                   * dz(1) + dz(1)/2.0
1847    ENDIF
1848   
1849    IF ( number_stretch_level_start > 1 )  THEN
1850       DO  n = 2, number_stretch_level_start
1851          dz_stretch_level_start(n) = INT( dz_stretch_level_start(n) /      &
1852                                           dz(n) ) * dz(n)
1853       ENDDO
1854    ENDIF
1855   
1856    IF ( number_stretch_level_end /= 0 )  THEN
1857       DO  n = 1, number_stretch_level_end
1858          dz_stretch_level_end(n) = INT( dz_stretch_level_end(n) /          &
1859                                         dz(n+1) ) * dz(n+1)
1860       ENDDO
1861    ENDIF
1862 
1863!
1864!-- Determine stretching factor if necessary
1865    IF ( number_stretch_level_end >= 1 )  THEN
1866       CALL calculate_stretching_factor( number_stretch_level_end, dz,      & 
1867                                         dz_stretch_factor_array,           &   
1868                                         dz_stretch_level_end,              &
1869                                         dz_stretch_level_start )
1870    ENDIF
1871
1872    z(1) = dz(1) * 0.5_wp
1873!
1874    dz_stretch_level_index = n
1875    dz_stretched = dz(1)
1876    DO  k = 2, n
1877
1878       IF ( dz_stretch_level <= z(k-1)  .AND.  dz_stretched < dz_max )  THEN
1879
1880          dz_stretched = dz_stretched * dz_stretch_factor
1881          dz_stretched = MIN( dz_stretched, dz_max )
1882
1883          IF ( dz_stretch_level_index == n )  dz_stretch_level_index = k-1
1884
1885       ENDIF
1886
1887       z(k) = z(k-1) + dz_stretched
1888
1889    ENDDO
1890!-- Determine u and v height levels considering the possibility of grid
1891!-- stretching in several heights.
1892    n = 1
1893    dz_stretch_level_start_index(:) = UBOUND(z, 1)
1894    dz_stretch_level_end_index(:) = UBOUND(z, 1)
1895    dz_stretched = dz(1)
1896
1897!-- The default value of dz_stretch_level_start is negative, thus the first
1898!-- condition is always true. Hence, the second condition is necessary.
1899    DO  k = 2, UBOUND(z, 1)
1900       IF ( dz_stretch_level_start(n) <= z(k-1) .AND.                      &
1901            dz_stretch_level_start(n) /= -9999999.9_wp )  THEN
1902          dz_stretched = dz_stretched * dz_stretch_factor_array(n)
1903         
1904          IF ( dz(n) > dz(n+1) )  THEN
1905             dz_stretched = MAX( dz_stretched, dz(n+1) ) !Restrict dz_stretched to the user-specified (higher) dz
1906          ELSE
1907             dz_stretched = MIN( dz_stretched, dz(n+1) ) !Restrict dz_stretched to the user-specified (lower) dz
1908          ENDIF
1909         
1910          IF ( dz_stretch_level_start_index(n) == UBOUND(z, 1) )            &
1911          dz_stretch_level_start_index(n) = k-1
1912         
1913       ENDIF
1914       
1915       z(k) = z(k-1) + dz_stretched
1916       
1917!
1918!--    Make sure that the stretching ends exactly at dz_stretch_level_end
1919       dz_level_end = ABS( z(k) - dz_stretch_level_end(n) ) 
1920       
1921       IF ( dz_level_end < dz(n+1)/3.0 )  THEN
1922          z(k) = dz_stretch_level_end(n)
1923          dz_stretched = dz(n+1)
1924          dz_stretch_level_end_index(n) = k
1925          n = n + 1             
1926       ENDIF
1927    ENDDO
1928
1929    DEALLOCATE( min_dz_stretch_level_end )
1930
1931 END SUBROUTINE stretched_z
1932
1933
1934!------------------------------------------------------------------------------!
1935! Description: [PALM subroutine]
1936! -----------------------------------------------------------------------------!
1937!> Calculation of the stretching factor through an iterative method. Ideas were
1938!> taken from the paper "Regional stretched grid generation and its application
1939!> to the NCAR RegCM (1999)". Normally, no analytic solution exists because the
1940!> system of equations has two variables (r,l) but four requirements
1941!> (l=integer, r=[0,88;1,2], Eq(6), Eq(5) starting from index j=1) which
1942!> results into an overdetermined system.
1943!------------------------------------------------------------------------------!
1944 SUBROUTINE calculate_stretching_factor( number_end, dz,                       &
1945                                         dz_stretch_factor_array,              &   
1946                                         dz_stretch_level_end,                 &
1947                                         dz_stretch_level_start )
1948 
1949    REAL(wp), DIMENSION(:), INTENT(IN)    ::  dz
1950    REAL(wp), DIMENSION(:), INTENT(INOUT) ::  dz_stretch_factor_array
1951    REAL(wp), DIMENSION(:), INTENT(IN)    ::  dz_stretch_level_end, dz_stretch_level_start
1952 
1953    INTEGER ::  iterations  !< number of iterations until stretch_factor_lower/upper_limit is reached 
1954    INTEGER ::  l_rounded   !< after l_rounded grid levels dz(n) is strechted to dz(n+1) with stretch_factor_2
1955    INTEGER ::  n           !< loop variable for stretching
1956   
1957    INTEGER, INTENT(IN) ::  number_end !< number of user-specified end levels for stretching
1958       
1959    REAL(wp) ::  delta_l               !< absolute difference between l and l_rounded
1960    REAL(wp) ::  delta_stretch_factor  !< absolute difference between stretch_factor_1 and stretch_factor_2
1961    REAL(wp) ::  delta_total_new       !< sum of delta_l and delta_stretch_factor for the next iteration (should be as small as possible)
1962    REAL(wp) ::  delta_total_old       !< sum of delta_l and delta_stretch_factor for the last iteration
1963    REAL(wp) ::  distance              !< distance between dz_stretch_level_start and dz_stretch_level_end (stretching region)
1964    REAL(wp) ::  l                     !< value that fulfil Eq. (5) in the paper mentioned above together with stretch_factor_1 exactly
1965    REAL(wp) ::  numerator             !< numerator of the quotient
1966    REAL(wp) ::  stretch_factor_1      !< stretching factor that fulfil Eq. (5) togehter with l exactly
1967    REAL(wp) ::  stretch_factor_2      !< stretching factor that fulfil Eq. (6) togehter with l_rounded exactly
1968   
1969    REAL(wp) ::  dz_stretch_factor_array_2(9) = 1.08_wp  !< Array that contains all stretch_factor_2 that belongs to stretch_factor_1
1970   
1971    REAL(wp), PARAMETER ::  stretch_factor_interval = 1.0E-06  !< interval for sampling possible stretching factors
1972    REAL(wp), PARAMETER ::  stretch_factor_lower_limit = 0.88  !< lowest possible stretching factor
1973    REAL(wp), PARAMETER ::  stretch_factor_upper_limit = 1.12  !< highest possible stretching factor
1974 
1975 
1976    l = 0
1977    DO  n = 1, number_end
1978   
1979       iterations = 1
1980       stretch_factor_1 = 1.0 
1981       stretch_factor_2 = 1.0
1982       delta_total_old = 1.0
1983       
1984       IF ( dz(n) > dz(n+1) )  THEN
1985          DO  WHILE ( stretch_factor_1 >= stretch_factor_lower_limit ) 
1986             
1987             stretch_factor_1 = 1.0 - iterations * stretch_factor_interval
1988             distance = ABS( dz_stretch_level_end(n) -                   &
1989                        dz_stretch_level_start(n) ) 
1990             numerator = distance*stretch_factor_1/dz(n) +               &
1991                         stretch_factor_1 - distance/dz(n)
1992             
1993             IF ( numerator > 0.0 )  THEN
1994                l = LOG( numerator ) / LOG( stretch_factor_1 ) - 1.0
1995                l_rounded = NINT( l )
1996                delta_l = ABS( l_rounded - l ) / l
1997             ENDIF
1998             
1999             stretch_factor_2 = EXP( LOG( dz(n+1)/dz(n) ) / (l_rounded) )
2000             
2001             delta_stretch_factor = ABS( stretch_factor_1 -              &
2002                                         stretch_factor_2 ) /            &
2003                                    stretch_factor_2
2004             
2005             delta_total_new = delta_l + delta_stretch_factor
2006
2007!
2008!--                stretch_factor_1 is taken to guarantee that the stretching
2009!--                procedure ends as close as possible to dz_stretch_level_end.
2010!--                stretch_factor_2 would guarantee that the stretched dz(n) is
2011!--                equal to dz(n+1) after l_rounded grid levels.
2012             IF (delta_total_new < delta_total_old)  THEN
2013                dz_stretch_factor_array(n) = stretch_factor_1
2014                dz_stretch_factor_array_2(n) = stretch_factor_2
2015                delta_total_old = delta_total_new
2016             ENDIF
2017             
2018             iterations = iterations + 1
2019           
2020          ENDDO
2021             
2022       ELSEIF ( dz(n) < dz(n+1) )  THEN
2023          DO  WHILE ( stretch_factor_1 <= stretch_factor_upper_limit )
2024                     
2025             stretch_factor_1 = 1.0 + iterations * stretch_factor_interval
2026             distance = ABS( dz_stretch_level_end(n) -                      &
2027                        dz_stretch_level_start(n) ) 
2028             numerator = distance*stretch_factor_1/dz(n) +                  &
2029                         stretch_factor_1 - distance/dz(n)
2030             
2031             l = LOG( numerator ) / LOG( stretch_factor_1 ) - 1.0
2032             l_rounded = NINT( l )
2033             delta_l = ABS( l_rounded - l ) / l
2034             
2035             stretch_factor_2 = EXP( LOG( dz(n+1)/dz(n) ) / (l_rounded) )
2036
2037             delta_stretch_factor = ABS( stretch_factor_1 -                 &
2038                                        stretch_factor_2 ) /                &
2039                                        stretch_factor_2
2040             
2041             delta_total_new = delta_l + delta_stretch_factor
2042             
2043!
2044!--          stretch_factor_1 is taken to guarantee that the stretching
2045!--          procedure ends as close as possible to dz_stretch_level_end.
2046!--          stretch_factor_2 would guarantee that the stretched dz(n) is
2047!--          equal to dz(n+1) after l_rounded grid levels.
2048             IF (delta_total_new < delta_total_old)  THEN
2049                dz_stretch_factor_array(n) = stretch_factor_1
2050                dz_stretch_factor_array_2(n) = stretch_factor_2
2051                delta_total_old = delta_total_new
2052             ENDIF
2053             
2054             iterations = iterations + 1
2055          ENDDO
2056         
2057       ELSE
2058          message = 'Two adjacent values of dz must be different'
2059          CALL inifor_abort( 'calculate_stretching_factor', message)
2060       ENDIF
2061
2062!
2063!--    Check if also the second stretching factor fits into the allowed
2064!--    interval. If not, print a warning for the user.
2065       IF ( dz_stretch_factor_array_2(n) < stretch_factor_lower_limit .OR.     & 
2066            dz_stretch_factor_array_2(n) > stretch_factor_upper_limit )  THEN
2067          WRITE( message, * ) 'stretch_factor_2 = ',                           &
2068                                     dz_stretch_factor_array_2(n), ' which is',&
2069                                     ' responsible for exactly reaching& dz =',&
2070                                      dz(n+1), 'after a specific amount of',   & 
2071                                     ' grid levels& exceeds the upper',        &
2072                                     ' limit =', stretch_factor_upper_limit,   &
2073                                     ' &or lower limit = ',                    &
2074                                     stretch_factor_lower_limit
2075          CALL inifor_abort( 'calculate_stretching_factor', message )
2076           
2077       ENDIF
2078    ENDDO
2079       
2080 END SUBROUTINE calculate_stretching_factor
2081
2082!------------------------------------------------------------------------------!
2083! Description:
2084! ------------
2085!> Computes the midpoint between two neighbouring coordinates in the given
2086!> coordinate vector 'z' and stores it in 'zw'.
2087!------------------------------------------------------------------------------!
2088 SUBROUTINE midpoints(z, zw)
2089
2090     REAL(wp), INTENT(IN)  ::  z(0:)
2091     REAL(wp), INTENT(OUT) ::  zw(1:)
2092
2093     INTEGER ::  k
2094
2095     DO  k = 1, UBOUND(zw, 1)
2096        zw(k) = 0.5_wp * (z(k-1) + z(k))
2097     ENDDO
2098
2099 END SUBROUTINE midpoints
2100
2101!------------------------------------------------------------------------------!
2102! Description:
2103! ------------
2104!> Defines INFOR's IO groups.
2105!------------------------------------------------------------------------------!
2106 SUBROUTINE setup_io_groups()
2107
2108    INTEGER ::  ngroups
2109
2110    ngroups = 16
2111    ALLOCATE( io_group_list(ngroups) )
2112
2113!
2114!-- soil temp
2115    io_group_list(1) = init_io_group(                                       &
2116       in_files = soil_files,                                               &
2117       out_vars = output_var_table(1:1),                                    &
2118       in_var_list = input_var_table(1:1),                                  &
2119       kind = 'soil-temperature'                                            &
2120    )
2121
2122!
2123!-- soil water
2124    io_group_list(2) = init_io_group(                                       &
2125       in_files = soil_files,                                               &
2126       out_vars = output_var_table(2:2),                                    &
2127       in_var_list = input_var_table(2:2),                                  &
2128       kind = 'soil-water'                                                  &
2129    )
2130
2131!
2132!-- potential temperature, surface pressure, specific humidity including
2133!-- nudging and subsidence, and geostrophic winds ug, vg
2134    io_group_list(3) = init_io_group(                                       &
2135       in_files = flow_files,                                               &
2136       out_vars = [output_var_table(56:64),                                 & ! internal averaged density and pressure profiles
2137                   output_var_table(3:8), output_var_table(9:14),           &
2138                   output_var_table(42:42), output_var_table(43:44),        &
2139                   output_var_table(49:51), output_var_table(52:54)],       &
2140       in_var_list = (/input_var_table(3), input_var_table(17),             & ! T, P, QV
2141                       input_var_table(4) /),                               &
2142       kind = 'thermodynamics',                                             &
2143       n_output_quantities = 4                                              & ! P, Theta, Rho, qv
2144    )
2145
2146!
2147!-- Moved to therodynamic io_group
2148    !io_group_list(4) = init_io_group(                                       &
2149    !   in_files = flow_files,                                               &
2150    !   out_vars = [output_var_table(9:14), output_var_table(52:54)],        &
2151    !   in_var_list = input_var_table(4:4),                                  &
2152    !   kind = 'scalar'                                                      &
2153    !)
2154
2155!
2156!-- u and v velocity
2157    io_group_list(5) = init_io_group(                                       &
2158       in_files = flow_files,                                               &
2159       out_vars = [output_var_table(15:26), output_var_table(45:46)],       &
2160       !out_vars = output_var_table(15:20),                                  &
2161       in_var_list = input_var_table(5:6),                                  &
2162       !in_var_list = input_var_table(5:5),                                  &
2163       kind = 'velocities'                                                  &
2164    )
2165   
2166!
2167!-- v velocity, deprecated!
2168    !io_group_list(6) = init_io_group(                                       &
2169    !   in_files = flow_files,                                               &
2170    !   out_vars = output_var_table(21:26),                                  &
2171    !   in_var_list = input_var_table(6:6),                                  &
2172    !   kind = 'horizontal velocity'                                         &
2173    !)
2174    !io_group_list(6)%to_be_processed = .FALSE.
2175   
2176!
2177!-- w velocity and subsidence and w nudging
2178    io_group_list(7) = init_io_group(                                       &
2179       in_files = flow_files,                                               &
2180       out_vars = [output_var_table(27:32), output_var_table(47:48)],       &
2181       in_var_list = input_var_table(7:7),                                  &
2182       kind = 'scalar'                                                      &
2183    )
2184!
2185!-- rain
2186    io_group_list(8) = init_io_group(                                       &
2187       in_files = soil_moisture_files,                                      &
2188       out_vars = output_var_table(33:33),                                  &
2189       in_var_list = input_var_table(8:8),                                  &
2190       kind = 'accumulated'                                                 &
2191    )
2192    io_group_list(8)%to_be_processed = .FALSE.
2193!
2194!-- snow
2195    io_group_list(9) = init_io_group(                                       &
2196       in_files = soil_moisture_files,                                      &
2197       out_vars = output_var_table(34:34),                                  &
2198       in_var_list = input_var_table(9:9),                                  &
2199       kind = 'accumulated'                                                 &
2200    )
2201    io_group_list(9)%to_be_processed = .FALSE.
2202!
2203!-- graupel
2204    io_group_list(10) = init_io_group(                                      &
2205       in_files = soil_moisture_files,                                      &
2206       out_vars = output_var_table(35:35),                                  &
2207       in_var_list = input_var_table(10:10),                                &
2208       kind = 'accumulated'                                                 &
2209    )
2210    io_group_list(10)%to_be_processed = .FALSE.
2211!
2212!-- evapotranspiration
2213    io_group_list(11) = init_io_group(                                      &
2214       in_files = soil_moisture_files,                                      &
2215       out_vars = output_var_table(37:37),                                  &
2216       in_var_list = input_var_table(11:11),                                &
2217       kind = 'accumulated'                                                 &
2218    )
2219    io_group_list(11)%to_be_processed = .FALSE.
2220!
2221!-- 2m air temperature
2222    io_group_list(12) = init_io_group(                                      &
2223       in_files = soil_moisture_files,                                      &
2224       out_vars = output_var_table(36:36),                                  &
2225       in_var_list = input_var_table(12:12),                                &
2226       kind = 'surface'                                                     &
2227    )
2228    io_group_list(12)%to_be_processed = .FALSE.
2229!
2230!-- incoming diffusive sw flux
2231    io_group_list(13) = init_io_group(                                      &
2232       in_files = radiation_files,                                          &
2233       out_vars = output_var_table(38:38),                                  &
2234       in_var_list = input_var_table(13:13),                                &
2235       kind = 'running average'                                             &
2236    )
2237    io_group_list(13)%to_be_processed = .FALSE.
2238!
2239!-- incoming direct sw flux
2240    io_group_list(14) = init_io_group(                                      &
2241       in_files = radiation_files,                                          &
2242       out_vars = output_var_table(39:39),                                  &
2243       in_var_list = input_var_table(14:14),                                &
2244       kind = 'running average'                                             &
2245    )
2246    io_group_list(14)%to_be_processed = .FALSE.
2247!
2248!-- sw radiation balance
2249    io_group_list(15) = init_io_group(                                      &
2250       in_files = radiation_files,                                          &
2251       out_vars = output_var_table(40:40),                                  &
2252       in_var_list = input_var_table(15:15),                                &
2253       kind = 'running average'                                             &
2254    )
2255    io_group_list(15)%to_be_processed = .FALSE.
2256!
2257!-- lw radiation balance
2258    io_group_list(16) = init_io_group(                                      &
2259       in_files = radiation_files,                                          &
2260       out_vars = output_var_table(41:41),                                  &
2261       in_var_list = input_var_table(16:16),                                &
2262       kind = 'running average'                                             &
2263    )
2264    io_group_list(16)%to_be_processed = .FALSE.
2265
2266 END SUBROUTINE setup_io_groups
2267
2268
2269!------------------------------------------------------------------------------!
2270! Description:
2271! ------------
2272!> Initializes an IO group with the list of input files, the lists of
2273!> input and output variables, and sets the number of output quantities.
2274!>
2275!> In this context, output quantities refers to the physical quantities required
2276!> to interpolate netCDF output variables, not the output variables itsleft. For
2277!> instance, the output variables ls_forcing_left/_right/_north/_south all rely
2278!> on the output quantity Theta.
2279!------------------------------------------------------------------------------!
2280 FUNCTION init_io_group(in_files, out_vars, in_var_list, kind,              &
2281                        n_output_quantities) RESULT(group)
2282    CHARACTER(LEN=PATH), INTENT(IN) ::  in_files(:)
2283    CHARACTER(LEN=*), INTENT(IN)    ::  kind
2284    TYPE(nc_var), INTENT(IN)        ::  out_vars(:)
2285    TYPE(nc_var), INTENT(IN)        ::  in_var_list(:)
2286    INTEGER, OPTIONAL               ::  n_output_quantities
2287
2288    TYPE(io_group)                  ::  group
2289
2290    group%nt = SIZE(in_files)
2291    group%nv = SIZE(out_vars)
2292    group%n_inputs = SIZE(in_var_list)
2293    group%kind = TRIM(kind)
2294!
2295!-- For the 'thermodynamics' IO group, one quantity more than input variables
2296!-- is needed to compute all output variables of the IO group. Concretely, in
2297!-- preprocess() the density is computed from T,P or PP,QV in adddition to
2298!-- the variables Theta, p, qv. In read_input_variables(),
2299!-- n_output_quantities is used to allocate the correct number of input
2300!-- buffers.
2301    IF ( PRESENT(n_output_quantities) )  THEN
2302       group%n_output_quantities = n_output_quantities
2303    ELSE
2304       group%n_output_quantities = group%n_inputs
2305    ENDIF
2306
2307    ALLOCATE(group%in_var_list(group%n_inputs))
2308    ALLOCATE(group%in_files(group%nt))
2309    ALLOCATE(group%out_vars(group%nv))
2310
2311    group%in_var_list = in_var_list
2312    group%in_files = in_files
2313    group%out_vars = out_vars
2314    group%to_be_processed = .TRUE.
2315
2316 END FUNCTION init_io_group
2317
2318
2319!------------------------------------------------------------------------------!
2320! Description:
2321! ------------
2322!> Deallocates all allocated variables.
2323!------------------------------------------------------------------------------!
2324 SUBROUTINE fini_grids()
2325
2326    CALL report('fini_grids', 'Deallocating grids', cfg%debug)
2327   
2328    DEALLOCATE(x, y, z, xu, yv, zw, z_column, zw_column)
2329
2330    DEALLOCATE(palm_grid%x,  palm_grid%y,  palm_grid%z,                        &
2331               palm_grid%xu, palm_grid%yv, palm_grid%zw,                       &
2332               palm_grid%clon,  palm_grid%clat,                                &
2333               palm_grid%clonu, palm_grid%clatu)
2334
2335    DEALLOCATE(palm_intermediate%x,  palm_intermediate%y,  palm_intermediate%z, &
2336               palm_intermediate%xu, palm_intermediate%yv, palm_intermediate%zw,&
2337               palm_intermediate%clon,  palm_intermediate%clat,                & 
2338               palm_intermediate%clonu, palm_intermediate%clatu)
2339
2340    DEALLOCATE(cosmo_grid%lon,  cosmo_grid%lat,                                &
2341               cosmo_grid%lonu, cosmo_grid%latv,                               &
2342               cosmo_grid%hfl)
2343
2344 END SUBROUTINE fini_grids
2345
2346
2347!------------------------------------------------------------------------------!
2348! Description:
2349! ------------
2350!> Initializes the variable list.
2351!------------------------------------------------------------------------------!
2352 SUBROUTINE setup_variable_tables(ic_mode)
2353    CHARACTER(LEN=*), INTENT(IN) ::  ic_mode
2354    INTEGER                      ::  n_invar = 0  !< number of variables in the input variable table
2355    INTEGER                      ::  n_outvar = 0 !< number of variables in the output variable table
2356    TYPE(nc_var), POINTER        ::  var
2357
2358    IF (TRIM(cfg%start_date) == '')  THEN
2359       message = 'Simulation start date has not been set.'
2360       CALL inifor_abort('setup_variable_tables', message)
2361    ENDIF
2362
2363    nc_source_text = 'COSMO-DE analysis from ' // TRIM(cfg%start_date)
2364
2365    n_invar = 17
2366    n_outvar = 64
2367    ALLOCATE( input_var_table(n_invar) )
2368    ALLOCATE( output_var_table(n_outvar) )
2369
2370!
2371!------------------------------------------------------------------------------
2372!- Section 1: NetCDF input variables
2373!------------------------------------------------------------------------------
2374    var => input_var_table(1)
2375    var%name = 'T_SO'
2376    var%to_be_processed = .TRUE.
2377    var%is_upside_down = .FALSE.
2378
2379    var => input_var_table(2)
2380    var%name = 'W_SO'
2381    var%to_be_processed = .TRUE.
2382    var%is_upside_down = .FALSE.
2383
2384    var => input_var_table(3)
2385    var%name = 'T'
2386    var%to_be_processed = .TRUE.
2387    var%is_upside_down = .TRUE.
2388
2389    var => input_var_table(4)
2390    var%name = 'QV'
2391    var%to_be_processed = .TRUE.
2392    var%is_upside_down = .TRUE.
2393
2394    var => input_var_table(5)
2395    var%name = 'U'
2396    var%to_be_processed = .TRUE.
2397    var%is_upside_down = .TRUE.
2398
2399    var => input_var_table(6)
2400    var%name = 'V'
2401    var%to_be_processed = .TRUE.
2402    var%is_upside_down = .TRUE.
2403
2404    var => input_var_table(7)
2405    var%name = 'W'
2406    var%to_be_processed = .TRUE.
2407    var%is_upside_down = .TRUE.
2408
2409    var => input_var_table(8)
2410    var%name = 'RAIN_GSP'
2411    var%to_be_processed = .TRUE.
2412    var%is_upside_down = .FALSE.
2413
2414    var => input_var_table(9)
2415    var%name = 'SNOW_GSP'
2416    var%to_be_processed = .TRUE.
2417    var%is_upside_down = .FALSE.
2418
2419    var => input_var_table(10)
2420    var%name = 'GRAU_GSP'
2421    var%to_be_processed = .TRUE.
2422    var%is_upside_down = .FALSE.
2423
2424    var => input_var_table(11)
2425    var%name = 'AEVAP_S'
2426    var%to_be_processed = .TRUE.
2427    var%is_upside_down = .FALSE.
2428
2429    var => input_var_table(12)
2430    var%name = 'T_2M'
2431    var%to_be_processed = .TRUE.
2432    var%is_upside_down = .FALSE.
2433
2434    var => input_var_table(13)
2435    var%name = 'ASWDIFD_S'
2436    var%to_be_processed = .TRUE.
2437    var%is_upside_down = .FALSE.
2438
2439    var => input_var_table(14)
2440    var%name = 'ASWDIR_S'
2441    var%to_be_processed = .TRUE.
2442    var%is_upside_down = .FALSE.
2443
2444    var => input_var_table(15)
2445    var%name = 'ASOB_S'
2446    var%to_be_processed = .TRUE.
2447    var%is_upside_down = .FALSE.
2448
2449    var => input_var_table(16)
2450    var%name = 'ATHB_S'
2451    var%to_be_processed = .TRUE.
2452    var%is_upside_down = .FALSE.
2453
2454    var => input_var_table(17)
2455    var%name = 'P'
2456    var%to_be_processed = .TRUE.
2457    var%is_upside_down = .TRUE.
2458
2459!
2460!------------------------------------------------------------------------------
2461!- Section 2: NetCDF output variables
2462!------------------------------------------------------------------------------
2463!
2464!------------------------------------------------------------------------------
2465! Section 2.1: Realistic forcings, i.e. 3D initial and boundary conditions
2466!------------------------------------------------------------------------------
2467    output_var_table(1) = init_nc_var(                                      &
2468       name              = 'init_soil_t',                                   &
2469       std_name          = "",                                              &
2470       long_name         = "initial soil temperature",                      &
2471       units             = "K",                                             &
2472       kind              = "init soil",                                     &
2473       input_id          = 1,                                               &
2474       output_file       = output_file,                                     &
2475       grid              = palm_grid,                                       &
2476       intermediate_grid = palm_intermediate                                &
2477    )
2478
2479    output_var_table(2) = init_nc_var(                                      &
2480       name              = 'init_soil_m',                                   &
2481       std_name          = "",                                              &
2482       long_name         = "initial soil moisture",                         &
2483       units             = "m^3/m^3",                                       &
2484       kind              = "init soil",                                     &
2485       input_id          = 1,                                               &
2486       output_file       = output_file,                                     &
2487       grid              = palm_grid,                                       &
2488       intermediate_grid = palm_intermediate                                &
2489    )
2490
2491    output_var_table(3) = init_nc_var(                                      &
2492       name              = 'init_atmosphere_pt',                            &
2493       std_name          = "",                                              &
2494       long_name         = "initial potential temperature",                 &
2495       units             = "K",                                             &
2496       kind              = "init scalar",                                   &
2497       input_id          = 1,                                               & ! first in (T, p) IO group
2498       output_file       = output_file,                                     &
2499       grid              = palm_grid,                                       &
2500       intermediate_grid = palm_intermediate,                               &
2501       is_profile = (TRIM(ic_mode) == 'profile')                            &
2502    )
2503    IF (TRIM(ic_mode) == 'profile')  THEN
2504       output_var_table(3)%averaging_grid => averaged_initial_scalar_profile
2505    ENDIF
2506
2507    output_var_table(4) = init_nc_var(                                      &
2508       name              = 'ls_forcing_left_pt',                            &
2509       std_name          = "",                                              &
2510       long_name         = "large-scale forcing for left model boundary for the potential temperature", &
2511       units             = "K",                                             &
2512       kind              = "left scalar",                                   &
2513       input_id          = 1,                                               &
2514       grid              = scalars_west_grid,                               &
2515       intermediate_grid = scalars_west_intermediate,                       &
2516       output_file = output_file                                            &
2517    )
2518
2519    output_var_table(5) = init_nc_var(                                      &
2520       name              = 'ls_forcing_right_pt',                           &
2521       std_name          = "",                                              &
2522       long_name         = "large-scale forcing for right model boundary for the potential temperature", &
2523       units             = "K",                                             &
2524       kind              = "right scalar",                                  &
2525       input_id          = 1,                                               &
2526       grid              = scalars_east_grid,                               &
2527       intermediate_grid = scalars_east_intermediate,                       &
2528       output_file = output_file                                            &
2529    )
2530
2531    output_var_table(6) = init_nc_var(                                      &
2532       name              = 'ls_forcing_north_pt',                           &
2533       std_name          = "",                                              &
2534       long_name         = "large-scale forcing for north model boundary for the potential temperature", &
2535       units             = "K",                                             &
2536       kind              = "north scalar",                                  &
2537       input_id          = 1,                                               &
2538       grid              = scalars_north_grid,                              &
2539       intermediate_grid = scalars_north_intermediate,                      &
2540       output_file = output_file                                            &
2541    )
2542
2543    output_var_table(7) = init_nc_var(                                      &
2544       name              = 'ls_forcing_south_pt',                           &
2545       std_name          = "",                                              &
2546       long_name         = "large-scale forcing for south model boundary for the potential temperature", &
2547       units             = "K",                                             &
2548       kind              = "south scalar",                                  &
2549       input_id          = 1,                                               &
2550       grid              = scalars_south_grid,                              &
2551       intermediate_grid = scalars_south_intermediate,                      &
2552       output_file = output_file                                            &
2553    )
2554
2555    output_var_table(8) = init_nc_var(                                      &
2556       name              = 'ls_forcing_top_pt',                             &
2557       std_name          = "",                                              &
2558       long_name         = "large-scale forcing for top model boundary for the potential temperature", &
2559       units             = "K",                                             &
2560       kind              = "top scalar",                                    &
2561       input_id          = 1,                                               &
2562       grid              = scalars_top_grid,                                &
2563       intermediate_grid = scalars_top_intermediate,                        &
2564       output_file = output_file                                            &
2565    )
2566
2567    output_var_table(9) = init_nc_var(                                      &
2568       name              = 'init_atmosphere_qv',                            &
2569       std_name          = "",                                              &
2570       long_name         = "initial specific humidity",                     &
2571       units             = "kg/kg",                                         &
2572       kind              = "init scalar",                                   &
2573       input_id          = 3,                                               &
2574       output_file       = output_file,                                     &
2575       grid              = palm_grid,                                       &
2576       intermediate_grid = palm_intermediate,                               &
2577       is_profile = (TRIM(ic_mode) == 'profile')                            &
2578    )
2579    IF (TRIM(ic_mode) == 'profile')  THEN
2580       output_var_table(9)%averaging_grid => averaged_initial_scalar_profile
2581    ENDIF
2582
2583    output_var_table(10) = init_nc_var(                                     &
2584       name              = 'ls_forcing_left_qv',                            &
2585       std_name          = "",                                              &
2586       long_name         = "large-scale forcing for left model boundary for the specific humidity", &
2587       units             = "kg/kg",                                         &
2588       kind              = "left scalar",                                   &
2589       input_id          = 3,                                               &
2590       output_file       = output_file,                                     &
2591       grid              = scalars_west_grid,                               &
2592       intermediate_grid = scalars_west_intermediate                        &
2593    )
2594
2595    output_var_table(11) = init_nc_var(                                     &
2596       name              = 'ls_forcing_right_qv',                           &
2597       std_name          = "",                                              &
2598       long_name         = "large-scale forcing for right model boundary for the specific humidity", &
2599       units             = "kg/kg",                                         &
2600       kind              = "right scalar",                                  &
2601       input_id          = 3,                                               &
2602       output_file       = output_file,                                     &
2603       grid              = scalars_east_grid,                               &
2604       intermediate_grid = scalars_east_intermediate                        &
2605    )
2606
2607    output_var_table(12) = init_nc_var(                                     &
2608       name              = 'ls_forcing_north_qv',                           &
2609       std_name          = "",                                              &
2610       long_name         = "large-scale forcing for north model boundary for the specific humidity", &
2611       units             = "kg/kg",                                         &
2612       kind              = "north scalar",                                  &
2613       input_id          = 3,                                               &
2614       output_file       = output_file,                                     &
2615       grid              = scalars_north_grid,                              &
2616       intermediate_grid = scalars_north_intermediate                       &
2617    )
2618
2619    output_var_table(13) = init_nc_var(                                     &
2620       name              = 'ls_forcing_south_qv',                           &
2621       std_name          = "",                                              &
2622       long_name         = "large-scale forcing for south model boundary for the specific humidity", &
2623       units             = "kg/kg",                                         &
2624       kind              = "south scalar",                                  &
2625       input_id          = 3,                                               &
2626       output_file       = output_file,                                     &
2627       grid              = scalars_south_grid,                              &
2628       intermediate_grid = scalars_south_intermediate                       &
2629    )
2630
2631    output_var_table(14) = init_nc_var(                                     &
2632       name              = 'ls_forcing_top_qv',                             &
2633       std_name          = "",                                              &
2634       long_name         = "large-scale forcing for top model boundary for the specific humidity", &
2635       units             = "kg/kg",                                         &
2636       kind              = "top scalar",                                    &
2637       input_id          = 3,                                               &
2638       output_file       = output_file,                                     &
2639       grid              = scalars_top_grid,                                &
2640       intermediate_grid = scalars_top_intermediate                         &
2641    )
2642
2643    output_var_table(15) = init_nc_var(                                     &
2644       name              = 'init_atmosphere_u',                             &
2645       std_name          = "",                                              &
2646       long_name         = "initial wind component in x direction",         &
2647       units             = "m/s",                                           &
2648       kind              = "init u",                                        &
2649       input_id          = 1,                                               & ! first in (U, V) I/O group
2650       output_file       = output_file,                                     &
2651       grid              = u_initial_grid,                                  &
2652       intermediate_grid = u_initial_intermediate,                          &
2653       is_profile = (TRIM(ic_mode) == 'profile')                            &
2654    )
2655    IF (TRIM(ic_mode) == 'profile')  THEN
2656       output_var_table(15)%averaging_grid => averaged_initial_scalar_profile
2657    ENDIF
2658
2659    output_var_table(16) = init_nc_var(                                     &
2660       name              = 'ls_forcing_left_u',                             &
2661       std_name          = "",                                              &
2662       long_name         = "large-scale forcing for left model boundary for the wind component in x direction", &
2663       units             = "m/s",                                           &
2664       kind              = "left u",                                        &
2665       input_id          = 1,                                               & ! first in (U, V) I/O group
2666       output_file       = output_file,                                     &
2667       grid              = u_west_grid,                                     &
2668       intermediate_grid = u_west_intermediate                              &
2669    )
2670
2671    output_var_table(17) = init_nc_var(                                     &
2672       name              = 'ls_forcing_right_u',                            &
2673       std_name          = "",                                              &
2674       long_name         = "large-scale forcing for right model boundary for the wind component in x direction", &
2675       units             = "m/s",                                           &
2676       kind              = "right u",                                       &
2677       input_id          = 1,                                               & ! first in (U, V) I/O group
2678       output_file       = output_file,                                     &
2679       grid              = u_east_grid,                                     &
2680       intermediate_grid = u_east_intermediate                              &
2681    )
2682
2683    output_var_table(18) = init_nc_var(                                     &
2684       name              = 'ls_forcing_north_u',                            &
2685       std_name          = "",                                              &
2686       long_name         = "large-scale forcing for north model boundary for the wind component in x direction", &
2687       units             = "m/s",                                           &
2688       kind              = "north u",                                       &
2689       input_id          = 1,                                               & ! first in (U, V) I/O group
2690       output_file       = output_file,                                     &
2691       grid              = u_north_grid,                                    &
2692       intermediate_grid = u_north_intermediate                             &
2693    )
2694
2695    output_var_table(19) = init_nc_var(                                     &
2696       name              = 'ls_forcing_south_u',                            &
2697       std_name          = "",                                              &
2698       long_name         = "large-scale forcing for south model boundary for the wind component in x direction", &
2699       units             = "m/s",                                           &
2700       kind              = "south u",                                       &
2701       input_id          = 1,                                               & ! first in (U, V) I/O group
2702       output_file       = output_file,                                     &
2703       grid              = u_south_grid,                                    &
2704       intermediate_grid = u_south_intermediate                             &
2705    )
2706
2707    output_var_table(20) = init_nc_var(                                     &
2708       name              = 'ls_forcing_top_u',                              &
2709       std_name          = "",                                              &
2710       long_name         = "large-scale forcing for top model boundary for the wind component in x direction", &
2711       units             = "m/s",                                           &
2712       kind              = "top u",                                         &
2713       input_id          = 1,                                               & ! first in (U, V) I/O group
2714       output_file       = output_file,                                     &
2715       grid              = u_top_grid,                                      &
2716       intermediate_grid = u_top_intermediate                               &
2717    )
2718
2719    output_var_table(21) = init_nc_var(                                     &
2720       name              = 'init_atmosphere_v',                             &
2721       std_name          = "",                                              &
2722       long_name         = "initial wind component in y direction",         &
2723       units             = "m/s",                                           &
2724       kind              = "init v",                                        &
2725       input_id          = 2,                                               & ! second in (U, V) I/O group
2726       output_file       = output_file,                                     &
2727       grid              = v_initial_grid,                                  &
2728       intermediate_grid = v_initial_intermediate,                          &
2729       is_profile = (TRIM(ic_mode) == 'profile')                            &
2730    )
2731    IF (TRIM(ic_mode) == 'profile')  THEN
2732       output_var_table(21)%averaging_grid => averaged_initial_scalar_profile
2733    ENDIF
2734
2735    output_var_table(22) = init_nc_var(                                     &
2736       name              = 'ls_forcing_left_v',                             &
2737       std_name          = "",                                              &
2738       long_name         = "large-scale forcing for left model boundary for the wind component in y direction", &
2739       units             = "m/s",                                           &
2740       kind              = "right v",                                       &
2741       input_id          = 2,                                               & ! second in (U, V) I/O group
2742       output_file       = output_file,                                     &
2743       grid              = v_west_grid,                                     &
2744       intermediate_grid = v_west_intermediate                              &
2745    )
2746
2747    output_var_table(23) = init_nc_var(                                     &
2748       name              = 'ls_forcing_right_v',                            &
2749       std_name          = "",                                              &
2750       long_name         = "large-scale forcing for right model boundary for the wind component in y direction", &
2751       units             = "m/s",                                           &
2752       kind              = "right v",                                       &
2753       input_id          = 2,                                               & ! second in (U, V) I/O group
2754       output_file       = output_file,                                     &
2755       grid              = v_east_grid,                                     &
2756       intermediate_grid = v_east_intermediate                              &
2757    )
2758
2759    output_var_table(24) = init_nc_var(                                     &
2760       name              = 'ls_forcing_north_v',                            &
2761       std_name          = "",                                              &
2762       long_name         = "large-scale forcing for north model boundary for the wind component in y direction", &
2763       units             = "m/s",                                           &
2764       kind              = "north v",                                       &
2765       input_id          = 2,                                               & ! second in (U, V) I/O group
2766       output_file       = output_file,                                     &
2767       grid              = v_north_grid,                                    &
2768       intermediate_grid = v_north_intermediate                             &
2769    )
2770
2771    output_var_table(25) = init_nc_var(                                     &
2772       name              = 'ls_forcing_south_v',                            &
2773       std_name          = "",                                              &
2774       long_name         = "large-scale forcing for south model boundary for the wind component in y direction", &
2775       units             = "m/s",                                           &
2776       kind              = "south v",                                       &
2777       input_id          = 2,                                               & ! second in (U, V) I/O group
2778       output_file       = output_file,                                     &
2779       grid              = v_south_grid,                                    &
2780       intermediate_grid = v_south_intermediate                             &
2781    )
2782
2783    output_var_table(26) = init_nc_var(                                     &
2784       name              = 'ls_forcing_top_v',                              &
2785       std_name          = "",                                              &
2786       long_name         = "large-scale forcing for top model boundary for the wind component in y direction", &
2787       units             = "m/s",                                           &
2788       kind              = "top v",                                         &
2789       input_id          = 2,                                               & ! second in (U, V) I/O group
2790       output_file       = output_file,                                     &
2791       grid              = v_top_grid,                                      &
2792       intermediate_grid = v_top_intermediate                               &
2793    )
2794
2795    output_var_table(27) = init_nc_var(                                     &
2796       name              = 'init_atmosphere_w',                             &
2797       std_name          = "",                                              &
2798       long_name         = "initial wind component in z direction",         &
2799       units             = "m/s",                                           &
2800       kind              = "init w",                                        &
2801       input_id          = 1,                                               &
2802       output_file       = output_file,                                     &
2803       grid              = w_initial_grid,                                  &
2804       intermediate_grid = w_initial_intermediate,                          &
2805       is_profile = (TRIM(ic_mode) == 'profile')                            &
2806    )
2807    IF (TRIM(ic_mode) == 'profile')  THEN
2808       output_var_table(27)%averaging_grid => averaged_initial_w_profile
2809    ENDIF
2810
2811    output_var_table(28) = init_nc_var(                                     &
2812       name              = 'ls_forcing_left_w',                             &
2813       std_name          = "",                                              &
2814       long_name         = "large-scale forcing for left model boundary for the wind component in z direction", &
2815       units             = "m/s",                                           &
2816       kind              = "left w",                                        &
2817       input_id          = 1,                                               &
2818       output_file       = output_file,                                     &
2819       grid              = w_west_grid,                                     &
2820       intermediate_grid = w_west_intermediate                              &
2821    )
2822
2823    output_var_table(29) = init_nc_var(                                     &
2824       name              = 'ls_forcing_right_w',                            &
2825       std_name          = "",                                              &
2826       long_name         = "large-scale forcing for right model boundary for the wind component in z direction", &
2827       units             = "m/s",                                           &
2828       kind              = "right w",                                       &
2829       input_id          = 1,                                               &
2830       output_file       = output_file,                                     &
2831       grid              = w_east_grid,                                     &
2832       intermediate_grid = w_east_intermediate                              &
2833    )
2834
2835    output_var_table(30) = init_nc_var(                                     &
2836       name              = 'ls_forcing_north_w',                            &
2837       std_name          = "",                                              &
2838       long_name         = "large-scale forcing for north model boundary for the wind component in z direction", &
2839       units             = "m/s",                                           &
2840       kind              = "north w",                                       &
2841       input_id          = 1,                                               &
2842       output_file       = output_file,                                     &
2843       grid              = w_north_grid,                                    &
2844       intermediate_grid = w_north_intermediate                             &
2845    )
2846
2847    output_var_table(31) = init_nc_var(                                     &
2848       name              = 'ls_forcing_south_w',                            &
2849       std_name          = "",                                              &
2850       long_name         = "large-scale forcing for south model boundary for the wind component in z direction", &
2851       units             = "m/s",                                           &
2852       kind              = "south w",                                       &
2853       input_id          = 1,                                               &
2854       output_file       = output_file,                                     &
2855       grid              = w_south_grid,                                    &
2856       intermediate_grid = w_south_intermediate                             &
2857    )
2858
2859    output_var_table(32) = init_nc_var(                                     &
2860       name              = 'ls_forcing_top_w',                              &
2861       std_name          = "",                                              &
2862       long_name         = "large-scale forcing for top model boundary for the wind component in z direction", &
2863       units             = "m/s",                                           &
2864       kind              = "top w",                                         &
2865       input_id          = 1,                                               &
2866       output_file       = output_file,                                     &
2867       grid              = w_top_grid,                                      &
2868       intermediate_grid = w_top_intermediate                               &
2869    )
2870
2871    output_var_table(33) = init_nc_var(                                     &
2872       name              = 'ls_forcing_soil_rain',                          &
2873       std_name          = "",                                              &
2874       long_name         = "large-scale forcing rain",                      &
2875       units             = "kg/m2",                                         &
2876       kind              = "surface forcing",                               &
2877       input_id          = 1,                                               &
2878       output_file       = output_file,                                     &
2879       grid              = palm_grid,                                       &
2880       intermediate_grid = palm_intermediate                                &
2881    )
2882
2883    output_var_table(34) = init_nc_var(                                     &
2884       name              = 'ls_forcing_soil_snow',                          &
2885       std_name          = "",                                              &
2886       long_name         = "large-scale forcing snow",                      &
2887       units             = "kg/m2",                                         &
2888       kind              = "surface forcing",                               &
2889       input_id          = 1,                                               &
2890       output_file       = output_file,                                     &
2891       grid              = palm_grid,                                       &
2892       intermediate_grid = palm_intermediate                                &
2893    )
2894
2895    output_var_table(35) = init_nc_var(                                     &
2896       name              = 'ls_forcing_soil_graupel',                       &
2897       std_name          = "",                                              &
2898       long_name         = "large-scale forcing graupel",                   &
2899       units             = "kg/m2",                                         &
2900       kind              = "surface forcing",                               &
2901       input_id          = 1,                                               &
2902       output_file       = output_file,                                     &
2903       grid              = palm_grid,                                       &
2904       intermediate_grid = palm_intermediate                                &
2905    )
2906
2907    output_var_table(36) = init_nc_var(                                     &
2908       name              = 'ls_forcing_soil_t_2m',                          &
2909       std_name          = "",                                              &
2910       long_name         = "large-scale forcing 2m air temperature",        &
2911       units             = "kg/m2",                                         &
2912       kind              = "surface forcing",                               &
2913       input_id          = 1,                                               &
2914       output_file       = output_file,                                     &
2915       grid              = palm_grid,                                       &
2916       intermediate_grid = palm_intermediate                                &
2917    )
2918
2919    output_var_table(37) = init_nc_var(                                     &
2920       name              = 'ls_forcing_soil_evap',                          &
2921       std_name          = "",                                              &
2922       long_name         = "large-scale forcing evapo-transpiration",       &
2923       units             = "kg/m2",                                         &
2924       kind              = "surface forcing",                               &
2925       input_id          = 1,                                               &
2926       output_file       = output_file,                                     &
2927       grid              = palm_grid,                                       &
2928       intermediate_grid = palm_intermediate                                &
2929    )
2930
2931    output_var_table(38) = init_nc_var(                                     &
2932       name              = 'rad_swd_dif_0',                                 &
2933       std_name          = "",                                              &
2934       long_name         = "incoming diffuse shortwave radiative flux at the surface", &
2935       units             = "W/m2",                                          &
2936       kind              = "surface forcing",                               &
2937       input_id          = 1,                                               &
2938       output_file       = output_file,                                     &
2939       grid              = palm_grid,                                       &
2940       intermediate_grid = palm_intermediate                                &
2941    )
2942
2943    output_var_table(39) = init_nc_var(                                     &
2944       name              = 'rad_swd_dir_0',                                 &
2945       std_name          = "",                                              &
2946       long_name         = "incoming direct shortwave radiative flux at the surface", &
2947       units             = "W/m2",                                          &
2948       kind              = "surface forcing",                               &
2949       input_id          = 1,                                               &
2950       output_file       = output_file,                                     &
2951       grid              = palm_grid,                                       &
2952       intermediate_grid = palm_intermediate                                &
2953    )
2954
2955    output_var_table(40) = init_nc_var(                                     &
2956       name              = 'rad_sw_bal_0',                                  &
2957       std_name          = "",                                              &
2958       long_name         = "shortwave radiation balance at the surface",    &
2959       units             = "W/m2",                                          &
2960       kind              = "surface forcing",                               &
2961       input_id          = 1,                                               &
2962       output_file       = output_file,                                     &
2963       grid              = palm_grid,                                       &
2964       intermediate_grid = palm_intermediate                                &
2965    )
2966
2967    output_var_table(41) = init_nc_var(                                     &
2968       name              = 'rad_lw_bal_0',                                  &
2969       std_name          = "",                                              &
2970       long_name         = "longwave radiation balance at the surface",     &
2971       units             = "W/m2",                                          &
2972       kind              = "surface forcing",                               &
2973       input_id          = 1,                                               &
2974       output_file       = output_file,                                     &
2975       grid              = palm_grid,                                       &
2976       intermediate_grid = palm_intermediate                                &
2977    )
2978!
2979!------------------------------------------------------------------------------
2980! Section 2.2: Idealized large-scale forcings
2981!------------------------------------------------------------------------------
2982    output_var_table(42) = init_nc_var(                                     &
2983       name              = 'surface_forcing_surface_pressure',              &
2984       std_name          = "",                                              &
2985       long_name         = "surface pressure",                              &
2986       units             = "Pa",                                            &
2987       kind              = "time series",                                   &
2988       input_id          = 2,                                               & ! second in (T, p) I/O group
2989       output_file       = output_file,                                     &
2990       grid              = palm_grid,                                       &
2991       intermediate_grid = palm_intermediate                                &
2992    )
2993    output_var_table(42)%averaging_grid => averaged_scalar_profile
2994
2995    output_var_table(43) = init_nc_var(                                     &
2996       name              = 'ls_forcing_ug',                                 &
2997       std_name          = "",                                              &
2998       long_name         = "geostrophic wind (u component)",                &
2999       units             = "m/s",                                           &
3000       kind              = "geostrophic",                                   &
3001       input_id          = 1,                                               &
3002       output_file       = output_file,                                     &
3003       grid              = averaged_scalar_profile,                         &
3004       intermediate_grid = averaged_scalar_profile                          &
3005    )
3006
3007    output_var_table(44) = init_nc_var(                                     &
3008       name              = 'ls_forcing_vg',                                 &
3009       std_name          = "",                                              &
3010       long_name         = "geostrophic wind (v component)",                &
3011       units             = "m/s",                                           &
3012       kind              = "geostrophic",                                   &
3013       input_id          = 1,                                               &
3014       output_file       = output_file,                                     &
3015       grid              = averaged_scalar_profile,                         &
3016       intermediate_grid = averaged_scalar_profile                          &
3017    )
3018
3019    output_var_table(45) = init_nc_var(                                     &
3020       name              = 'nudging_u',                                     &
3021       std_name          = "",                                              &
3022       long_name         = "wind component in x direction",                 &
3023       units             = "m/s",                                           &
3024       kind              = "geostrophic",                                   &
3025       input_id          = 1,                                               &
3026       output_file       = output_file,                                     &
3027       grid              = averaged_scalar_profile,                         &
3028       intermediate_grid = averaged_scalar_profile                          &
3029    )
3030    output_var_table(45)%to_be_processed = ls_forcing_variables_required
3031
3032    output_var_table(46) = init_nc_var(                                     &
3033       name              = 'nudging_v',                                     &
3034       std_name          = "",                                              &
3035       long_name         = "wind component in y direction",                 &
3036       units             = "m/s",                                           &
3037       kind              = "large-scale scalar forcing",                    &
3038       input_id          = 1,                                               &
3039       output_file       = output_file,                                     &
3040       grid              = averaged_scalar_profile,                         &
3041       intermediate_grid = averaged_scalar_profile                          &
3042    )
3043    output_var_table(46)%to_be_processed = ls_forcing_variables_required
3044
3045    output_var_table(47) = init_nc_var(                                     &
3046       name              = 'ls_forcing_sub_w',                              &
3047       std_name          = "",                                              &
3048       long_name         = "subsidence velocity of w",                      &
3049       units             = "m/s",                                           &
3050       kind              = "large-scale w forcing",                         &
3051       input_id          = 1,                                               &
3052       output_file       = output_file,                                     &
3053       grid              = averaged_scalar_profile,                         &
3054       intermediate_grid = averaged_scalar_profile                          &
3055    )
3056    output_var_table(47)%to_be_processed = ls_forcing_variables_required
3057
3058    output_var_table(48) = init_nc_var(                                     &
3059       name              = 'nudging_w',                                     &
3060       std_name          = "",                                              &
3061       long_name         = "wind component in w direction",                 &
3062       units             = "m/s",                                           &
3063       kind              = "large-scale w forcing",                         &
3064       input_id          = 1,                                               &
3065       output_file       = output_file,                                     &
3066       grid              = averaged_w_profile,                              &
3067       intermediate_grid = averaged_w_profile                               &
3068    )
3069    output_var_table(48)%to_be_processed = ls_forcing_variables_required
3070
3071
3072    output_var_table(49) = init_nc_var(                                     &
3073       name              = 'ls_forcing_adv_pt',                             &
3074       std_name          = "",                                              &
3075       long_name         = "advection of potential temperature",            &
3076       units             = "K/s",                                           &
3077       kind              = "large-scale scalar forcing",                    &
3078       input_id          = 1,                                               &
3079       output_file       = output_file,                                     &
3080       grid              = averaged_scalar_profile,                         &
3081       intermediate_grid = averaged_scalar_profile                          &
3082    )
3083    output_var_table(49)%to_be_processed = ls_forcing_variables_required
3084
3085    output_var_table(50) = init_nc_var(                                     &
3086       name              = 'ls_forcing_sub_pt',                             &
3087       std_name          = "",                                              &
3088       long_name         = "subsidence velocity of potential temperature",  &
3089       units             = "K/s",                                           &
3090       kind              = "large-scale scalar forcing",                    &
3091       input_id          = 1,                                               &
3092       output_file       = output_file,                                     &
3093       grid              = averaged_scalar_profile,                         &
3094       intermediate_grid = averaged_scalar_profile                          &
3095    )
3096    output_var_table(50)%to_be_processed = ls_forcing_variables_required
3097
3098    output_var_table(51) = init_nc_var(                                     &
3099       name              = 'nudging_pt',                                    &
3100       std_name          = "",                                              &
3101       long_name         = "potential temperature",                         &
3102       units             = "K",                                             &
3103       kind              = "large-scale scalar forcing",                    &
3104       input_id          = 1,                                               &
3105       output_file       = output_file,                                     &
3106       grid              = averaged_scalar_profile,                         &
3107       intermediate_grid = averaged_scalar_profile                          &
3108    )
3109    output_var_table(51)%to_be_processed = ls_forcing_variables_required
3110
3111    output_var_table(52) = init_nc_var(                                     &
3112       name              = 'ls_forcing_adv_qv',                             &
3113       std_name          = "",                                              &
3114       long_name         = "advection of specific humidity",                &
3115       units             = "kg/kg/s",                                       &
3116       kind              = "large-scale scalar forcing",                    &
3117       input_id          = 3,                                               &
3118       output_file       = output_file,                                     &
3119       grid              = averaged_scalar_profile,                         &
3120       intermediate_grid = averaged_scalar_profile                          &
3121    )
3122    output_var_table(52)%to_be_processed = ls_forcing_variables_required
3123
3124
3125    output_var_table(53) = init_nc_var(                                     &
3126       name              = 'ls_forcing_sub_qv',                             &
3127       std_name          = "",                                              &
3128       long_name         = "subsidence velocity of specific humidity",      &
3129       units             = "kg/kg/s",                                       &
3130       kind              = "large-scale scalar forcing",                    &
3131       input_id          = 3,                                               &
3132       output_file       = output_file,                                     &
3133       grid              = averaged_scalar_profile,                         &
3134       intermediate_grid = averaged_scalar_profile                          &
3135    )
3136    output_var_table(53)%to_be_processed = ls_forcing_variables_required
3137
3138    output_var_table(54) = init_nc_var(                                     &
3139       name              = 'nudging_qv',                                    &
3140       std_name          = "",                                              &
3141       long_name         = "specific humidity",                             &
3142       units             = "kg/kg",                                         &
3143       kind              = "large-scale scalar forcing",                    &
3144       input_id          = 3,                                               &
3145       output_file       = output_file,                                     &
3146       grid              = averaged_scalar_profile,                         &
3147       intermediate_grid = averaged_scalar_profile                          &
3148    )
3149    output_var_table(54)%to_be_processed = ls_forcing_variables_required
3150
3151    output_var_table(55) = init_nc_var(                                     &
3152       name              = 'nudging_tau',                                   &
3153       std_name          = "",                                              &
3154       long_name         = "nudging relaxation time scale",                 &
3155       units             = "s",                                             &
3156       kind              = "constant scalar profile",                       &
3157       input_id          = 1,                                               &
3158       output_file       = output_file,                                     &
3159       grid              = averaged_scalar_profile,                         &
3160       intermediate_grid = averaged_scalar_profile                          &
3161    )
3162    output_var_table(55)%to_be_processed = ls_forcing_variables_required
3163
3164
3165    output_var_table(56) = init_nc_var(                                     &
3166       name              = 'internal_density_centre',                              &
3167       std_name          = "",                                              &
3168       long_name         = "",                                              &
3169       units             = "",                                              &
3170       kind              = "internal profile",                              &
3171       input_id          = 4,                                               &
3172       output_file       = output_file,                                     &
3173       grid              = averaged_scalar_profile,                         &
3174       intermediate_grid = averaged_scalar_profile                          &
3175    )
3176    output_var_table(56)%averaging_grid => averaged_scalar_profile
3177
3178
3179    output_var_table(57) = init_nc_var(                                     &
3180       name              = 'internal_density_north',                       &
3181       std_name          = "",                                              &
3182       long_name         = "",                                              &
3183       units             = "",                                              &
3184       kind              = "internal profile",                              &
3185       input_id          = 4,                                               &
3186       output_file       = output_file,                                     &
3187       grid              = north_averaged_scalar_profile,                   &
3188       intermediate_grid = north_averaged_scalar_profile                    &
3189    )
3190    output_var_table(57)%averaging_grid => north_averaged_scalar_profile
3191    output_var_table(57)%to_be_processed = .NOT. cfg%ug_defined_by_user
3192
3193
3194    output_var_table(58) = init_nc_var(                                     &
3195       name              = 'internal_density_south',                       &
3196       std_name          = "",                                              &
3197       long_name         = "",                                              &
3198       units             = "",                                              &
3199       kind              = "internal profile",                              &
3200       input_id          = 4,                                               &
3201       output_file       = output_file,                                     &
3202       grid              = south_averaged_scalar_profile,                   &
3203       intermediate_grid = south_averaged_scalar_profile                    &
3204    )
3205    output_var_table(58)%averaging_grid => south_averaged_scalar_profile
3206    output_var_table(58)%to_be_processed = .NOT. cfg%ug_defined_by_user
3207
3208
3209    output_var_table(59) = init_nc_var(                                     &
3210       name              = 'internal_density_east',                        &
3211       std_name          = "",                                              &
3212       long_name         = "",                                              &
3213       units             = "",                                              &
3214       kind              = "internal profile",                              &
3215       input_id          = 4,                                               &
3216       output_file       = output_file,                                     &
3217       grid              = east_averaged_scalar_profile,                    &
3218       intermediate_grid = east_averaged_scalar_profile                     &
3219    )
3220    output_var_table(59)%averaging_grid => east_averaged_scalar_profile
3221    output_var_table(59)%to_be_processed = .NOT. cfg%ug_defined_by_user
3222
3223
3224    output_var_table(60) = init_nc_var(                                     &
3225       name              = 'internal_density_west',                        &
3226       std_name          = "",                                              &
3227       long_name         = "",                                              &
3228       units             = "",                                              &
3229       kind              = "internal profile",                              &
3230       input_id          = 4,                                               &
3231       output_file       = output_file,                                     &
3232       grid              = west_averaged_scalar_profile,                    &
3233       intermediate_grid = west_averaged_scalar_profile                     &
3234    )
3235    output_var_table(60)%averaging_grid => west_averaged_scalar_profile
3236    output_var_table(60)%to_be_processed = .NOT. cfg%ug_defined_by_user
3237
3238    output_var_table(61) = init_nc_var(                                     &
3239       name              = 'internal_pressure_north',                       &
3240       std_name          = "",                                              &
3241       long_name         = "",                                              &
3242       units             = "",                                              &
3243       kind              = "internal profile",                              &
3244       input_id          = 2,                                               &
3245       output_file       = output_file,                                     &
3246       grid              = north_averaged_scalar_profile,                   &
3247       intermediate_grid = north_averaged_scalar_profile                    &
3248    )
3249    output_var_table(61)%averaging_grid => north_averaged_scalar_profile
3250    output_var_table(61)%to_be_processed = .NOT. cfg%ug_defined_by_user
3251
3252
3253    output_var_table(62) = init_nc_var(                                     &
3254       name              = 'internal_pressure_south',                       &
3255       std_name          = "",                                              &
3256       long_name         = "",                                              &
3257       units             = "",                                              &
3258       kind              = "internal profile",                              &
3259       input_id          = 2,                                               &
3260       output_file       = output_file,                                     &
3261       grid              = south_averaged_scalar_profile,                   &
3262       intermediate_grid = south_averaged_scalar_profile                    &
3263    )
3264    output_var_table(62)%averaging_grid => south_averaged_scalar_profile
3265    output_var_table(62)%to_be_processed = .NOT. cfg%ug_defined_by_user
3266
3267
3268    output_var_table(63) = init_nc_var(                                     &
3269       name              = 'internal_pressure_east',                        &
3270       std_name          = "",                                              &
3271       long_name         = "",                                              &
3272       units             = "",                                              &
3273       kind              = "internal profile",                              &
3274       input_id          = 2,                                               &
3275       output_file       = output_file,                                     &
3276       grid              = east_averaged_scalar_profile,                    &
3277       intermediate_grid = east_averaged_scalar_profile                     &
3278    )
3279    output_var_table(63)%averaging_grid => east_averaged_scalar_profile
3280    output_var_table(63)%to_be_processed = .NOT. cfg%ug_defined_by_user
3281
3282
3283    output_var_table(64) = init_nc_var(                                     &
3284       name              = 'internal_pressure_west',                        &
3285       std_name          = "",                                              &
3286       long_name         = "",                                              &
3287       units             = "",                                              &
3288       kind              = "internal profile",                              &
3289       input_id          = 2,                                               &
3290       output_file       = output_file,                                     &
3291       grid              = west_averaged_scalar_profile,                    &
3292       intermediate_grid = west_averaged_scalar_profile                     &
3293    )
3294    output_var_table(64)%averaging_grid => west_averaged_scalar_profile
3295    output_var_table(64)%to_be_processed = .NOT. cfg%ug_defined_by_user
3296
3297!
3298!-- Attributes shared among all variables
3299    output_var_table(:)%source = nc_source_text
3300
3301
3302 END SUBROUTINE setup_variable_tables
3303
3304
3305!------------------------------------------------------------------------------!
3306! Description:
3307! ------------
3308!> Initializes an nc_var varible with the given parameters. The 'kind'
3309!> parameter is used to infer the correct netCDF IDs and the level of detail,
3310!> 'lod', as defined by the PALM-4U input data standard.
3311!------------------------------------------------------------------------------!
3312 FUNCTION init_nc_var(name, std_name, long_name, units, kind, input_id,     &
3313                      grid, intermediate_grid, output_file, is_profile)     &
3314    RESULT(var)
3315
3316    CHARACTER(LEN=*), INTENT(IN)      ::  name, std_name, long_name, units, kind
3317    INTEGER, INTENT(IN)               ::  input_id
3318    TYPE(grid_definition), INTENT(IN), TARGET ::  grid, intermediate_grid
3319    TYPE(nc_file), INTENT(IN)         ::  output_file
3320    LOGICAL, INTENT(IN), OPTIONAL     ::  is_profile
3321
3322    CHARACTER(LEN=LNAME)              ::  out_var_kind
3323    TYPE(nc_var)                      ::  var
3324
3325    out_var_kind = TRIM(kind)
3326
3327    IF (PRESENT(is_profile))  THEN
3328       IF (is_profile)  out_var_kind = TRIM(kind) // ' profile'
3329    ENDIF
3330
3331    var%name              = name
3332    var%standard_name     = std_name
3333    var%long_name         = long_name
3334    var%units             = units
3335    var%kind              = TRIM(out_var_kind)
3336    var%input_id          = input_id
3337    var%nt                = SIZE (output_file%time)
3338    var%grid              => grid
3339    var%intermediate_grid => intermediate_grid
3340
3341    SELECT CASE( TRIM(out_var_kind) )
3342
3343!
3344!--    TODO: Using global module variables 'init_variables_required' and
3345!--    TODO: 'boundary_variables_required'. Encapsulate in settings type
3346!--    TODO: and pass into init_nc_var.
3347       CASE( 'init soil' )
3348          var%nt              = 1
3349          var%lod             = 2
3350          var%ndim            = 3
3351          var%dimids(1:3)     = output_file%dimids_soil
3352          var%dimvarids(1:3)  = output_file%dimvarids_soil
3353          var%to_be_processed = init_variables_required
3354          var%is_internal     = .FALSE.
3355          var%task            = "interpolate_2d"
3356
3357       CASE( 'init scalar' )
3358          var%nt              = 1
3359          var%lod             = 2
3360          var%ndim            = 3
3361          var%dimids(1:3)     = output_file%dimids_scl
3362          var%dimvarids(1:3)  = output_file%dimvarids_scl
3363          var%to_be_processed = init_variables_required
3364          var%is_internal     = .FALSE.
3365          var%task            = "interpolate_3d"
3366
3367       CASE( 'init u' )
3368          var%nt              = 1
3369          var%lod             = 2
3370          var%ndim            = 3
3371          var%dimids(1)       = output_file%dimids_vel(1)
3372          var%dimids(2)       = output_file%dimids_scl(2)
3373          var%dimids(3)       = output_file%dimids_scl(3)
3374          var%dimvarids(1)    = output_file%dimvarids_vel(1)
3375          var%dimvarids(2)    = output_file%dimvarids_scl(2)
3376          var%dimvarids(3)    = output_file%dimvarids_scl(3)
3377          var%to_be_processed = init_variables_required
3378          var%is_internal     = .FALSE.
3379          var%task            = "interpolate_3d"
3380
3381       CASE( 'init v' )
3382          var%nt              = 1
3383          var%lod             = 2
3384          var%ndim            = 3
3385          var%dimids(1)       = output_file%dimids_scl(1)
3386          var%dimids(2)       = output_file%dimids_vel(2)
3387          var%dimids(3)       = output_file%dimids_scl(3)
3388          var%dimvarids(1)    = output_file%dimvarids_scl(1)
3389          var%dimvarids(2)    = output_file%dimvarids_vel(2)
3390          var%dimvarids(3)    = output_file%dimvarids_scl(3)
3391          var%to_be_processed = init_variables_required
3392          var%is_internal     = .FALSE.
3393          var%task            = "interpolate_3d"
3394
3395       CASE( 'init w' )
3396          var%nt              = 1
3397          var%lod             = 2
3398          var%ndim            = 3
3399          var%dimids(1)       = output_file%dimids_scl(1)
3400          var%dimids(2)       = output_file%dimids_scl(2)
3401          var%dimids(3)       = output_file%dimids_vel(3)
3402          var%dimvarids(1)    = output_file%dimvarids_scl(1)
3403          var%dimvarids(2)    = output_file%dimvarids_scl(2)
3404          var%dimvarids(3)    = output_file%dimvarids_vel(3)
3405          var%to_be_processed = init_variables_required
3406          var%is_internal     = .FALSE.
3407          var%task            = "interpolate_3d"
3408
3409       CASE( 'init scalar profile', 'init u profile', 'init v profile')
3410          var%nt              = 1
3411          var%lod             = 1
3412          var%ndim            = 1
3413          var%dimids(1)       = output_file%dimids_scl(3)    !z
3414          var%dimvarids(1)    = output_file%dimvarids_scl(3) !z
3415          var%to_be_processed = init_variables_required
3416          var%is_internal     = .FALSE.
3417          var%task            = "average profile"
3418
3419       CASE( 'init w profile')
3420          var%nt              = 1
3421          var%lod             = 1
3422          var%ndim            = 1
3423          var%dimids(1)       = output_file%dimids_vel(3)    !z
3424          var%dimvarids(1)    = output_file%dimvarids_vel(3) !z
3425          var%to_be_processed = init_variables_required
3426          var%is_internal     = .FALSE.
3427          var%task            = "average profile"
3428
3429       CASE( 'surface forcing' )
3430          var%lod             = -1
3431          var%ndim            = 3
3432          var%dimids(3)       = output_file%dimid_time
3433          var%dimids(1:2)     = output_file%dimids_soil(1:2)
3434          var%dimvarids(3)    = output_file%dimvarid_time
3435          var%dimvarids(1:2)  = output_file%dimvarids_soil(1:2)
3436          var%to_be_processed = surface_forcing_required
3437          var%is_internal     = .FALSE.
3438          var%task            = "interpolate_2d"
3439
3440       CASE( 'left scalar', 'right scalar')
3441          var%lod             = -1
3442          var%ndim            = 3
3443          var%dimids(3)       = output_file%dimid_time
3444          var%dimids(1)       = output_file%dimids_scl(2)
3445          var%dimids(2)       = output_file%dimids_scl(3)
3446          var%dimvarids(3)    = output_file%dimvarid_time
3447          var%dimvarids(1)    = output_file%dimvarids_scl(2)
3448          var%dimvarids(2)    = output_file%dimvarids_scl(3)
3449          var%to_be_processed = boundary_variables_required
3450          var%is_internal     = .FALSE.
3451          var%task            = "interpolate_3d"
3452
3453       CASE( 'north scalar', 'south scalar')
3454          var%lod             = -1
3455          var%ndim            = 3
3456          var%dimids(3)       = output_file%dimid_time
3457          var%dimids(1)       = output_file%dimids_scl(1)
3458          var%dimids(2)       = output_file%dimids_scl(3)
3459          var%dimvarids(3)    = output_file%dimvarid_time
3460          var%dimvarids(1)    = output_file%dimvarids_scl(1)
3461          var%dimvarids(2)    = output_file%dimvarids_scl(3)
3462          var%to_be_processed = boundary_variables_required
3463          var%is_internal     = .FALSE.
3464          var%task            = "interpolate_3d"
3465
3466       CASE( 'top scalar', 'top w' )
3467          var%lod             = -1
3468          var%ndim            = 3
3469          var%dimids(3)       = output_file%dimid_time
3470          var%dimids(1)       = output_file%dimids_scl(1)
3471          var%dimids(2)       = output_file%dimids_scl(2)
3472          var%dimvarids(3)    = output_file%dimvarid_time
3473          var%dimvarids(1)    = output_file%dimvarids_scl(1)
3474          var%dimvarids(2)    = output_file%dimvarids_scl(2)
3475          var%to_be_processed = boundary_variables_required
3476          var%is_internal     = .FALSE.
3477          var%task            = "interpolate_3d"
3478
3479       CASE( 'left u', 'right u' )
3480          var%lod             = -1
3481          var%ndim            = 3
3482          var%dimids(3)       = output_file%dimid_time
3483          var%dimids(1)       = output_file%dimids_scl(2)
3484          var%dimids(2)       = output_file%dimids_scl(3)
3485          var%dimvarids(3)    = output_file%dimvarid_time
3486          var%dimvarids(1)    = output_file%dimvarids_scl(2)
3487          var%dimvarids(2)    = output_file%dimvarids_scl(3)
3488          var%to_be_processed = boundary_variables_required
3489          var%is_internal     = .FALSE.
3490          var%task            = "interpolate_3d"
3491
3492       CASE( 'north u', 'south u' )
3493          var%lod             = -1
3494          var%ndim            = 3
3495          var%dimids(3)       = output_file%dimid_time    !t
3496          var%dimids(1)       = output_file%dimids_vel(1) !x
3497          var%dimids(2)       = output_file%dimids_scl(3) !z
3498          var%dimvarids(3)    = output_file%dimvarid_time
3499          var%dimvarids(1)    = output_file%dimvarids_vel(1)
3500          var%dimvarids(2)    = output_file%dimvarids_scl(3)
3501          var%to_be_processed = boundary_variables_required
3502          var%is_internal     = .FALSE.
3503          var%task            = "interpolate_3d"
3504
3505       CASE( 'top u' )
3506          var%lod             = -1
3507          var%ndim            = 3
3508          var%dimids(3)       = output_file%dimid_time    !t
3509          var%dimids(1)       = output_file%dimids_vel(1) !x
3510          var%dimids(2)       = output_file%dimids_scl(2) !z
3511          var%dimvarids(3)    = output_file%dimvarid_time
3512          var%dimvarids(1)    = output_file%dimvarids_vel(1)
3513          var%dimvarids(2)    = output_file%dimvarids_scl(2)
3514          var%to_be_processed = boundary_variables_required
3515          var%is_internal     = .FALSE.
3516          var%task            = "interpolate_3d"
3517
3518       CASE( 'left v', 'right v' )
3519          var%lod             = -1
3520          var%ndim            = 3
3521          var%dimids(3)       = output_file%dimid_time
3522          var%dimids(1)       = output_file%dimids_vel(2)
3523          var%dimids(2)       = output_file%dimids_scl(3)
3524          var%dimvarids(3)    = output_file%dimvarid_time
3525          var%dimvarids(1)    = output_file%dimvarids_vel(2)
3526          var%dimvarids(2)    = output_file%dimvarids_scl(3)
3527          var%to_be_processed = boundary_variables_required
3528          var%is_internal     = .FALSE.
3529          var%task            = "interpolate_3d"
3530
3531       CASE( 'north v', 'south v' )
3532          var%lod             = -1
3533          var%ndim            = 3
3534          var%dimids(3)       = output_file%dimid_time    !t
3535          var%dimids(1)       = output_file%dimids_scl(1) !x
3536          var%dimids(2)       = output_file%dimids_scl(3) !z
3537          var%dimvarids(3)    = output_file%dimvarid_time
3538          var%dimvarids(1)    = output_file%dimvarids_scl(1)
3539          var%dimvarids(2)    = output_file%dimvarids_scl(3)
3540          var%to_be_processed = boundary_variables_required
3541          var%is_internal     = .FALSE.
3542          var%task            = "interpolate_3d"
3543
3544       CASE( 'top v' )
3545          var%lod             = -1
3546          var%ndim            = 3
3547          var%dimids(3)       = output_file%dimid_time    !t
3548          var%dimids(1)       = output_file%dimids_scl(1) !x
3549          var%dimids(2)       = output_file%dimids_vel(2) !z
3550          var%dimvarids(3)    = output_file%dimvarid_time
3551          var%dimvarids(1)    = output_file%dimvarids_scl(1)
3552          var%dimvarids(2)    = output_file%dimvarids_vel(2)
3553          var%to_be_processed = boundary_variables_required
3554          var%is_internal     = .FALSE.
3555          var%task            = "interpolate_3d"
3556
3557       CASE( 'left w', 'right w')
3558          var%lod             = -1
3559          var%ndim            = 3
3560          var%dimids(3)       = output_file%dimid_time
3561          var%dimids(1)       = output_file%dimids_scl(2)
3562          var%dimids(2)       = output_file%dimids_vel(3)
3563          var%dimvarids(3)    = output_file%dimvarid_time
3564          var%dimvarids(1)    = output_file%dimvarids_scl(2)
3565          var%dimvarids(2)    = output_file%dimvarids_vel(3)
3566          var%to_be_processed = boundary_variables_required
3567          var%is_internal     = .FALSE.
3568          var%task            = "interpolate_3d"
3569
3570       CASE( 'north w', 'south w' )
3571          var%lod             = -1
3572          var%ndim            = 3
3573          var%dimids(3)       = output_file%dimid_time    !t
3574          var%dimids(1)       = output_file%dimids_scl(1) !x
3575          var%dimids(2)       = output_file%dimids_vel(3) !z
3576          var%dimvarids(3)    = output_file%dimvarid_time
3577          var%dimvarids(1)    = output_file%dimvarids_scl(1)
3578          var%dimvarids(2)    = output_file%dimvarids_vel(3)
3579          var%to_be_processed = boundary_variables_required
3580          var%is_internal     = .FALSE.
3581          var%task            = "interpolate_3d"
3582
3583       CASE( 'time series' )
3584          var%lod             = 0
3585          var%ndim            = 1
3586          var%dimids(1)       = output_file%dimid_time    !t
3587          var%dimvarids(1)    = output_file%dimvarid_time
3588          var%to_be_processed = .TRUE.
3589          var%is_internal     = .FALSE.
3590          var%task            = "average profile"
3591
3592       CASE( 'constant scalar profile' )
3593          var%lod             = -1
3594          var%ndim            = 2
3595          var%dimids(2)       = output_file%dimid_time    !t
3596          var%dimids(1)       = output_file%dimids_scl(3) !z
3597          var%dimvarids(2)    = output_file%dimvarid_time
3598          var%dimvarids(1)    = output_file%dimvarids_scl(3)
3599          var%to_be_processed = .TRUE.
3600          var%is_internal     = .FALSE.
3601          var%task            = "set profile"
3602
3603       CASE( 'large-scale scalar forcing' )
3604          var%lod             = -1
3605          var%ndim            = 2
3606          var%dimids(2)       = output_file%dimid_time    !t
3607          var%dimids(1)       = output_file%dimids_scl(3) !z
3608          var%dimvarids(2)    = output_file%dimvarid_time
3609          var%dimvarids(1)    = output_file%dimvarids_scl(3)
3610          var%to_be_processed = ls_forcing_variables_required
3611          var%is_internal     = .FALSE.
3612          var%task            = "average large-scale profile"
3613
3614       CASE( 'geostrophic' )
3615          var%lod             = -1
3616          var%ndim            = 2
3617          var%dimids(2)       = output_file%dimid_time    !t
3618          var%dimids(1)       = output_file%dimids_scl(3) !z
3619          var%dimvarids(2)    = output_file%dimvarid_time
3620          var%dimvarids(1)    = output_file%dimvarids_scl(3)
3621          var%to_be_processed = .TRUE.
3622          var%is_internal     = .FALSE.
3623          var%task            = "geostrophic winds"
3624
3625       CASE( 'large-scale w forcing' )
3626          var%lod             = -1
3627          var%ndim            = 2
3628          var%dimids(2)       = output_file%dimid_time    !t
3629          var%dimids(1)       = output_file%dimids_vel(3) !z
3630          var%dimvarids(2)    = output_file%dimvarid_time
3631          var%dimvarids(1)    = output_file%dimvarids_vel(3)
3632          var%to_be_processed = ls_forcing_variables_required
3633          var%is_internal     = .FALSE.
3634          var%task            = "average large-scale profile"
3635
3636       CASE( 'internal profile' )
3637          var%lod             = -1
3638          var%ndim            = 2
3639          var%dimids(2)       = output_file%dimid_time    !t
3640          var%dimids(1)       = output_file%dimids_scl(3) !z
3641          var%dimvarids(2)    = output_file%dimvarid_time
3642          var%dimvarids(1)    = output_file%dimvarids_scl(3)
3643          var%to_be_processed = .TRUE.
3644          var%is_internal     = .TRUE.
3645          var%task            = "internal profile"
3646
3647       CASE DEFAULT
3648           message = "Variable kind '" // TRIM(kind) // "' not recognized."
3649           CALL inifor_abort ('init_nc_var', message)
3650
3651    END SELECT
3652
3653 END FUNCTION init_nc_var
3654
3655
3656 SUBROUTINE fini_variables()
3657
3658    CALL report('fini_variables', 'Deallocating variable table', cfg%debug)
3659    DEALLOCATE( input_var_table )
3660
3661 END SUBROUTINE fini_variables
3662
3663
3664 SUBROUTINE fini_io_groups()
3665
3666    CALL report('fini_io_groups', 'Deallocating IO groups', cfg%debug)
3667    DEALLOCATE( io_group_list )
3668
3669 END SUBROUTINE fini_io_groups
3670
3671
3672 SUBROUTINE fini_file_lists()
3673   
3674    CALL report('fini_file_lists', 'Deallocating file lists', cfg%debug)
3675    DEALLOCATE( flow_files, soil_files, radiation_files, soil_moisture_files )
3676
3677 END SUBROUTINE fini_file_lists
3678
3679
3680!------------------------------------------------------------------------------!
3681! Description:
3682! ------------
3683!> Carries out any physical conversion of the quantities in the given input
3684!> buffer needed to obtain the quantity required by PALM-4U. For instance,
3685!> velocities are rotated to the PALM-4U coordinate system and the potential
3686!> temperature is computed from the absolute temperature and pressure.
3687!>
3688!> Note, that the preprocessing does not include any grid change. The result
3689!> array will match a COSMO-DE scalar array.
3690!------------------------------------------------------------------------------!
3691 SUBROUTINE preprocess(group, input_buffer, cosmo_grid, iter)
3692
3693    TYPE(io_group), INTENT(INOUT), TARGET       ::  group
3694    TYPE(container), INTENT(INOUT), ALLOCATABLE ::  input_buffer(:)
3695    TYPE(grid_definition), INTENT(IN)           ::  cosmo_grid
3696    INTEGER, INTENT(IN)                         ::  iter
3697   
3698    REAL(wp), ALLOCATABLE                       ::  basic_state_pressure(:)
3699    TYPE(container), ALLOCATABLE                ::  preprocess_buffer(:)
3700    INTEGER                                     ::  hour, dt
3701    INTEGER                                     ::  i, j, k
3702    INTEGER                                     ::  nx, ny, nz
3703   
3704    input_buffer(:)%is_preprocessed = .FALSE.
3705     
3706    SELECT CASE( group%kind )
3707       
3708       CASE( 'velocities' )
3709!
3710!--       Allocate a compute buffer with the same number of arrays as the input
3711          ALLOCATE( preprocess_buffer( SIZE(input_buffer) ) )
3712
3713!
3714!--       Allocate u and v arrays with scalar dimensions
3715          nx = SIZE(input_buffer(1)%array, 1)
3716          ny = SIZE(input_buffer(1)%array, 2)
3717          nz = SIZE(input_buffer(1)%array, 3)
3718          ALLOCATE( preprocess_buffer(1)%array(nx, ny, nz) ) ! u buffer
3719          ALLOCATE( preprocess_buffer(2)%array(nx, ny, nz) ) ! v buffer
3720
3721          CALL log_runtime('time', 'alloc')
3722
3723!
3724!--       interpolate U and V to centres
3725          CALL centre_velocities( u_face = input_buffer(1)%array,            &
3726                                  v_face = input_buffer(2)%array,            &
3727                                  u_centre = preprocess_buffer(1)%array,     &
3728                                  v_centre = preprocess_buffer(2)%array )
3729         
3730          cfg%rotation_method = 'rotated-pole'
3731          SELECT CASE(cfg%rotation_method)
3732
3733             CASE('rotated-pole')
3734!           
3735!--             rotate U and V to PALM-4U orientation and overwrite U and V with
3736!--             rotated velocities
3737                DO  k = 1, nz
3738                DO  j = 1, ny
3739                DO  i = 1, nx
3740                   CALL uv2uvrot( urot = preprocess_buffer(1)%array(i,j,k),     &
3741                                  vrot = preprocess_buffer(2)%array(i,j,k),     &
3742                                  rlat = cosmo_grid%lat(j-1),                   &
3743                                  rlon = cosmo_grid%lon(i-1),                   &
3744                                  pollat = phi_cn,                                &
3745                                  pollon = lambda_cn,                             &
3746                                  u = input_buffer(1)%array(i,j,k),             &
3747                                  v = input_buffer(2)%array(i,j,k) )
3748                ENDDO
3749                ENDDO
3750                ENDDO
3751             
3752             CASE DEFAULT
3753                message = "Rotation method '" // TRIM(cfg%rotation_method) //   &
3754                   "' not recognized."
3755                CALL inifor_abort('preprocess', message)
3756
3757          END SELECT
3758
3759          input_buffer(1)%array(1,:,:) = 0.0_wp
3760          input_buffer(2)%array(1,:,:) = 0.0_wp
3761          input_buffer(1)%array(:,1,:) = 0.0_wp
3762          input_buffer(2)%array(:,1,:) = 0.0_wp
3763
3764          input_buffer(1:2)%is_preprocessed = .TRUE.
3765          CALL log_runtime('time', 'comp')
3766
3767          DEALLOCATE( preprocess_buffer )
3768          CALL log_runtime('time', 'alloc')
3769
3770          message = "Input buffers for group '" // TRIM(group%kind) // "'"//&
3771             " preprocessed sucessfully."
3772          CALL report('preprocess', message)
3773       
3774       CASE( 'thermodynamics' ) ! T, P, QV
3775          nx = SIZE(input_buffer(1)%array, 1)
3776          ny = SIZE(input_buffer(1)%array, 2)
3777          nz = SIZE(input_buffer(1)%array, 3)
3778
3779!
3780!--       Compute absolute pressure if presure perturbation has been read in.
3781          IF ( TRIM(group%in_var_list(2)%name) == 'PP' )  THEN
3782             message = "Absolute pressure, P, not available, " //              &
3783                       "computing from pressure preturbation PP."
3784             CALL report('preprocess', message)
3785
3786             ALLOCATE( basic_state_pressure(1:nz) )
3787             CALL log_runtime('time', 'alloc')
3788
3789             DO  j = 1, ny
3790             DO  i = 1, nx
3791
3792                CALL get_basic_state(cosmo_grid%hfl(i,j,:), BETA, P_SL, T_SL,&
3793                                     RD, G, basic_state_pressure)
3794
3795!
3796!--             Overwrite pressure perturbation with absolute pressure. HECTO
3797!--             converts pressure perturbation from hPa to Pa.
3798                input_buffer (2)%array(i,j,:) =                              &
3799                   HECTO * input_buffer (2)%array(i,j,:) +                   &
3800                   basic_state_pressure(:)
3801
3802             ENDDO
3803             ENDDO
3804             CALL log_runtime('time', 'comp')
3805
3806             DEALLOCATE( basic_state_pressure )
3807             CALL log_runtime('time', 'alloc')
3808
3809             group%in_var_list(2)%name = 'P'
3810
3811          ENDIF
3812!
3813!--       mark pressure as preprocessed
3814          input_buffer(2)%is_preprocessed = .TRUE.
3815
3816!
3817!--       Copy temperature to the last input buffer array
3818          ALLOCATE(                                                            &
3819              input_buffer( group%n_output_quantities )%array (nx, ny, nz) &
3820          )
3821          CALL log_runtime('time', 'alloc')
3822          input_buffer(group%n_output_quantities)%array(:,:,:) =           &
3823              input_buffer(1)%array(:,:,:)
3824
3825!
3826!--       Convert absolute in place to potential temperature
3827          CALL potential_temperature(                                          &
3828             t = input_buffer(1)%array(:,:,:),                               &
3829             p = input_buffer(2)%array(:,:,:),                               &
3830             p_ref = P_REF,                                                    &
3831             r = RD_PALM,                                                      &
3832             cp = CP_PALM                                                      &
3833          )
3834
3835!
3836!--       mark potential temperature as preprocessed
3837          input_buffer(1)%is_preprocessed = .TRUE.
3838
3839!
3840!--       Convert temperature copy to density
3841          CALL moist_density(                                                  &
3842             t_rho = input_buffer(group%n_output_quantities)%array(:,:,:), &
3843             p = input_buffer(2)%array(:,:,:),                               &
3844             qv = input_buffer(3)%array(:,:,:),                              &
3845             rd = RD,                                                          &
3846             rv = RV                                                           &
3847          )
3848
3849!
3850!--       mark qv as preprocessed
3851          input_buffer(3)%is_preprocessed = .TRUE.
3852
3853!
3854!--       mark density as preprocessed
3855          input_buffer(group%n_output_quantities)%is_preprocessed = .TRUE.
3856
3857
3858          message = "Input buffers for group '" // TRIM(group%kind) // "'"//&
3859             " preprocessed sucessfully."
3860          CALL report('preprocess', message)
3861       
3862       CASE( 'scalar' ) ! S or W
3863          input_buffer(:)%is_preprocessed = .TRUE.
3864
3865       CASE( 'soil-temperature' ) !
3866         
3867          CALL fill_water_cells(soiltyp, input_buffer(1)%array, &
3868                                SIZE(input_buffer(1)%array, 3), &
3869                                FILL_ITERATIONS)
3870          input_buffer(:)%is_preprocessed = .TRUE.
3871
3872       CASE( 'soil-water' ) !
3873
3874          CALL fill_water_cells(soiltyp, input_buffer(1)%array, &
3875                                SIZE(input_buffer(1)%array, 3), &
3876                                FILL_ITERATIONS)
3877
3878          nx = SIZE(input_buffer(1)%array, 1)
3879          ny = SIZE(input_buffer(1)%array, 2)
3880          nz = SIZE(input_buffer(1)%array, 3)
3881
3882          DO  k = 1, nz
3883          DO  j = 1, ny
3884          DO  i = 1, nx
3885             input_buffer(1)%array(i,j,k) =                                  &
3886                 input_buffer(1)%array(i,j,k) * d_depth_rho_inv(k)
3887          ENDDO
3888          ENDDO
3889          ENDDO
3890
3891          message = "Converted soil water from [kg/m^2] to [m^3/m^3]"
3892          CALL report('preprocess', message)
3893
3894          input_buffer(:)%is_preprocessed = .TRUE.
3895
3896       CASE( 'surface' ) !
3897          input_buffer(:)%is_preprocessed = .TRUE.
3898
3899       CASE( 'accumulated' ) !
3900          message = "De-accumulating '" // TRIM(group%in_var_list(1)%name) //&
3901                    "' in iteration " // TRIM(str(iter)) 
3902          CALL report('preprocess', message)
3903
3904          hour = iter - 1
3905          dt = MODULO(hour, 3) + 1 ! averaging period
3906          SELECT CASE(dt)
3907
3908!
3909!--          input has been accumulated over one hour. Leave as is
3910!--          input_buffer(1)%array(:,:,:) carrries one-hour integral
3911             CASE(1)
3912             
3913!           
3914!--          input has been accumulated over two hours. Subtract previous step
3915!--          input_buffer(1)%array(:,:,:) carrries one-hour integral
3916!--          input_buffer(2)%array(:,:,:) carrries two-hour integral
3917             CASE(2)
3918                CALL deaverage(                                                   &
3919                         avg_1 = input_buffer(1)%array(:,:,:), t1 = 1.0_wp,     &
3920                         avg_2 = input_buffer(2)%array(:,:,:), t2 = 1.0_wp,     &
3921                         avg_3 = input_buffer(1)%array(:,:,:), t3 = 1.0_wp )
3922!           
3923!--             input_buffer(1)%array(:,:,:) carrries one-hour integral of second hour
3924             
3925!           
3926!--          input has been accumulated over three hours. Subtract previous step
3927!--          input_buffer(1)%array(:,:,:) carrries three-hour integral
3928!--          input_buffer(2)%array(:,:,:) still carrries two-hour integral
3929             CASE(3)
3930                CALL deaverage(                                                   &
3931                        avg_1 = input_buffer(2)%array(:,:,:), t1 = 1.0_wp,      &
3932                        avg_2 = input_buffer(1)%array(:,:,:), t2 = 1.0_wp,      &
3933                        avg_3 = input_buffer(1)%array(:,:,:), t3 = 1.0_wp )
3934!           
3935!--             input_buffer(1)%array(:,:,:) carrries one-hour integral of third hourA
3936             
3937             CASE DEFAULT
3938                message = "Invalid averaging period '" // TRIM(str(dt)) // " hours"
3939                CALL inifor_abort('preprocess', message)
3940
3941          END SELECT
3942          input_buffer(:)%is_preprocessed = .TRUE.
3943
3944       CASE( 'running average' ) !
3945          message = "De-averaging '" // TRIM(group%in_var_list(1)%name) //   &
3946                    "' in iteration " // TRIM(str(iter)) 
3947          CALL report('preprocess', message)
3948
3949          hour = iter - 1
3950!
3951!--       averaging period
3952          dt = MODULO(hour, 3) + 1
3953          SELECT CASE(dt)
3954!
3955!--          input has been accumulated over one hour. Leave as is
3956!--          input_buffer(1)%array(:,:,:) carrries one-hour integral
3957             CASE(1)
3958             
3959!           
3960!--          input has been accumulated over two hours. Subtract previous step
3961!--          input_buffer(1)%array(:,:,:) carrries one-hour integral
3962!--          input_buffer(2)%array(:,:,:) carrries two-hour integral
3963             CASE(2)
3964                CALL deaverage( input_buffer(1)%array(:,:,:), 1.0_wp,           &
3965                                input_buffer(2)%array(:,:,:), 2.0_wp,           &
3966                                input_buffer(1)%array(:,:,:), 1.0_wp)
3967!           
3968!--             input_buffer(1)%array(:,:,:) carrries one-hour integral of second hour
3969             
3970!           
3971!--          input has been accumulated over three hours. Subtract previous step
3972!--          input_buffer(1)%array(:,:,:) carrries three-hour integral
3973!--          input_buffer(2)%array(:,:,:) still carrries two-hour integral
3974             CASE(3)
3975                CALL deaverage( input_buffer(2)%array(:,:,:), 2.0_wp,           &
3976                                input_buffer(1)%array(:,:,:), 3.0_wp,           &
3977                                input_buffer(1)%array(:,:,:), 1.0_wp)
3978!           
3979!--             input_buffer(1)%array(:,:,:) carrries one-hour integral of third hourA
3980             
3981             CASE DEFAULT
3982                message = "Invalid averaging period '" // TRIM(str(dt)) // " hours"
3983                CALL inifor_abort('preprocess', message)
3984
3985          END SELECT
3986          input_buffer(:)%is_preprocessed = .TRUE.
3987
3988       CASE DEFAULT
3989          message = "IO group kind '" // TRIM(group%kind) // "' is not supported."
3990          CALL inifor_abort('prerpocess', message)
3991
3992    END SELECT
3993    CALL log_runtime('time', 'comp')
3994
3995 END SUBROUTINE preprocess
3996
3997
3998!------------------------------------------------------------------------------!
3999! Description:
4000! ------------
4001!> Computes average soil values in COSMO-DE water cells from neighbouring
4002!> non-water cells. This is done as a preprocessing step for the COSMO-DE
4003!> soil input arrays, which contain unphysical values for water cells.
4004!>
4005!> This routine computes the average of up to all nine neighbouring cells
4006!> or keeps the original value, if not at least one non-water neightbour
4007!> is available.
4008!>
4009!> By repeatedly applying this step, soil data can be extrapolated into
4010!> 'water' regions occupying multiple cells, one cell per iteration.
4011!>
4012!> Input parameters:
4013!> -----------------
4014!> soiltyp : 2d map of COSMO-DE soil types
4015!> nz : number of layers in the COSMO-DE soil
4016!> niter : number iterations
4017!>
4018!> Output parameters:
4019!> ------------------
4020!> array : the soil array (i.e. water content or temperature)
4021!------------------------------------------------------------------------------!
4022 SUBROUTINE fill_water_cells(soiltyp, array, nz, niter)
4023    INTEGER(iwp), DIMENSION(:,:,:), INTENT(IN) :: soiltyp
4024    REAL(wp), DIMENSION(:,:,:), INTENT(INOUT)  :: array
4025    INTEGER, INTENT(IN)                        :: nz, niter
4026
4027    REAL(wp), DIMENSION(nz)                    :: column
4028    INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE  :: old_soiltyp, new_soiltyp
4029    INTEGER                                    :: l, i, j, nx, ny, n_cells, ii, jj, iter
4030    INTEGER, DIMENSION(8)                      :: di, dj
4031
4032    nx = SIZE(array, 1)
4033    ny = SIZE(array, 2)
4034    di = (/ -1, -1, -1, 0,  0,  1, 1, 1 /)
4035    dj = (/ -1,  0,  1, -1, 1, -1, 0, 1 /)
4036
4037    ALLOCATE(old_soiltyp(SIZE(soiltyp,1), &
4038                         SIZE(soiltyp,2) ))
4039
4040    ALLOCATE(new_soiltyp(SIZE(soiltyp,1), &
4041                         SIZE(soiltyp,2) ))
4042
4043    old_soiltyp(:,:) = soiltyp(:,:,1)
4044    new_soiltyp(:,:) = soiltyp(:,:,1)
4045
4046    DO  iter = 1, niter
4047
4048       DO  j = 1, ny
4049       DO  i = 1, nx
4050       
4051          IF (old_soiltyp(i,j) == WATER_ID)  THEN
4052
4053             n_cells = 0
4054             column(:) = 0.0_wp
4055             DO  l = 1, SIZE(di)
4056
4057                ii = MIN(nx, MAX(1, i + di(l)))
4058                jj = MIN(ny, MAX(1, j + dj(l)))
4059
4060                IF (old_soiltyp(ii,jj) .NE. WATER_ID)  THEN
4061                   n_cells = n_cells + 1
4062                   column(:) = column(:) + array(ii,jj,:)
4063                ENDIF
4064
4065             ENDDO
4066
4067!
4068!--          Overwrite if at least one non-water neighbour cell is available
4069             IF (n_cells > 0)  THEN
4070                array(i,j,:) = column(:) / n_cells
4071                new_soiltyp(i,j) = 0
4072             ENDIF
4073
4074          ENDIF
4075
4076       ENDDO
4077       ENDDO
4078
4079       old_soiltyp(:,:) = new_soiltyp(:,:)
4080
4081    ENDDO
4082
4083    DEALLOCATE(old_soiltyp, new_soiltyp)
4084
4085 END SUBROUTINE fill_water_cells
4086
4087
4088 END MODULE inifor_grid
4089#endif
Note: See TracBrowser for help on using the repository browser.