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

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

inifor: Average initial profiles over the PALM, not the geostrophic, region

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