Home

Context Navigation

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

Last change on this file since 3790 was 3785, checked in by eckhard, 6 years ago
inifor: Removed unused variables, improved coding style
Property svn:keywords set to `Id`
File size: 212.3 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |