Home

Context Navigation

inifor_grid.f90 @ 3613

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