Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |