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

Last change on this file since 3678 was 3678, checked in by eckhard, 3 years ago

inifor: bugfix: avoid empty averaging regions, check if all input files are present

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