Home

Context Navigation

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

Last change on this file since 4011 was 3997, checked in by eckhard, 6 years ago
inifor: Read origin_z from static driver if given; alert user to warnings
Property svn:keywords set to `Id`
File size: 206.1 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |