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

Last change on this file since 3182 was 3182, checked in by suehring, 3 years ago

New Inifor features: grid stretching, improved command-interface, support start dates in different formats in both YYYYMMDD and YYYYMMDDHH, Ability to manually control input file prefixes (--radiation-prefix, --soil-preifx, --flow-prefix, --soilmoisture-prefix) for compatiblity with DWD forcast naming scheme; GNU-style short and long option; Prepared output of large-scale forcing profiles (no computation yet); Added preprocessor flag netcdf4 to switch output format between netCDF 3 and 4; Updated netCDF variable names and attributes to comply with PIDS v1.9; Inifor bugfixes: Improved compatibility with older Intel Intel compilers by avoiding implicit array allocation; Added origin_lon/_lat values and correct reference time in dynamic driver global attributes; corresponding PALM changes: adjustments to revised Inifor; variables names in dynamic driver adjusted; enable geostrophic forcing also in offline nested mode; variable names in LES-LES and COSMO offline nesting changed; lateral boundary flags for nesting, in- and outflow conditions renamed

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