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

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

inifor: bugfix: removed dependency on soilmoisture input files; added netcdf preprocessor flag

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