Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |