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

Last change on this file since 3997 was 3997, checked in by eckhard, 2 years ago

inifor: Read origin_z from static driver if given; alert user to warnings

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