Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |