Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |