Home

Context Navigation

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

Last change on this file since 4659 was 4659, checked in by eckhard, 4 years ago
inifor: Support for COSMO cloud water and precipitation
Property svn:keywords set to `Id`
File size: 231.8 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |