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

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

inifor: Removed unused variables, improved coding style

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