Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |