Home

Context Navigation

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

Last change on this file since 4523 was 4523, checked in by eckhard, 4 years ago
fixed constant-density pressure extrapolation, respect integer working precision
Property svn:keywords set to `Id`
File size: 204.4 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |