Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |