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

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

inifor: Read COSMO rotated pole from dataset, check if PALM domain exceeds COSMO domain

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