Home

Context Navigation

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

Last change on this file since 2704 was 2696, checked in by kanani, 7 years ago
Merge of branch palm4u into trunk
Property svn:keywords set to `Id`
File size: 153.5 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |