Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |