Home

Context Navigation

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

Last change on this file since 3759 was 3716, checked in by eckhard, 6 years ago
inifor: bugfix: removed dependency on soilmoisture input files; added netcdf preprocessor flag
Property svn:keywords set to `Id`
File size: 211.6 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |