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

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

inifor: Updated documentation

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