Home

Context Navigation

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

Last change on this file since 4538 was 4538, checked in by eckhard, 4 years ago
inifor: Fix issue where --elevation/-z option was ignored, make it mandatory
Property svn:keywords set to `Id`
File size: 203.4 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |