Home

Context Navigation

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

Last change on this file since 3182 was 3182, checked in by suehring, 6 years ago
New Inifor features: grid stretching, improved command-interface, support start dates in different formats in both YYYYMMDD and YYYYMMDDHH, Ability to manually control input file prefixes (--radiation-prefix, --soil-preifx, --flow-prefix, --soilmoisture-prefix) for compatiblity with DWD forcast naming scheme; GNU-style short and long option; Prepared output of large-scale forcing profiles (no computation yet); Added preprocessor flag netcdf4 to switch output format between netCDF 3 and 4; Updated netCDF variable names and attributes to comply with PIDS v1.9; Inifor bugfixes: Improved compatibility with older Intel Intel compilers by avoiding implicit array allocation; Added origin_lon/_lat values and correct reference time in dynamic driver global attributes; corresponding PALM changes: adjustments to revised Inifor; variables names in dynamic driver adjusted; enable geostrophic forcing also in offline nested mode; variable names in LES-LES and COSMO offline nesting changed; lateral boundary flags for nesting, in- and outflow conditions renamed
Property svn:keywords set to `Id`
File size: 178.7 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |