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

Last change on this file since 3456 was 3456, checked in by eckhard, 4 years ago

inifor: Removed surface forcing and internal arrays from netCDF output

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