Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |