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

Last change on this file since 2718 was 2718, checked in by maronga, 4 years ago

deleting of deprecated files; headers updated where needed

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