Home

Context Navigation

grid.f90 @ 3401

Last change on this file since 3401 was 3395, checked in by eckhard, 6 years ago
inifor: Added computation of geostrophic winds from COSMO input
Property svn:keywords set to `Id`
File size: 204.1 KB

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

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

Context Navigation

source: palm/trunk/UTIL/inifor/src/grid.f90 @ 3401

Download in other formats: