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

Last change on this file since 3255 was 3183, checked in by suehring, 6 years ago

last commit documented

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