Home

Context Navigation

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

Last change on this file since 4725 was 4715, checked in by eckhard, 4 years ago
inifor: Make LODs for top boundary conditions consistent with PALM
Property svn:keywords set to `Id`
File size: 239.0 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |