Home

Context Navigation

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

Last change on this file since 4698 was 4675, checked in by eckhard, 4 years ago
Support for homogeneous (domain-averaged) boundary conditions and soil profile initialization
Property svn:keywords set to `Id`
File size: 238.9 KB

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

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |