Home

Context Navigation

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

Last change on this file since 3556 was 3537, checked in by eckhard, 6 years ago
inifor: COSMO-D2 support
Property svn:keywords set to `Id`
File size: 206.3 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |