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

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

inifor: Prefixed all INIFOR modules with inifor_ and removed unused variables

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