Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |