Home

Context Navigation

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

Last change on this file since 3615 was 3615, checked in by raasch, 5 years ago
bugfix for last commit: abort replaced by inifor_abort
Property svn:keywords set to `Id`
File size: 212.4 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |