Home

Context Navigation

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

Last change on this file since 3557 was 3557, checked in by eckhard, 5 years ago
inifor: Updated documentation
Property svn:keywords set to `Id`
File size: 209.9 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |