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

Last change on this file since 4019 was 4019, checked in by eckhard, 2 years ago

inifor: bugfix: Fix wrong z_top height with gfortran -O2/-O3

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