Home

Context Navigation

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

Last change on this file since 3801 was 3801, checked in by eckhard, 5 years ago
inifor: Read COSMO rotated pole from dataset, check if PALM domain exceeds COSMO domain
Property svn:keywords set to `Id`
File size: 211.5 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |