source: palm/trunk/UTIL/inifor/src/grid.f90 @ 3395

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

inifor: Added computation of geostrophic winds from COSMO input

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