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

Last change on this file was 4843, checked in by raasch, 4 years ago

local namelist parameter added to switch off the module although the respective module namelist appears in the namelist file, further copyright updates

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