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

Last change on this file since 3802 was 3802, checked in by raasch, 3 years ago

unused variables removed, unused subroutines commented out, type conversion added to avoid compiler warning about constant integer division truncation, script document_changes made bash compatible

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