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

Last change on this file since 3698 was 3680, checked in by knoop, 6 years ago

Added cpp-option netcdf to inifor

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