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

Last change on this file since 3614 was 3614, checked in by raasch, 5 years ago

unused variables removed, abort renamed inifor_abort to avoid intrinsic problem in Fortran

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