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

Last change on this file since 3537 was 3537, checked in by eckhard, 6 years ago

inifor: COSMO-D2 support

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