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

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

inifor: Support for COSMO cloud water and precipitation

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