Home

Context Navigation

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

Last change on this file since 4372 was 4074, checked in by eckhard, 5 years ago
inifor: Changed initialization mode from 'volume' to 'profile'
Property svn:keywords set to `Id`
File size: 206.4 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |