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

Last change on this file since 4475 was 4475, checked in by gronemeier, 18 months ago

bugfixes/changes to INIFOR:

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