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

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

inifor: Changed initialization mode from 'volume' to 'profile'

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