Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |