Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |