Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |