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

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

inifor: Renamed source files for compatibilty with PALM build system

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