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

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

Support for homogeneous (domain-averaged) boundary conditions and soil profile initialization

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