Home

Context Navigation

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

Last change on this file since 4481 was 4481, checked in by maronga, 4 years ago
Bugfix for copyright updates in document_changes; copyright update applied to all files
Property svn:keywords set to `Id`
File size: 206.7 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |