Home

Context Navigation

inifor_grid.f90 @ 4866

Last change on this file since 4866 was 4843, checked in by raasch, 4 years ago
local namelist parameter added to switch off the module although the respective module namelist appears in the namelist file, further copyright updates
Property svn:keywords set to `Id`
File size: 239.4 KB

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

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

Context Navigation

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

Download in other formats: