source: palm/trunk/SOURCE/netcdf_data_input_mod.f90 @ 4598

Last change on this file since 4598 was 4507, checked in by gronemeier, 12 months ago

netcdf_data_input_mod:

  • bugfix: check terrain height for fill values directly after reading (PA0550)
  • changes:
    • remove check for negative zt (PA0551)
    • add reference height from input file to PALM reference height (origin_z)

init_grid:

  • update origin_z with shifting height of orography (oro_min)
  • remove message PA0505
  • Property svn:keywords set to Id
File size: 257.5 KB
Line 
1!> @file netcdf_data_input_mod.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 1997-2020 Leibniz Universitaet Hannover
18!------------------------------------------------------------------------------!
19!
20! Current revisions:
21! -----------------
22!
23!
24! Former revisions:
25! -----------------
26! $Id: netcdf_data_input_mod.f90 4507 2020-04-22 18:21:45Z suehring $
27! - bugfix: check terrain height for fill values directly after reading
28! - changes:
29!   - remove check for negative zt
30!   - add reference height from input file upon PALM reference height (origin_z)
31!
32! 4457 2020-03-11 14:20:43Z raasch
33! use statement for exchange horiz added,
34! bugfixes for calls of exchange horiz 2d
35!
36! 4435 2020-03-03 10:38:41Z raasch
37! temporary bugfix to avoid compile problems with older NetCDFD libraries on IMUK machines
38!
39! 4434 2020-03-03 10:02:18Z oliver.maas
40! added optional netcdf data input for wtm array input parameters
41!
42! 4404 2020-02-12 17:01:53Z suehring
43! Fix misplaced preprocessor directives.
44!
45! 4401 2020-02-11 16:19:09Z suehring
46! Define a default list of coordinate reference system variables used when
47! no static driver input is available
48!
49! 4400 2020-02-10 20:32:41Z suehring
50! - Routine to inquire default fill values added
51! - netcdf_data_input_att and netcdf_data_input_var routines removed
52!
53! 4392 2020-01-31 16:14:57Z pavelkrc
54! (resler) Decrease length of reading buffer (fix problem of ifort/icc compilers)
55!
56! 4389 2020-01-29 08:22:42Z raasch
57! Error messages refined for reading ASCII topo file, also reading of topo file revised so that
58! statement labels and goto statements are not required any more
59!
60! 4388 2020-01-28 16:36:55Z raasch
61! bugfix for error messages while reading ASCII topo file
62!
63! 4387 2020-01-28 11:44:20Z banzhafs
64! Added subroutine get_variable_string_generic ( )
65! and added to interface get_variable to circumvent
66! unknown application-specific restrictions
67! in existing function get_variable_string ( ),
68! which is retained for backward compatibility (ECC)
69!
70! 4370 2020-01-10 14:00:44Z raasch
71! collective read switched off on NEC Aurora to avoid hang situations
72!
73! 4362 2020-01-07 17:15:02Z suehring
74! Input of plant canopy variables from static driver moved to plant-canopy
75! model
76!
77! 4360 2020-01-07 11:25:50Z suehring
78! Correct single message calls, local checks must be given by the respective
79! mpi rank.
80!
81! 4346 2019-12-18 11:55:56Z motisi
82! Introduction of wall_flags_total_0, which currently sets bits based on static
83! topography information used in wall_flags_static_0
84!
85! 4329 2019-12-10 15:46:36Z motisi
86! Renamed wall_flags_0 to wall_flags_static_0
87!
88! 4321 2019-12-04 10:26:38Z pavelkrc
89! Further revise check for surface fractions
90!
91! 4313 2019-11-27 14:07:00Z suehring
92! Checks for surface fractions revised
93!
94! 4312 2019-11-27 14:06:25Z suehring
95! Open input files with read-only attribute instead of write attribute.
96!
97! 4280 2019-10-29 14:34:15Z monakurppa
98! Remove id_emis flags from get_variable_4d_to_3d_real and
99! get_variable_5d_to_4d_real
100!
101! 4258 2019-10-07 13:29:08Z suehring
102! - Migrate input of soil temperature and moisture to land-surface model.
103! - Remove interpolate routines and move the only required subroutine to
104!   land-surface model.
105!
106! 4247 2019-09-30 10:18:24Z pavelkrc
107! Add reading and processing of building_surface_pars
108!
109! 4226 2019-09-10 17:03:24Z suehring
110! - Netcdf input routine for dimension length renamed
111! - Move offline-nesting-specific checks to nesting_offl_mod
112! - Module-specific input of boundary data for offline nesting moved to
113!   nesting_offl_mod
114! - Define module specific data type for offline nesting in nesting_offl_mod
115!
116! 4190 2019-08-27 15:42:37Z suehring
117! type real_1d changed to real_1d_3d
118!
119! 4186 2019-08-23 16:06:14Z suehring
120! Minor formatting adjustments
121!
122! 4182 2019-08-22 15:20:23Z scharf
123! Corrected "Former revisions" section
124!
125! 4178 2019-08-21 11:13:06Z suehring
126! Implement input of external radiation forcing. Therefore, provide public
127! subroutines and variables.
128!
129! 4150 2019-08-08 20:00:47Z suehring
130! Some variables are given the public attribute, in order to call netcdf input
131! from single routines
132!
133! 4125 2019-07-29 13:31:44Z suehring
134! To enable netcdf-parallel access for lateral boundary data (dynamic input),
135! zero number of elements are passed to the respective get_variable routine
136! for non-boundary cores.
137!
138! 4100 2019-07-17 08:11:29Z forkel
139! Made check for input_pids_dynamic and 'inifor' more general
140!
141! 4012 2019-05-31 15:19:05Z monakurppa
142!
143! 3994 2019-05-22 18:08:09Z suehring
144! Remove single location message
145!
146! 3976 2019-05-15 11:02:34Z hellstea
147! Remove unused variables from last commit
148!
149! 3969 2019-05-13 12:14:33Z suehring
150! - clean-up index notations for emission_values to eliminate magic numbers
151! - introduce temporary variable dum_var_5d as well as subroutines
152!   get_var_5d_real and get_var_5d_real_dynamic
153! - remove emission-specific code in generic get_variable routines
154! - in subroutine netcdf_data_input_chemistry_data change netCDF LOD 1
155!   (default) emission_values to the following index order:
156!   z, y, x, species, category
157! - in subroutine netcdf_data_input_chemistry_data
158!   changed netCDF LOD 2 pre-processed emission_values to the following index
159!   order: time, z, y, x, species
160! - in type chem_emis_att_type replace nspec with n_emiss_species
161!   but retained nspec for backward compatibility with salsa_mod. (E.C. Chan)
162!
163! 3961 2019-05-08 16:12:31Z suehring
164! Revise checks for building IDs and types
165!
166! 3943 2019-05-02 09:50:41Z maronga
167! Temporarily disabled some (faulty) checks for static driver.
168!
169! 3942 2019-04-30 13:08:30Z kanani
170! Fix: increase LEN of all NetCDF attribute values (caused crash in
171! netcdf_create_global_atts due to insufficient length)
172!
173! 3941 2019-04-30 09:48:33Z suehring
174! Move check for grid dimension to an earlier point in time when first array
175! is read.
176! Improve checks for building types / IDs with respect to 2D/3D buildings.
177!
178! 3885 2019-04-11 11:29:34Z kanani
179! Changes related to global restructuring of location messages and introduction
180! of additional debug messages
181!
182! 3864 2019-04-05 09:01:56Z monakurppa
183! get_variable_4d_to_3d_real modified to enable read in data of type
184! data(t,y,x,n) one timestep at a time + some routines made public
185!
186! 3855 2019-04-03 10:00:59Z suehring
187! Typo removed
188!
189! 3854 2019-04-02 16:59:33Z suehring
190! Bugfix in one of the checks. Typo removed.
191!
192! 3744 2019-02-15 18:38:58Z suehring
193! Enable mesoscale offline nesting for chemistry variables as well as
194! initialization of chemistry via dynamic input file.
195!
196! 3705 2019-01-29 19:56:39Z suehring
197! Interface for attribute input of 8-bit and 32-bit integer
198!
199! 3704 2019-01-29 19:51:41Z suehring
200! unused variables removed
201!
202! 2696 2017-12-14 17:12:51Z kanani
203! Initial revision (suehring)
204!
205! Authors:
206! --------
207! @author Matthias Suehring
208! @author Edward C. Chan
209! @author Emanuele Russo
210!
211! Description:
212! ------------
213!> Modulue contains routines to input data according to Palm input data
214!> standart using dynamic and static input files.
215!> @todo - Chemistry: revise reading of netcdf file and ajdust formatting
216!>         according to standard!!! (ecc/done)
217!> @todo - Order input alphabetically
218!> @todo - Revise error messages and error numbers
219!> @todo - Input of missing quantities (chemical species, emission rates)
220!> @todo - Defninition and input of still missing variable attributes
221!>         (ecc/what are they?)
222!> @todo - Input of initial geostrophic wind profiles with cyclic conditions.
223!> @todo - remove z dimension from default_emission_data nad preproc_emission_data
224!          and correpsonding subroutines get_var_5d_real and get_var_5d_dynamic (ecc)
225!> @todo - decpreciate chem_emis_att_type@nspec (ecc)
226!> @todo - depreciate subroutines get_variable_4d_to_3d_real and
227!>         get_variable_5d_to_4d_real (ecc)
228!> @todo - introduce useful debug_message(s)
229!------------------------------------------------------------------------------!
230 MODULE netcdf_data_input_mod
231
232    USE control_parameters,                                                    &
233        ONLY:  coupling_char, io_blocks, io_group
234
235    USE cpulog,                                                                &
236        ONLY:  cpu_log, log_point_s
237
238    USE indices,                                                               &
239        ONLY:  nbgp
240
241    USE kinds
242
243#if defined ( __netcdf )
244    USE NETCDF
245#endif
246
247    USE pegrid
248
249    USE surface_mod,                                                           &
250        ONLY:  ind_pav_green, ind_veg_wall, ind_wat_win
251!
252!-- Define type for dimensions.
253    TYPE dims_xy
254       INTEGER(iwp) :: nx                             !< dimension length in x
255       INTEGER(iwp) :: ny                             !< dimension length in y
256       INTEGER(iwp) :: nz                             !< dimension length in z
257       REAL(wp), DIMENSION(:), ALLOCATABLE :: x       !< dimension array in x
258       REAL(wp), DIMENSION(:), ALLOCATABLE :: y       !< dimension array in y
259       REAL(wp), DIMENSION(:), ALLOCATABLE :: z       !< dimension array in z
260    END TYPE dims_xy
261    TYPE init_type
262
263       CHARACTER(LEN=16) ::  init_char = 'init_atmosphere_'          !< leading substring for init variables
264       CHARACTER(LEN=23) ::  origin_time = '2000-01-01 00:00:00 +00' !< reference time of input data
265       CHARACTER(LEN=100), DIMENSION(:), ALLOCATABLE ::  var_names_chem !< list of chemistry variable names that can potentially be on file
266
267       INTEGER(iwp) ::  lod_msoil !< level of detail - soil moisture
268       INTEGER(iwp) ::  lod_pt    !< level of detail - pt
269       INTEGER(iwp) ::  lod_q     !< level of detail - q
270       INTEGER(iwp) ::  lod_tsoil !< level of detail - soil temperature
271       INTEGER(iwp) ::  lod_u     !< level of detail - u-component
272       INTEGER(iwp) ::  lod_v     !< level of detail - v-component
273       INTEGER(iwp) ::  lod_w     !< level of detail - w-component
274       INTEGER(iwp) ::  nx        !< number of scalar grid points along x in dynamic input file
275       INTEGER(iwp) ::  nxu       !< number of u grid points along x in dynamic input file
276       INTEGER(iwp) ::  ny        !< number of scalar grid points along y in dynamic input file
277       INTEGER(iwp) ::  nyv       !< number of v grid points along y in dynamic input file
278       INTEGER(iwp) ::  nzs       !< number of vertical soil levels in dynamic input file
279       INTEGER(iwp) ::  nzu       !< number of vertical levels on scalar grid in dynamic input file
280       INTEGER(iwp) ::  nzw       !< number of vertical levels on w grid in dynamic input file
281       
282       INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  lod_chem !< level of detail - chemistry variables
283
284       LOGICAL ::  from_file_msoil  = .FALSE. !< flag indicating whether soil moisture is already initialized from file
285       LOGICAL ::  from_file_pt     = .FALSE. !< flag indicating whether pt is already initialized from file
286       LOGICAL ::  from_file_q      = .FALSE. !< flag indicating whether q is already initialized from file
287       LOGICAL ::  from_file_tsoil  = .FALSE. !< flag indicating whether soil temperature is already initialized from file
288       LOGICAL ::  from_file_u      = .FALSE. !< flag indicating whether u is already initialized from file
289       LOGICAL ::  from_file_ug     = .FALSE. !< flag indicating whether ug is already initialized from file
290       LOGICAL ::  from_file_v      = .FALSE. !< flag indicating whether v is already initialized from file
291       LOGICAL ::  from_file_vg     = .FALSE. !< flag indicating whether ug is already initialized from file
292       LOGICAL ::  from_file_w      = .FALSE. !< flag indicating whether w is already initialized from file
293       
294       LOGICAL, DIMENSION(:), ALLOCATABLE ::  from_file_chem !< flag indicating whether chemistry variable is read from file
295
296       REAL(wp) ::  fill_msoil              !< fill value for soil moisture
297       REAL(wp) ::  fill_pt                 !< fill value for pt
298       REAL(wp) ::  fill_q                  !< fill value for q
299       REAL(wp) ::  fill_tsoil              !< fill value for soil temperature
300       REAL(wp) ::  fill_u                  !< fill value for u
301       REAL(wp) ::  fill_v                  !< fill value for v
302       REAL(wp) ::  fill_w                  !< fill value for w
303       REAL(wp) ::  latitude = 0.0_wp       !< latitude of the lower left corner
304       REAL(wp) ::  longitude = 0.0_wp      !< longitude of the lower left corner
305       REAL(wp) ::  origin_x = 500000.0_wp  !< UTM easting of the lower left corner
306       REAL(wp) ::  origin_y = 0.0_wp       !< UTM northing of the lower left corner
307       REAL(wp) ::  origin_z = 0.0_wp       !< reference height of input data
308       REAL(wp) ::  rotation_angle = 0.0_wp !< rotation angle of input data
309
310       REAL(wp), DIMENSION(:), ALLOCATABLE ::  fill_chem    !< fill value - chemistry variables
311       REAL(wp), DIMENSION(:), ALLOCATABLE ::  msoil_1d     !< initial vertical profile of soil moisture
312       REAL(wp), DIMENSION(:), ALLOCATABLE ::  pt_init      !< initial vertical profile of pt
313       REAL(wp), DIMENSION(:), ALLOCATABLE ::  q_init       !< initial vertical profile of q
314       REAL(wp), DIMENSION(:), ALLOCATABLE ::  tsoil_1d     !< initial vertical profile of soil temperature
315       REAL(wp), DIMENSION(:), ALLOCATABLE ::  u_init       !< initial vertical profile of u
316       REAL(wp), DIMENSION(:), ALLOCATABLE ::  ug_init      !< initial vertical profile of ug
317       REAL(wp), DIMENSION(:), ALLOCATABLE ::  v_init       !< initial vertical profile of v
318       REAL(wp), DIMENSION(:), ALLOCATABLE ::  vg_init      !< initial vertical profile of ug
319       REAL(wp), DIMENSION(:), ALLOCATABLE ::  w_init       !< initial vertical profile of w
320       REAL(wp), DIMENSION(:), ALLOCATABLE ::  z_soil       !< vertical levels in soil in dynamic input file, used for interpolation
321       REAL(wp), DIMENSION(:), ALLOCATABLE ::  zu_atmos     !< vertical levels at scalar grid in dynamic input file, used for interpolation
322       REAL(wp), DIMENSION(:), ALLOCATABLE ::  zw_atmos     !< vertical levels at w grid in dynamic input file, used for interpolation
323       
324       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  chem_init  !< initial vertical profiles of chemistry variables
325
326
327       REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  msoil_3d !< initial 3d soil moisture provide by Inifor and interpolated onto soil grid
328       REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  tsoil_3d !< initial 3d soil temperature provide by Inifor and interpolated onto soil grid
329
330    END TYPE init_type
331
332!-- Data type for the general information of chemistry emissions, do not dependent on the particular chemical species
333    TYPE chem_emis_att_type
334
335       !-DIMENSIONS
336       
337       INTEGER(iwp)                                 :: nspec=0            !< no of chem species provided in emission_values
338       INTEGER(iwp)                                 :: n_emiss_species=0  !< no of chem species provided in emission_values
339                                                                          !< same function as nspec, which will be depreciated (ecc)
340                                                                                 
341       INTEGER(iwp)                                 :: ncat=0             !< number of emission categories
342       INTEGER(iwp)                                 :: nvoc=0             !< number of VOC components
343       INTEGER(iwp)                                 :: npm=0              !< number of PM components
344       INTEGER(iwp)                                 :: nnox=2             !< number of NOx components: NO and NO2
345       INTEGER(iwp)                                 :: nsox=2             !< number of SOX components: SO and SO4
346       INTEGER(iwp)                                 :: nhoursyear         !< number of hours of a specific year in the HOURLY mode
347                                                                          !< of the default mode
348       INTEGER(iwp)                                 :: nmonthdayhour      !< number of month days and hours in the MDH mode
349                                                                          !< of the default mode
350       INTEGER(iwp)                                 :: dt_emission        !< Number of emissions timesteps for one year
351                                                                          !< in the pre-processed emissions case
352       !-- 1d emission input variables
353       CHARACTER (LEN=25),ALLOCATABLE, DIMENSION(:) :: pm_name       !< Names of PM components
354       CHARACTER (LEN=25),ALLOCATABLE, DIMENSION(:) :: cat_name      !< Emission category names
355       CHARACTER (LEN=25),ALLOCATABLE, DIMENSION(:) :: species_name  !< Names of emission chemical species
356       CHARACTER (LEN=25),ALLOCATABLE, DIMENSION(:) :: voc_name      !< Names of VOCs components
357       CHARACTER (LEN=25)                           :: units         !< Units
358
359       INTEGER(iwp)                                 :: i_hour         !< indices for assigning emission values at different timesteps
360       INTEGER(iwp),ALLOCATABLE, DIMENSION(:)       :: cat_index      !< Indices for emission categories
361       INTEGER(iwp),ALLOCATABLE, DIMENSION(:)       :: species_index  !< Indices for emission chem species
362
363       REAL(wp),ALLOCATABLE, DIMENSION(:)           :: xm             !< Molecular masses of emission chem species
364
365       !-- 2d emission input variables
366       REAL(wp),ALLOCATABLE, DIMENSION(:,:)         :: hourly_emis_time_factor  !< Time factors for HOURLY emissions (DEFAULT mode)
367       REAL(wp),ALLOCATABLE, DIMENSION(:,:)         :: mdh_emis_time_factor     !< Time factors for MDH emissions (DEFAULT mode)
368       REAL(wp),ALLOCATABLE, DIMENSION(:,:)         :: nox_comp                 !< Composition of NO and NO2
369       REAL(wp),ALLOCATABLE, DIMENSION(:,:)         :: sox_comp                 !< Composition of SO2 and SO4
370       REAL(wp),ALLOCATABLE, DIMENSION(:,:)         :: voc_comp                 !< Composition of VOC components (not fixed)
371
372       !-- 3d emission input variables
373       REAL(wp),ALLOCATABLE, DIMENSION(:,:,:)       :: pm_comp                  !< Composition of PM components (not fixed)
374 
375    END TYPE chem_emis_att_type
376
377
378!-- Data type for the values of chemistry emissions
379    TYPE chem_emis_val_type
380
381       !REAL(wp),ALLOCATABLE, DIMENSION(:,:)     :: stack_height           !< stack height (ecc / to be implemented)
382       REAL(wp),ALLOCATABLE, DIMENSION(:,:,:)    :: default_emission_data  !< Emission input values for LOD1 (DEFAULT mode)
383       REAL(wp),ALLOCATABLE, DIMENSION(:,:,:,:)  :: preproc_emission_data  !< Emission input values for LOD2 (PRE-PROCESSED mode)
384
385    END TYPE chem_emis_val_type
386
387!
388!-- Define data structures for different input data types.
389!-- 8-bit Integer 2D
390    TYPE int_2d_8bit
391       INTEGER(KIND=1) ::  fill = -127                      !< fill value
392       INTEGER(KIND=1), DIMENSION(:,:), ALLOCATABLE ::  var !< respective variable
393
394       LOGICAL ::  from_file = .FALSE.  !< flag indicating whether an input variable is available and read from file or default values are used
395    END TYPE int_2d_8bit
396!
397!-- 8-bit Integer 3D
398    TYPE int_3d_8bit
399       INTEGER(KIND=1) ::  fill = -127                           !< fill value
400       INTEGER(KIND=1), DIMENSION(:,:,:), ALLOCATABLE ::  var_3d !< respective variable
401
402       LOGICAL ::  from_file = .FALSE.  !< flag indicating whether an input variable is available and read from file or default values are used
403    END TYPE int_3d_8bit
404!
405!-- 32-bit Integer 2D
406    TYPE int_2d_32bit
407       INTEGER(iwp) ::  fill = -9999                      !< fill value
408       INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE ::  var  !< respective variable
409
410       LOGICAL ::  from_file = .FALSE. !< flag indicating whether an input variable is available and read from file or default values are used
411    END TYPE int_2d_32bit
412!
413!-- Define data type to read 1D or 3D real variables.
414    TYPE real_1d_3d
415       LOGICAL ::  from_file = .FALSE.  !< flag indicating whether an input variable is available and read from file or default values are used
416
417       INTEGER(iwp) ::  lod = -1        !< level-of-detail
418       
419       REAL(wp) ::  fill = -9999.9_wp                  !< fill value
420       
421       REAL(wp), DIMENSION(:),     ALLOCATABLE ::  var1d     !< respective 1D variable
422       REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  var3d     !< respective 3D variable
423    END TYPE real_1d_3d   
424!
425!-- Define data type to read 2D real variables
426    TYPE real_2d
427       LOGICAL ::  from_file = .FALSE.  !< flag indicating whether an input variable is available and read from file or default values are used
428
429       INTEGER(iwp) ::  lod             !< level-of-detail
430       
431       REAL(wp) ::  fill = -9999.9_wp                !< fill value
432       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  var !< respective variable
433    END TYPE real_2d
434
435!
436!-- Define data type to read 3D real variables
437    TYPE real_3d
438       LOGICAL ::  from_file = .FALSE.  !< flag indicating whether an input variable is available and read from file or default values are used
439
440       INTEGER(iwp) ::  nz   !< number of grid points along vertical dimension
441
442       REAL(wp) ::  fill = -9999.9_wp                  !< fill value
443       REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  var !< respective variable
444    END TYPE real_3d
445!
446!-- Define data structure where the dimension and type of the input depends
447!-- on the given level of detail.
448!-- For buildings, the input is either 2D float, or 3d byte.
449    TYPE build_in
450       INTEGER(iwp)    ::  lod = 1                               !< level of detail
451       INTEGER(KIND=1) ::  fill2 = -127                          !< fill value for lod = 2
452       INTEGER(iwp)    ::  nz                                    !< number of vertical layers in file
453       INTEGER(KIND=1), DIMENSION(:,:,:), ALLOCATABLE ::  var_3d !< 3d variable (lod = 2)
454
455       REAL(wp), DIMENSION(:), ALLOCATABLE ::  z                 !< vertical coordinate for 3D building, used for consistency check
456
457       LOGICAL ::  from_file = .FALSE.  !< flag indicating whether an input variable is available and read from file or default values are used
458
459       REAL(wp)                              ::  fill1 = -9999.9_wp !< fill values for lod = 1
460       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  var_2d             !< 2d variable (lod = 1)
461       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  oro_max            !< terraing height under particular buildings
462    END TYPE build_in
463
464!
465!-- For soil_type, the input is either 2D or 3D one-byte integer.
466    TYPE soil_in
467       INTEGER(iwp)                                   ::  lod = 1      !< level of detail
468       INTEGER(KIND=1)                                ::  fill = -127  !< fill value for lod = 2
469       INTEGER(KIND=1), DIMENSION(:,:), ALLOCATABLE   ::  var_2d       !< 2d variable (lod = 1)
470       INTEGER(KIND=1), DIMENSION(:,:,:), ALLOCATABLE ::  var_3d       !< 3d variable (lod = 2)
471
472       LOGICAL ::  from_file = .FALSE.  !< flag indicating whether an input variable is available and read from file or default values are used
473    END TYPE soil_in
474
475!
476!-- Define data type for fractions between surface types
477    TYPE fracs
478       INTEGER(iwp)                            ::  nf             !< total number of fractions
479       INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  nfracs         !< dimension array for fraction
480
481       LOGICAL ::  from_file = .FALSE. !< flag indicating whether an input variable is available and read from file or default values are used
482
483       REAL(wp)                                ::  fill = -9999.9_wp !< fill value
484       REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  frac              !< respective fraction between different surface types
485    END TYPE fracs
486!
487!-- Data type for parameter lists, Depending on the given level of detail,
488!-- the input is 3D or 4D
489    TYPE pars
490       INTEGER(iwp)                            ::  lod = 1         !< level of detail
491       INTEGER(iwp)                            ::  np              !< total number of parameters
492       INTEGER(iwp)                            ::  nz              !< vertical dimension - number of soil layers
493       INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  layers          !< dimension array for soil layers
494       INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  pars            !< dimension array for parameters
495
496       LOGICAL ::  from_file = .FALSE.  !< flag indicating whether an input variable is available and read from file or default values are used
497
498       REAL(wp)                                  ::  fill = -9999.9_wp !< fill value
499       REAL(wp), DIMENSION(:,:,:), ALLOCATABLE   ::  pars_xy           !< respective parameters, level of detail = 1
500       REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE ::  pars_xyz          !< respective parameters, level of detail = 2
501    END TYPE pars
502!
503!-- Data type for surface parameter lists
504    TYPE pars_surf
505       INTEGER(iwp)                                ::  np          !< total number of parameters
506       INTEGER(iwp)                                ::  nsurf       !< number of local surfaces
507       INTEGER(iwp), DIMENSION(:,:,:), ALLOCATABLE ::  index_ji    !< index for beginning and end of surfaces at (j,i)
508       INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE   ::  coords      !< (k,j,i,norm_z,norm_y,norm_x)
509                                                                   !< k,j,i:                surface position
510                                                                   !< norm_z,norm_y,norm_x: surface normal vector
511
512       LOGICAL ::  from_file = .FALSE.  !< flag indicating whether an input variable is available and read from file or default values are used
513
514       REAL(wp)                              ::  fill = -9999.9_wp !< fill value
515       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  pars              !< respective parameters per surface
516    END TYPE pars_surf
517!
518!-- Define type for global file attributes
519!-- Please refer to the PALM data standard for a detailed description of each
520!-- attribute.
521    TYPE global_atts_type
522       CHARACTER(LEN=200) ::  acronym = ' '                      !< acronym of institution
523       CHARACTER(LEN=7)   ::  acronym_char = 'acronym'           !< name of attribute
524       CHARACTER(LEN=200) ::  author  = ' '                      !< first name, last name, email adress
525       CHARACTER(LEN=6)   ::  author_char = 'author'             !< name of attribute
526       CHARACTER(LEN=200) ::  campaign = 'PALM-4U'               !< name of campaign
527       CHARACTER(LEN=8)   ::  campaign_char = 'campaign'         !< name of attribute
528       CHARACTER(LEN=200) ::  comment = ' '                      !< comment to data
529       CHARACTER(LEN=7)   ::  comment_char = 'comment'           !< name of attribute
530       CHARACTER(LEN=200) ::  contact_person = ' '               !< first name, last name, email adress
531       CHARACTER(LEN=14)  ::  contact_person_char = 'contact_person'  !< name of attribute
532       CHARACTER(LEN=200) ::  conventions = 'CF-1.7'             !< netCDF convention
533       CHARACTER(LEN=11)  ::  conventions_char = 'Conventions'   !< name of attribute
534       CHARACTER(LEN=23 ) ::  creation_time = ' '                !< creation time of data set
535       CHARACTER(LEN=13)  ::  creation_time_char = 'creation_time'  !< name of attribute
536       CHARACTER(LEN=200) ::  data_content = ' '                 !< content of data set
537       CHARACTER(LEN=12)  ::  data_content_char = 'data_content' !< name of attribute
538       CHARACTER(LEN=200) ::  dependencies = ' '                 !< dependencies of data set
539       CHARACTER(LEN=12)  ::  dependencies_char = 'dependencies' !< name of attribute
540       CHARACTER(LEN=200) ::  history = ' '                      !< information about data processing
541       CHARACTER(LEN=7)   ::  history_char = 'history'           !< name of attribute
542       CHARACTER(LEN=200) ::  institution = ' '                  !< name of responsible institution
543       CHARACTER(LEN=11)  ::  institution_char = 'institution'   !< name of attribute
544       CHARACTER(LEN=200) ::  keywords = ' '                     !< keywords of data set
545       CHARACTER(LEN=8)   ::  keywords_char = 'keywords'         !< name of attribute
546       CHARACTER(LEN=200) ::  licence = ' '                      !< licence of data set
547       CHARACTER(LEN=7)   ::  licence_char = 'licence'           !< name of attribute
548       CHARACTER(LEN=200) ::  location = ' '                     !< place which refers to data set
549       CHARACTER(LEN=8)   ::  location_char = 'location'         !< name of attribute
550       CHARACTER(LEN=10)  ::  origin_lat_char = 'origin_lat'     !< name of attribute
551       CHARACTER(LEN=10)  ::  origin_lon_char = 'origin_lon'     !< name of attribute
552       CHARACTER(LEN=23 ) ::  origin_time = '2000-01-01 00:00:00 +00'  !< reference time
553       CHARACTER(LEN=11)  ::  origin_time_char = 'origin_time'   !< name of attribute
554       CHARACTER(LEN=8)   ::  origin_x_char = 'origin_x'         !< name of attribute
555       CHARACTER(LEN=8)   ::  origin_y_char = 'origin_y'         !< name of attribute
556       CHARACTER(LEN=8)   ::  origin_z_char = 'origin_z'         !< name of attribute
557       CHARACTER(LEN=12)  ::  palm_version_char = 'palm_version' !< name of attribute
558       CHARACTER(LEN=200) ::  references = ' '                   !< literature referring to data set
559       CHARACTER(LEN=10)  ::  references_char = 'references'     !< name of attribute
560       CHARACTER(LEN=14)  ::  rotation_angle_char = 'rotation_angle'  !< name of attribute
561       CHARACTER(LEN=200) ::  site = ' '                         !< name of model domain
562       CHARACTER(LEN=4)   ::  site_char = 'site'                 !< name of attribute
563       CHARACTER(LEN=200) ::  source = ' '                       !< source of data set
564       CHARACTER(LEN=6)   ::  source_char = 'source'             !< name of attribute
565       CHARACTER(LEN=200) ::  title = ' '                        !< title of data set
566       CHARACTER(LEN=5)   ::  title_char = 'title'               !< name of attribute
567       CHARACTER(LEN=7)   ::  version_char = 'version'           !< name of attribute
568
569       INTEGER(iwp) ::  version              !< version of data set
570
571       REAL(wp) ::  fillvalue = -9999.0      !< default fill value
572       REAL(wp) ::  origin_lat               !< latitude of lower left corner
573       REAL(wp) ::  origin_lon               !< longitude of lower left corner
574       REAL(wp) ::  origin_x                 !< easting (UTM coordinate) of lower left corner
575       REAL(wp) ::  origin_y                 !< northing (UTM coordinate) of lower left corner
576       REAL(wp) ::  origin_z                 !< reference height
577       REAL(wp) ::  palm_version             !< PALM version of data set
578       REAL(wp) ::  rotation_angle           !< rotation angle of coordinate system of data set
579    END TYPE global_atts_type
580!
581!-- Define type for coordinate reference system (crs)
582    TYPE crs_type
583       CHARACTER(LEN=200) ::  epsg_code = 'EPSG:25831'                   !< EPSG code
584       CHARACTER(LEN=200) ::  grid_mapping_name = 'transverse_mercator'  !< name of grid mapping
585       CHARACTER(LEN=200) ::  long_name = 'coordinate reference system'  !< name of variable crs
586       CHARACTER(LEN=200) ::  units = 'm'                                !< unit of crs
587
588       REAL(wp) ::  false_easting = 500000.0_wp                  !< false easting
589       REAL(wp) ::  false_northing = 0.0_wp                      !< false northing
590       REAL(wp) ::  inverse_flattening = 298.257223563_wp        !< 1/f (default for WGS84)
591       REAL(wp) ::  latitude_of_projection_origin = 0.0_wp       !< latitude of projection origin
592       REAL(wp) ::  longitude_of_central_meridian = 3.0_wp       !< longitude of central meridian of UTM zone (default: zone 31)
593       REAL(wp) ::  longitude_of_prime_meridian = 0.0_wp         !< longitude of prime meridian
594       REAL(wp) ::  scale_factor_at_central_meridian = 0.9996_wp !< scale factor of UTM coordinates
595       REAL(wp) ::  semi_major_axis = 6378137.0_wp               !< length of semi major axis (default for WGS84)
596    END TYPE crs_type
597
598!
599!-- Define variables
600    TYPE(crs_type)   ::  coord_ref_sys  !< coordinate reference system
601
602    TYPE(dims_xy)    ::  dim_static     !< data structure for x, y-dimension in static input file
603
604    TYPE(init_type) ::  init_3d    !< data structure for the initialization of the 3D flow and soil fields
605    TYPE(init_type) ::  init_model !< data structure for the initialization of the model
606
607!
608!-- Define 2D variables of type NC_BYTE
609    TYPE(int_2d_8bit)  ::  albedo_type_f     !< input variable for albedo type
610    TYPE(int_2d_8bit)  ::  building_type_f   !< input variable for building type
611    TYPE(int_2d_8bit)  ::  pavement_type_f   !< input variable for pavenment type
612    TYPE(int_2d_8bit)  ::  street_crossing_f !< input variable for water type
613    TYPE(int_2d_8bit)  ::  street_type_f     !< input variable for water type
614    TYPE(int_2d_8bit)  ::  vegetation_type_f !< input variable for vegetation type
615    TYPE(int_2d_8bit)  ::  water_type_f      !< input variable for water type
616!
617!-- Define 3D variables of type NC_BYTE
618    TYPE(int_3d_8bit)  ::  building_obstruction_f    !< input variable for building obstruction
619    TYPE(int_3d_8bit)  ::  building_obstruction_full !< input variable for building obstruction
620!
621!-- Define 2D variables of type NC_INT
622    TYPE(int_2d_32bit) ::  building_id_f     !< input variable for building ID
623!
624!-- Define 2D variables of type NC_FLOAT
625    TYPE(real_2d) ::  terrain_height_f       !< input variable for terrain height
626    TYPE(real_2d) ::  uvem_irradiance_f      !< input variable for uvem irradiance lookup table
627    TYPE(real_2d) ::  uvem_integration_f     !< input variable for uvem integration
628!
629!-- Define 3D variables of type NC_FLOAT
630    TYPE(real_3d) ::  root_area_density_lsm_f !< input variable for root area density - parametrized vegetation
631    TYPE(real_3d) ::  uvem_radiance_f         !< input variable for uvem radiance lookup table
632    TYPE(real_3d) ::  uvem_projarea_f         !< input variable for uvem projection area lookup table
633!
634!-- Define input variable for buildings
635    TYPE(build_in) ::  buildings_f           !< input variable for buildings
636!
637!-- Define input variables for soil_type
638    TYPE(soil_in)  ::  soil_type_f           !< input variable for soil type
639
640    TYPE(fracs) ::  surface_fraction_f       !< input variable for surface fraction
641
642    TYPE(pars)  ::  albedo_pars_f              !< input variable for albedo parameters
643    TYPE(pars)  ::  building_pars_f            !< input variable for building parameters
644    TYPE(pars)  ::  pavement_pars_f            !< input variable for pavement parameters
645    TYPE(pars)  ::  pavement_subsurface_pars_f !< input variable for pavement parameters
646    TYPE(pars)  ::  soil_pars_f                !< input variable for soil parameters
647    TYPE(pars)  ::  vegetation_pars_f          !< input variable for vegetation parameters
648    TYPE(pars)  ::  water_pars_f               !< input variable for water parameters
649
650    TYPE(pars_surf)  ::  building_surface_pars_f  !< input variable for building surface parameters
651
652    TYPE(chem_emis_att_type)                             ::  chem_emis_att    !< Input Information of Chemistry Emission Data from netcdf 
653    TYPE(chem_emis_val_type), ALLOCATABLE, DIMENSION(:)  ::  chem_emis        !< Input Chemistry Emission Data from netcdf 
654
655    CHARACTER(LEN=3)  ::  char_lod  = 'lod'         !< name of level-of-detail attribute in NetCDF file
656
657    CHARACTER(LEN=10) ::  char_fill = '_FillValue'        !< name of fill value attribute in NetCDF file
658
659    CHARACTER(LEN=100) ::  input_file_static  = 'PIDS_STATIC'  !< Name of file which comprises static input data
660    CHARACTER(LEN=100) ::  input_file_dynamic = 'PIDS_DYNAMIC' !< Name of file which comprises dynamic input data
661    CHARACTER(LEN=100) ::  input_file_chem    = 'PIDS_CHEM'    !< Name of file which comprises chemistry input data
662    CHARACTER(LEN=100) ::  input_file_uvem    = 'PIDS_UVEM'    !< Name of file which comprises static uv_exposure model input data
663    CHARACTER(LEN=100) ::  input_file_vm      = 'PIDS_VM'      !< Name of file which comprises virtual measurement data
664    CHARACTER(LEN=100) ::  input_file_wtm     = 'PIDS_WTM'     !< Name of file which comprises wind turbine model input data
665       
666    CHARACTER(LEN=25), ALLOCATABLE, DIMENSION(:) ::  string_values  !< output of string variables read from netcdf input files
667    CHARACTER(LEN=50), DIMENSION(:), ALLOCATABLE ::  vars_pids      !< variable in input file
668
669    INTEGER(iwp)                                     ::  id_emis        !< NetCDF id of input file for chemistry emissions: TBD: It has to be removed
670
671    INTEGER(iwp) ::  nc_stat         !< return value of nf90 function call
672    INTEGER(iwp) ::  num_var_pids    !< number of variables in file
673    INTEGER(iwp) ::  pids_id         !< file id
674
675    LOGICAL ::  input_pids_static  = .FALSE.   !< Flag indicating whether Palm-input-data-standard file containing static information exists
676    LOGICAL ::  input_pids_dynamic = .FALSE.   !< Flag indicating whether Palm-input-data-standard file containing dynamic information exists
677    LOGICAL ::  input_pids_chem    = .FALSE.   !< Flag indicating whether Palm-input-data-standard file containing chemistry information exists
678    LOGICAL ::  input_pids_uvem    = .FALSE.   !< Flag indicating whether uv-expoure-model input file containing static information exists
679    LOGICAL ::  input_pids_vm      = .FALSE.   !< Flag indicating whether input file for virtual measurements exist
680    LOGICAL ::  input_pids_wtm     = .FALSE.   !< Flag indicating whether input file for wind turbine model exists
681
682    LOGICAL ::  collective_read = .FALSE.      !< Enable NetCDF collective read
683
684    REAL(wp), DIMENSION(8) ::  crs_list        !< list of coord_ref_sys values
685
686    TYPE(global_atts_type) ::  input_file_atts !< global attributes of input file
687
688    SAVE
689
690    PRIVATE
691
692    INTERFACE netcdf_data_input_check_dynamic
693       MODULE PROCEDURE netcdf_data_input_check_dynamic
694    END INTERFACE netcdf_data_input_check_dynamic
695
696    INTERFACE netcdf_data_input_check_static
697       MODULE PROCEDURE netcdf_data_input_check_static
698    END INTERFACE netcdf_data_input_check_static
699
700    INTERFACE netcdf_data_input_chemistry_data
701       MODULE PROCEDURE netcdf_data_input_chemistry_data
702    END INTERFACE netcdf_data_input_chemistry_data
703
704    INTERFACE get_dimension_length
705       MODULE PROCEDURE get_dimension_length
706    END INTERFACE get_dimension_length
707
708    INTERFACE inquire_fill_value
709       MODULE PROCEDURE inquire_fill_value_int
710       MODULE PROCEDURE inquire_fill_value_real
711    END INTERFACE inquire_fill_value
712
713    INTERFACE netcdf_data_input_inquire_file
714       MODULE PROCEDURE netcdf_data_input_inquire_file
715    END INTERFACE netcdf_data_input_inquire_file
716
717    INTERFACE netcdf_data_input_init
718       MODULE PROCEDURE netcdf_data_input_init
719    END INTERFACE netcdf_data_input_init
720
721    INTERFACE netcdf_data_input_init_3d
722       MODULE PROCEDURE netcdf_data_input_init_3d
723    END INTERFACE netcdf_data_input_init_3d
724   
725    INTERFACE netcdf_data_input_surface_data
726       MODULE PROCEDURE netcdf_data_input_surface_data
727    END INTERFACE netcdf_data_input_surface_data
728
729    INTERFACE netcdf_data_input_uvem
730       MODULE PROCEDURE netcdf_data_input_uvem
731    END INTERFACE netcdf_data_input_uvem
732
733    INTERFACE get_variable
734       MODULE PROCEDURE get_variable_1d_char
735       MODULE PROCEDURE get_variable_1d_int
736       MODULE PROCEDURE get_variable_1d_real
737       MODULE PROCEDURE get_variable_2d_int8
738       MODULE PROCEDURE get_variable_2d_int32
739       MODULE PROCEDURE get_variable_2d_real
740       MODULE PROCEDURE get_variable_3d_int8
741       MODULE PROCEDURE get_variable_3d_real
742       MODULE PROCEDURE get_variable_3d_real_dynamic
743       MODULE PROCEDURE get_variable_4d_to_3d_real
744       MODULE PROCEDURE get_variable_4d_real
745       MODULE PROCEDURE get_variable_5d_to_4d_real
746       MODULE PROCEDURE get_variable_5d_real           ! (ecc) temp subroutine 4 reading 5D NC arrays
747       MODULE PROCEDURE get_variable_5d_real_dynamic   ! 2B removed as z is out of emission_values
748       MODULE PROCEDURE get_variable_string
749       MODULE PROCEDURE get_variable_string_generic    ! (ecc) generic string function
750
751    END INTERFACE get_variable
752
753    INTERFACE get_variable_pr
754       MODULE PROCEDURE get_variable_pr
755    END INTERFACE get_variable_pr
756
757    INTERFACE get_attribute
758       MODULE PROCEDURE get_attribute_real
759       MODULE PROCEDURE get_attribute_int8
760       MODULE PROCEDURE get_attribute_int32
761       MODULE PROCEDURE get_attribute_string
762    END INTERFACE get_attribute
763
764!
765!-- Public data structures
766    PUBLIC real_1d_3d,                                                         &
767           real_2d,                                                            &
768           real_3d
769!
770!-- Public variables
771    PUBLIC albedo_pars_f, albedo_type_f, buildings_f,                          &
772           building_id_f, building_pars_f, building_surface_pars_f,            &
773           building_type_f,                                                    &
774           char_fill,                                                          &
775           char_lod,                                                           &
776           chem_emis, chem_emis_att, chem_emis_att_type, chem_emis_val_type,   &
777           coord_ref_sys,                                                      &
778           crs_list,                                                           &
779           init_3d, init_model, input_file_atts,                               &
780           input_file_dynamic,                                                 &
781           input_file_static,                                                  &
782           input_pids_static,                                                  &
783           input_pids_dynamic, input_pids_vm, input_file_vm,                   &
784           input_pids_wtm, input_file_wtm,                                     &
785           num_var_pids,                                                       &
786           pavement_pars_f, pavement_subsurface_pars_f, pavement_type_f,       &
787           pids_id,                                                            &
788           root_area_density_lsm_f, soil_pars_f,                               &
789           soil_type_f, street_crossing_f, street_type_f, surface_fraction_f,  &
790           terrain_height_f, vegetation_pars_f, vegetation_type_f,             &
791           vars_pids,                                                          &
792           water_pars_f, water_type_f
793!
794!-- Public uv exposure variables
795    PUBLIC building_obstruction_f, input_file_uvem, input_pids_uvem,           &
796           netcdf_data_input_uvem,                                             &
797           uvem_integration_f, uvem_irradiance_f,                              &
798           uvem_projarea_f, uvem_radiance_f
799
800!
801!-- Public subroutines
802    PUBLIC netcdf_data_input_check_dynamic,                                    &
803           netcdf_data_input_check_static,                                     &
804           netcdf_data_input_chemistry_data,                                   &
805           get_dimension_length,                                               &
806           netcdf_data_input_inquire_file,                                     &
807           netcdf_data_input_init,                                             &
808           netcdf_data_input_init_3d,                                          &
809           netcdf_data_input_surface_data,                                     &
810           netcdf_data_input_topo,                                             &
811           get_attribute,                                                      &
812           get_variable,                                                       &
813           get_variable_pr,                                                    &
814           open_read_file,                                                     &
815           check_existence,                                                    &
816           inquire_fill_value,                                                 &
817           inquire_num_variables,                                              &
818           inquire_variable_names,                                             &
819           close_input_file
820
821
822 CONTAINS
823
824!------------------------------------------------------------------------------!
825! Description:
826! ------------
827!> Inquires whether NetCDF input files according to Palm-input-data standard
828!> exist. Moreover, basic checks are performed.
829!------------------------------------------------------------------------------!
830    SUBROUTINE netcdf_data_input_inquire_file
831
832       USE control_parameters,                                                 &
833           ONLY:  topo_no_distinct
834
835       IMPLICIT NONE
836
837#if defined ( __netcdf )
838       INQUIRE( FILE = TRIM( input_file_static )   // TRIM( coupling_char ),   &
839                EXIST = input_pids_static  )
840       INQUIRE( FILE = TRIM( input_file_dynamic ) // TRIM( coupling_char ),    &
841                EXIST = input_pids_dynamic )
842       INQUIRE( FILE = TRIM( input_file_chem )    // TRIM( coupling_char ),    &
843                EXIST = input_pids_chem )
844       INQUIRE( FILE = TRIM( input_file_uvem )    // TRIM( coupling_char ),    &
845                EXIST = input_pids_uvem  )
846       INQUIRE( FILE = TRIM( input_file_vm )      // TRIM( coupling_char ),    &
847                EXIST = input_pids_vm )
848       INQUIRE( FILE = TRIM( input_file_wtm )     // TRIM( coupling_char ),    &
849                EXIST = input_pids_wtm )
850#endif
851
852!
853!--    As long as topography can be input via ASCII format, no distinction
854!--    between building and terrain can be made. This case, classify all
855!--    surfaces as default type. Same in case land-surface and urban-surface
856!--    model are not applied.
857       IF ( .NOT. input_pids_static )  THEN
858          topo_no_distinct = .TRUE.
859       ENDIF
860
861    END SUBROUTINE netcdf_data_input_inquire_file
862
863!------------------------------------------------------------------------------!
864! Description:
865! ------------
866!> Reads global attributes and coordinate reference system required for
867!> initialization of the model.
868!------------------------------------------------------------------------------!
869    SUBROUTINE netcdf_data_input_init
870
871       IMPLICIT NONE
872
873       INTEGER(iwp) ::  id_mod     !< NetCDF id of input file
874       INTEGER(iwp) ::  var_id_crs !< NetCDF id of variable crs
875
876!
877!--    Define default list of crs attributes. This is required for coordinate
878!--    transformation.
879       crs_list = (/ coord_ref_sys%semi_major_axis,                            &
880                     coord_ref_sys%inverse_flattening,                         &
881                     coord_ref_sys%longitude_of_prime_meridian,                &
882                     coord_ref_sys%longitude_of_central_meridian,              &
883                     coord_ref_sys%scale_factor_at_central_meridian,           &
884                     coord_ref_sys%latitude_of_projection_origin,              &
885                     coord_ref_sys%false_easting,                              &
886                     coord_ref_sys%false_northing /)
887
888       IF ( .NOT. input_pids_static )  RETURN
889
890#if defined ( __netcdf )
891!
892!--    Open file in read-only mode
893       CALL open_read_file( TRIM( input_file_static ) //                       &
894                            TRIM( coupling_char ), id_mod )
895!
896!--    Read global attributes
897       CALL get_attribute( id_mod, input_file_atts%origin_lat_char,            &
898                           input_file_atts%origin_lat, .TRUE. )
899
900       CALL get_attribute( id_mod, input_file_atts%origin_lon_char,            &
901                           input_file_atts%origin_lon, .TRUE. )
902
903       CALL get_attribute( id_mod, input_file_atts%origin_time_char,           &
904                           input_file_atts%origin_time, .TRUE. )
905
906       CALL get_attribute( id_mod, input_file_atts%origin_x_char,              &
907                           input_file_atts%origin_x, .TRUE. )
908
909       CALL get_attribute( id_mod, input_file_atts%origin_y_char,              &
910                           input_file_atts%origin_y, .TRUE. )
911
912       CALL get_attribute( id_mod, input_file_atts%origin_z_char,              &
913                           input_file_atts%origin_z, .TRUE. )
914
915       CALL get_attribute( id_mod, input_file_atts%rotation_angle_char,        &
916                           input_file_atts%rotation_angle, .TRUE. )
917
918       CALL get_attribute( id_mod, input_file_atts%author_char,                &
919                           input_file_atts%author, .TRUE., no_abort=.FALSE. )
920       CALL get_attribute( id_mod, input_file_atts%contact_person_char,        &
921                           input_file_atts%contact_person, .TRUE., no_abort=.FALSE. )
922       CALL get_attribute( id_mod, input_file_atts%institution_char,           &
923                           input_file_atts%institution,    .TRUE., no_abort=.FALSE. )
924       CALL get_attribute( id_mod, input_file_atts%acronym_char,               &
925                           input_file_atts%acronym,        .TRUE., no_abort=.FALSE. )
926
927       CALL get_attribute( id_mod, input_file_atts%campaign_char,              &
928                           input_file_atts%campaign, .TRUE., no_abort=.FALSE. )
929       CALL get_attribute( id_mod, input_file_atts%location_char,              &
930                           input_file_atts%location, .TRUE., no_abort=.FALSE. )
931       CALL get_attribute( id_mod, input_file_atts%site_char,                  &
932                           input_file_atts%site,     .TRUE., no_abort=.FALSE. )
933
934       CALL get_attribute( id_mod, input_file_atts%source_char,                &
935                           input_file_atts%source,     .TRUE., no_abort=.FALSE. )
936       CALL get_attribute( id_mod, input_file_atts%references_char,            &
937                           input_file_atts%references, .TRUE., no_abort=.FALSE. )
938       CALL get_attribute( id_mod, input_file_atts%keywords_char,              &
939                           input_file_atts%keywords,   .TRUE., no_abort=.FALSE. )
940       CALL get_attribute( id_mod, input_file_atts%licence_char,               &
941                           input_file_atts%licence,    .TRUE., no_abort=.FALSE. )
942       CALL get_attribute( id_mod, input_file_atts%comment_char,               &
943                           input_file_atts%comment,    .TRUE., no_abort=.FALSE. )
944!
945!--    Read coordinate reference system if available
946       nc_stat = NF90_INQ_VARID( id_mod, 'crs', var_id_crs )
947       IF ( nc_stat == NF90_NOERR )  THEN
948          CALL get_attribute( id_mod, 'epsg_code',                             &
949                              coord_ref_sys%epsg_code,                         &
950                              .FALSE., 'crs' )
951          CALL get_attribute( id_mod, 'false_easting',                         &
952                              coord_ref_sys%false_easting,                     &
953                              .FALSE., 'crs' )
954          CALL get_attribute( id_mod, 'false_northing',                        &
955                              coord_ref_sys%false_northing,                    &
956                              .FALSE., 'crs' )
957          CALL get_attribute( id_mod, 'grid_mapping_name',                     &
958                              coord_ref_sys%grid_mapping_name,                 &
959                              .FALSE., 'crs' )
960          CALL get_attribute( id_mod, 'inverse_flattening',                    &
961                              coord_ref_sys%inverse_flattening,                &
962                              .FALSE., 'crs' )
963          CALL get_attribute( id_mod, 'latitude_of_projection_origin',         &
964                              coord_ref_sys%latitude_of_projection_origin,     &
965                              .FALSE., 'crs' )
966          CALL get_attribute( id_mod, 'long_name',                             &
967                              coord_ref_sys%long_name,                         &
968                              .FALSE., 'crs' )
969          CALL get_attribute( id_mod, 'longitude_of_central_meridian',         &
970                              coord_ref_sys%longitude_of_central_meridian,     &
971                              .FALSE., 'crs' )
972          CALL get_attribute( id_mod, 'longitude_of_prime_meridian',           &
973                              coord_ref_sys%longitude_of_prime_meridian,       &
974                              .FALSE., 'crs' )
975          CALL get_attribute( id_mod, 'scale_factor_at_central_meridian',      &
976                              coord_ref_sys%scale_factor_at_central_meridian,  &
977                              .FALSE., 'crs' )
978          CALL get_attribute( id_mod, 'semi_major_axis',                       &
979                              coord_ref_sys%semi_major_axis,                   &
980                              .FALSE., 'crs' )
981          CALL get_attribute( id_mod, 'units',                                 &
982                              coord_ref_sys%units,                             &
983                              .FALSE., 'crs' )
984       ELSE
985!
986!--       Calculate central meridian from origin_lon
987          coord_ref_sys%longitude_of_central_meridian = &
988             CEILING( input_file_atts%origin_lon / 6.0_wp ) * 6.0_wp - 3.0_wp
989       ENDIF
990!
991!--    Finally, close input file
992       CALL close_input_file( id_mod )
993#endif
994!
995!--    Copy latitude, longitude, origin_z, rotation angle on init type
996!--    NOTE: A shifting height might have already been saved to orgin_z in
997!--    init_grid; therefore, do not override but add the reference height from
998!--    the input file.
999       init_model%latitude        = input_file_atts%origin_lat
1000       init_model%longitude       = input_file_atts%origin_lon
1001       init_model%origin_time     = input_file_atts%origin_time 
1002       init_model%origin_x        = input_file_atts%origin_x
1003       init_model%origin_y        = input_file_atts%origin_y
1004       init_model%origin_z        = init_model%origin_z + input_file_atts%origin_z
1005       init_model%rotation_angle  = input_file_atts%rotation_angle 
1006
1007!
1008!--    Update list of crs attributes. This is required for coordinate
1009!--    transformation.
1010       crs_list = (/ coord_ref_sys%semi_major_axis,                            &
1011                     coord_ref_sys%inverse_flattening,                         &
1012                     coord_ref_sys%longitude_of_prime_meridian,                &
1013                     coord_ref_sys%longitude_of_central_meridian,              &
1014                     coord_ref_sys%scale_factor_at_central_meridian,           &
1015                     coord_ref_sys%latitude_of_projection_origin,              &
1016                     coord_ref_sys%false_easting,                              &
1017                     coord_ref_sys%false_northing /)
1018!
1019!--    In case of nested runs, each model domain might have different longitude
1020!--    and latitude, which would result in different Coriolis parameters and
1021!--    sun-zenith angles. To avoid this, longitude and latitude in each model
1022!--    domain will be set to the values of the root model. Please note, this
1023!--    synchronization is required already here.
1024#if defined( __parallel )
1025       CALL MPI_BCAST( init_model%latitude,  1, MPI_REAL, 0,                   &
1026                       MPI_COMM_WORLD, ierr )
1027       CALL MPI_BCAST( init_model%longitude, 1, MPI_REAL, 0,                   &
1028                       MPI_COMM_WORLD, ierr )
1029#endif
1030
1031    END SUBROUTINE netcdf_data_input_init
1032
1033
1034!------------------------------------------------------------------------------!
1035! Description:
1036! ------------
1037!> Reads Chemistry NETCDF Input data, such as emission values, emission species, etc.
1038!------------------------------------------------------------------------------!
1039
1040    SUBROUTINE netcdf_data_input_chemistry_data(emt_att,emt)
1041
1042       USE chem_modules,                                       &
1043           ONLY:  emiss_lod, time_fac_type, surface_csflux_name
1044
1045       USE control_parameters,                                 &
1046           ONLY:  message_string
1047
1048       USE indices,                                            &
1049           ONLY:  nxl, nxr, nys, nyn
1050
1051       IMPLICIT NONE
1052
1053       TYPE(chem_emis_att_type), INTENT(INOUT)                             ::  emt_att
1054       TYPE(chem_emis_val_type), ALLOCATABLE, DIMENSION(:), INTENT(INOUT)  ::  emt
1055   
1056       INTEGER(iwp)  ::  i, j, k      !< generic counters
1057       INTEGER(iwp)  ::  ispec        !< index for number of emission species in input
1058       INTEGER(iwp)  ::  len_dims     !< Length of dimension
1059       INTEGER(iwp)  ::  num_vars     !< number of variables in netcdf input file
1060
1061!
1062!-- dum_var_4d are designed to read in emission_values from the chemistry netCDF file.
1063!-- Currently the vestigial "z" dimension in emission_values makes it a 5D array,
1064!-- hence the corresponding dum_var_5d array.  When the "z" dimension is removed
1065!-- completely, dum_var_4d will be used instead
1066!-- (ecc 20190425)
1067
1068!       REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:)    ::  dum_var_4d  !< temp array 4 4D chem emission data
1069       REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:,:)  ::  dum_var_5d  !< temp array 4 5D chem emission data
1070
1071!
1072!-- Start processing data
1073!
1074!-- Emission LOD 0 (Parameterized mode)
1075
1076        IF  ( emiss_lod == 0 )  THEN
1077
1078! for reference (ecc)
1079!       IF (TRIM(mode_emis) == "PARAMETERIZED" .OR. TRIM(mode_emis) == "parameterized") THEN
1080
1081           ispec=1
1082           emt_att%n_emiss_species = 0
1083
1084!
1085!-- number of species
1086
1087           DO  WHILE (TRIM( surface_csflux_name( ispec ) ) /= 'novalue' )
1088
1089             emt_att%n_emiss_species = emt_att%n_emiss_species + 1
1090             ispec=ispec+1
1091!
1092!-- followling line retained for compatibility with salsa_mod
1093!-- which still uses emt_att%nspec heavily (ecc)
1094
1095             emt_att%nspec = emt_att%nspec + 1
1096
1097           ENDDO
1098
1099!
1100!-- allocate emission values data type arrays
1101
1102          ALLOCATE ( emt(emt_att%n_emiss_species) )
1103
1104!
1105!-- Read EMISSION SPECIES NAMES
1106
1107!
1108!-- allocate space for strings
1109
1110          ALLOCATE (emt_att%species_name(emt_att%n_emiss_species) )
1111 
1112         DO ispec = 1, emt_att%n_emiss_species
1113            emt_att%species_name(ispec) = TRIM(surface_csflux_name(ispec))
1114         ENDDO
1115
1116!
1117!-- LOD 1 (default mode) and LOD 2 (pre-processed mode)
1118
1119       ELSE
1120
1121#if defined ( __netcdf )
1122
1123          IF ( .NOT. input_pids_chem )  RETURN
1124
1125!
1126!-- first we allocate memory space for the emission species and then
1127!-- we differentiate between LOD 1 (default mode) and LOD 2 (pre-processed mode)
1128
1129!
1130!-- open emission data file ( {palmcase}_chemistry )
1131
1132          CALL open_read_file ( TRIM(input_file_chem) // TRIM(coupling_char), id_emis )
1133
1134!
1135!-- inquire number of variables
1136
1137          CALL inquire_num_variables ( id_emis, num_vars )
1138
1139!
1140!-- Get General Dimension Lengths: only # species and # categories.
1141!-- Tther dimensions depend on the emission mode or specific components
1142
1143          CALL get_dimension_length ( id_emis, emt_att%n_emiss_species, 'nspecies' )
1144
1145!
1146!-- backward compatibility for salsa_mod (ecc)
1147
1148          emt_att%nspec = emt_att%n_emiss_species
1149
1150!
1151!-- Allocate emission values data type arrays
1152
1153          ALLOCATE ( emt(emt_att%n_emiss_species) )
1154
1155!
1156!-- READING IN SPECIES NAMES
1157
1158!
1159!-- Allocate memory for species names
1160
1161          ALLOCATE ( emt_att%species_name(emt_att%n_emiss_species) )
1162
1163!
1164!-- Retrieve variable name (again, should use n_emiss_strlen)
1165
1166          CALL get_variable( id_emis, 'emission_name',    &
1167                             string_values, emt_att%n_emiss_species )
1168          emt_att%species_name=string_values
1169
1170!
1171!-- dealocate string_values previously allocated in get_variable call
1172
1173          IF  ( ALLOCATED(string_values) )  DEALLOCATE (string_values)
1174
1175!
1176!-- READING IN SPECIES INDICES
1177
1178!
1179!-- Allocate memory for species indices
1180
1181          ALLOCATE ( emt_att%species_index(emt_att%n_emiss_species) )
1182
1183!
1184!-- Retrieve variable data
1185
1186          CALL get_variable( id_emis, 'emission_index', emt_att%species_index )
1187!
1188!-- Now the routine has to distinguish between chemistry emission
1189!-- LOD 1 (DEFAULT mode) and LOD 2 (PRE-PROCESSED mode)
1190
1191!
1192!-- START OF EMISSION LOD 1 (DEFAULT MODE)
1193
1194
1195          IF  ( emiss_lod == 1 )  THEN
1196
1197! for reference (ecc)
1198!          IF (TRIM(mode_emis) == "DEFAULT" .OR. TRIM(mode_emis) == "default") THEN
1199
1200!
1201!-- get number of emission categories
1202
1203             CALL get_dimension_length ( id_emis, emt_att%ncat, 'ncat' )
1204
1205!-- READING IN EMISSION CATEGORIES INDICES
1206
1207             ALLOCATE ( emt_att%cat_index(emt_att%ncat) )
1208
1209!
1210!-- Retrieve variable data
1211
1212             CALL get_variable( id_emis, 'emission_cat_index', emt_att%cat_index )
1213
1214
1215!
1216!-- Loop through individual species to get basic information on
1217!-- VOC/PM/NOX/SOX
1218
1219!------------------------------------------------------------------------------
1220!-- NOTE - CHECK ARRAY INDICES FOR READING IN NAMES AND SPECIES
1221!--        IN LOD1 (DEFAULT MODE) FOR THE VARIOUS MODE SPLITS
1222!--        AS ALL ID_EMIS CONDITIONALS HAVE BEEN REMOVED FROM GET_VAR
1223!--        FUNCTIONS.  IN THEORY THIS WOULD MEAN ALL ARRAYS SHOULD BE
1224!--        READ FROM 0 to N-1 (C CONVENTION) AS OPPOSED TO 1 to N
1225!--        (FORTRAN CONVENTION).  KEEP THIS IN MIND !!
1226!--        (ecc 20190424)
1227!------------------------------------------------------------------------------
1228 
1229             DO  ispec = 1, emt_att%n_emiss_species
1230
1231!
1232!-- VOC DATA (name and composition)
1233
1234                IF  ( TRIM(emt_att%species_name(ispec)) == "VOC" .OR.                  &
1235                      TRIM(emt_att%species_name(ispec)) == "voc" )  THEN
1236
1237!
1238!-- VOC name
1239                   CALL get_dimension_length ( id_emis, emt_att%nvoc, 'nvoc' )
1240                   ALLOCATE ( emt_att%voc_name(emt_att%nvoc) )
1241                   CALL get_variable ( id_emis,"emission_voc_name",  &
1242                                       string_values, emt_att%nvoc )
1243                   emt_att%voc_name = string_values
1244                   IF  ( ALLOCATED(string_values) )  DEALLOCATE (string_values)
1245
1246!
1247!-- VOC composition
1248
1249                   ALLOCATE ( emt_att%voc_comp(emt_att%ncat,emt_att%nvoc) )
1250                   CALL get_variable ( id_emis, "composition_voc", emt_att%voc_comp,     &
1251                                       1, emt_att%ncat, 1, emt_att%nvoc )
1252
1253                ENDIF  ! VOC
1254
1255!
1256!-- PM DATA (name and composition)
1257
1258                IF  ( TRIM(emt_att%species_name(ispec)) == "PM" .OR.                   &
1259                      TRIM(emt_att%species_name(ispec)) == "pm")  THEN
1260
1261!
1262!-- PM name
1263
1264                   CALL get_dimension_length ( id_emis, emt_att%npm, 'npm' )
1265                   ALLOCATE ( emt_att%pm_name(emt_att%npm) )
1266                   CALL get_variable ( id_emis, "pm_name", string_values, emt_att%npm )
1267                   emt_att%pm_name = string_values
1268                   IF  ( ALLOCATED(string_values) )  DEALLOCATE (string_values)     
1269
1270!
1271!-- PM composition (PM1, PM2.5 and PM10)
1272
1273                   len_dims = 3  ! PM1, PM2.5, PM10
1274                   ALLOCATE(emt_att%pm_comp(emt_att%ncat,emt_att%npm,len_dims))
1275                   CALL get_variable ( id_emis, "composition_pm", emt_att%pm_comp,       &
1276                                       1, emt_att%ncat, 1, emt_att%npm, 1, len_dims )
1277
1278                ENDIF  ! PM
1279
1280!
1281!-- NOX (NO and NO2)
1282
1283                IF  ( TRIM(emt_att%species_name(ispec)) == "NOX" .OR.                  &
1284                      TRIM(emt_att%species_name(ispec)) == "nox" )  THEN
1285
1286                   ALLOCATE ( emt_att%nox_comp(emt_att%ncat,emt_att%nnox) )
1287                   CALL get_variable ( id_emis, "composition_nox", emt_att%nox_comp,     &
1288                                       1, emt_att%ncat, 1, emt_att%nnox )
1289
1290                ENDIF  ! NOX
1291
1292!
1293!-- SOX (SO2 and SO4)
1294
1295                IF  ( TRIM(emt_att%species_name(ispec)) == "SOX" .OR.                  &
1296                      TRIM(emt_att%species_name(ispec)) == "sox" )  THEN
1297
1298                   ALLOCATE ( emt_att%sox_comp(emt_att%ncat,emt_att%nsox) )
1299                   CALL get_variable ( id_emis, "composition_sox", emt_att%sox_comp,     &
1300                                       1, emt_att%ncat, 1, emt_att%nsox )
1301
1302                ENDIF  ! SOX
1303
1304             ENDDO  ! do ispec
1305
1306!
1307!-- EMISSION TIME SCALING FACTORS (hourly and MDH data)
1308 
1309!     
1310!-- HOUR   
1311             IF  ( TRIM(time_fac_type) == "HOUR" .OR.                        &
1312                   TRIM(time_fac_type) == "hour" )  THEN
1313
1314                CALL get_dimension_length ( id_emis, emt_att%nhoursyear, 'nhoursyear' )
1315                ALLOCATE ( emt_att%hourly_emis_time_factor(emt_att%ncat,emt_att%nhoursyear) )
1316                CALL get_variable ( id_emis, "emission_time_factors",          &
1317                                    emt_att%hourly_emis_time_factor,           &
1318                                    1, emt_att%ncat, 1, emt_att%nhoursyear )
1319
1320!
1321!-- MDH
1322
1323             ELSE IF  ( TRIM(time_fac_type)  ==  "MDH" .OR.                  &
1324                        TRIM(time_fac_type)  ==  "mdh" )  THEN
1325
1326                CALL get_dimension_length ( id_emis, emt_att%nmonthdayhour, 'nmonthdayhour' )
1327                ALLOCATE ( emt_att%mdh_emis_time_factor(emt_att%ncat,emt_att%nmonthdayhour) )
1328                CALL get_variable ( id_emis, "emission_time_factors",          &
1329                                    emt_att%mdh_emis_time_factor,              &
1330                                    1, emt_att%ncat, 1, emt_att%nmonthdayhour )
1331
1332!
1333!-- ERROR (time factor undefined)
1334
1335             ELSE
1336
1337                message_string = 'We are in the DEFAULT chemistry emissions mode: '  //  &
1338                                 '     !no time-factor type specified!'              //  &
1339                                 'Please specify the value of time_fac_type:'        //  &
1340                                 '         either "MDH" or "HOUR"'                 
1341                CALL message( 'netcdf_data_input_chemistry_data', 'CM0200', 2, 2, 0, 6, 0 ) 
1342 
1343
1344             ENDIF  ! time_fac_type
1345
1346!
1347!-- read in default (LOD1) emissions from chemisty netCDF file per species
1348
1349!
1350!-- NOTE - at the moment the data is read in per species, but in the future it would
1351!--        be much more sensible to read in per species per time step to reduce
1352!--        memory consumption and, to a lesser degree, dimensionality of data exchange
1353!--        (I expect this will be necessary when the problem size is large)
1354
1355             DO ispec = 1, emt_att%n_emiss_species
1356
1357!
1358!-- allocate space for species specific emission values
1359!-- NOTE - this array is extended by 1 cell in each horizontal direction
1360!--        to compensate for an apparent linear offset.  The reason of this
1361!--        offset is not known but it has been determined to take place beyond the
1362!--        scope of this module, and has little to do with index conventions.
1363!--        That is, setting the array horizontal limit from nx0:nx1 to 1:(nx1-nx0+1)
1364!--        or nx0+1:nx1+1 did not result in correct or definite behavior
1365!--        This must be looked at at some point by the Hannover team but for now
1366!--        this workaround is deemed reasonable (ecc 20190417)
1367
1368                IF ( .NOT. ALLOCATED ( emt(ispec)%default_emission_data ) )  THEN
1369                    ALLOCATE ( emt(ispec)%default_emission_data(emt_att%ncat,nys:nyn+1,nxl:nxr+1) )
1370                ENDIF
1371!
1372!-- allocate dummy variable w/ index order identical to that shown in the netCDF header
1373
1374                ALLOCATE ( dum_var_5d(1,nys:nyn,nxl:nxr,1,emt_att%ncat) )
1375!
1376!-- get variable.  be very careful
1377!-- I am using get_variable_5d_real_dynamic (note logical argument at the end)
1378!-- 1) use Fortran index convention (i.e., 1 to N)
1379!-- 2) index order must be in reverse order from above allocation order
1380 
1381                CALL get_variable ( id_emis, "emission_values", dum_var_5d, &
1382                                    1,            ispec, nxl+1,     nys+1,     1,                    &
1383                                    emt_att%ncat, 1,     nxr-nxl+1, nyn-nys+1, emt_att%dt_emission,  &
1384                                    .FALSE. )
1385!
1386!-- assign temp array to data structure then deallocate temp array
1387!-- NOTE - indices are shifted from nx0:nx1 to nx0+1:nx1+1 to offset
1388!--        the emission data array to counter said domain offset
1389!--        (ecc 20190417)
1390
1391                DO k = 1, emt_att%ncat
1392                   DO j = nys+1, nyn+1
1393                      DO i = nxl+1, nxr+1
1394                         emt(ispec)%default_emission_data(k,j,i) = dum_var_5d(1,j-1,i-1,1,k)
1395                      ENDDO
1396                   ENDDO
1397                ENDDO
1398
1399                DEALLOCATE ( dum_var_5d )
1400
1401             ENDDO  ! ispec
1402!
1403!-- UNITS
1404
1405             CALL get_attribute(id_emis,"units",emt_att%units,.FALSE.,"emission_values")
1406
1407!
1408!-- END DEFAULT MODE
1409
1410
1411!
1412!-- START LOD 2 (PRE-PROCESSED MODE)
1413
1414          ELSE IF  ( emiss_lod == 2 )  THEN
1415
1416! for reference (ecc)
1417!          ELSE IF (TRIM(mode_emis) == "PRE-PROCESSED" .OR. TRIM(mode_emis) == "pre-processed") THEN
1418
1419!
1420!-- For LOD 2 only VOC and emission data need be read
1421
1422!------------------------------------------------------------------------------
1423!-- NOTE - CHECK ARRAY INDICES FOR READING IN NAMES AND SPECIES
1424!--        IN LOD2 (PRE-PROCESSED MODE) FOR THE VARIOUS MODE SPLITS
1425!--        AS ALL ID_EMIS CONDITIONALS HAVE BEEN REMOVED FROM GET_VAR
1426!--        FUNCTIONS.  IN THEORY THIS WOULD MEAN ALL ARRAYS SHOULD BE
1427!--        READ FROM 0 to N-1 (C CONVENTION) AS OPPOSED TO 1 to N
1428!--        (FORTRAN CONVENTION).  KEEP THIS IN MIND !!
1429!--        (ecc 20190424)
1430!------------------------------------------------------------------------------
1431
1432             DO ispec = 1, emt_att%n_emiss_species
1433
1434!
1435!-- VOC DATA (name and composition)
1436
1437                IF  ( TRIM(emt_att%species_name(ispec)) == "VOC" .OR.                  &
1438                      TRIM(emt_att%species_name(ispec)) == "voc" )  THEN
1439
1440!
1441!-- VOC name
1442                   CALL get_dimension_length ( id_emis, emt_att%nvoc, 'nvoc' )
1443                   ALLOCATE ( emt_att%voc_name(emt_att%nvoc) )
1444                   CALL get_variable ( id_emis, "emission_voc_name",                     &
1445                                       string_values, emt_att%nvoc)
1446                   emt_att%voc_name = string_values
1447                   IF  ( ALLOCATED(string_values) )  DEALLOCATE (string_values)
1448
1449!
1450!-- VOC composition
1451 
1452                   ALLOCATE ( emt_att%voc_comp(emt_att%ncat,emt_att%nvoc) )
1453                   CALL get_variable ( id_emis, "composition_voc", emt_att%voc_comp,     &
1454                                       1, emt_att%ncat, 1, emt_att%nvoc )
1455                ENDIF  ! VOC
1456 
1457             ENDDO  ! ispec
1458
1459!
1460!-- EMISSION DATA
1461
1462             CALL get_dimension_length ( id_emis, emt_att%dt_emission, 'time' )   
1463 
1464!
1465!-- read in pre-processed (LOD2) emissions from chemisty netCDF file per species
1466
1467!
1468!-- NOTE - at the moment the data is read in per species, but in the future it would
1469!--        be much more sensible to read in per species per time step to reduce
1470!--        memory consumption and, to a lesser degree, dimensionality of data exchange
1471!--        (I expect this will be necessary when the problem size is large)
1472
1473             DO ispec = 1, emt_att%n_emiss_species
1474
1475!
1476!-- allocate space for species specific emission values
1477!-- NOTE - this array is extended by 1 cell in each horizontal direction
1478!--        to compensate for an apparent linear offset.  The reason of this
1479!--        offset is not known but it has been determined to take place beyond the
1480!--        scope of this module, and has little to do with index conventions.
1481!--        That is, setting the array horizontal limit from nx0:nx1 to 1:(nx1-nx0+1)
1482!--        or nx0+1:nx1+1 did not result in correct or definite behavior
1483!--        This must be looked at at some point by the Hannover team but for now
1484!--        this workaround is deemed reasonable (ecc 20190417)
1485
1486                IF ( .NOT. ALLOCATED( emt(ispec)%preproc_emission_data ) )  THEN
1487                   ALLOCATE( emt(ispec)%preproc_emission_data(                           &
1488                             emt_att%dt_emission, 1, nys:nyn+1, nxl:nxr+1) )
1489                ENDIF
1490!
1491!-- allocate dummy variable w/ index order identical to that shown in the netCDF header
1492
1493                ALLOCATE ( dum_var_5d(emt_att%dt_emission,1,nys:nyn,nxl:nxr,1) )
1494!
1495!-- get variable.  be very careful
1496!-- I am using get_variable_5d_real_dynamic (note logical argument at the end)
1497!-- 1) use Fortran index convention (i.e., 1 to N)
1498!-- 2) index order must be in reverse order from above allocation order
1499
1500                CALL get_variable ( id_emis, "emission_values", dum_var_5d, &
1501                                    ispec, nxl+1,     nys+1,     1, 1,                   &
1502                                    1,     nxr-nxl+1, nyn-nys+1, 1, emt_att%dt_emission, &
1503                                    .FALSE. )
1504!
1505!-- assign temp array to data structure then deallocate temp array
1506!-- NOTE - indices are shifted from nx0:nx1 to nx0+1:nx1+1 to offset
1507!--        the emission data array to counter said unkonwn offset
1508!--        (ecc 20190417)
1509
1510                DO k = 1, emt_att%dt_emission
1511                   DO j = nys+1, nyn+1
1512                      DO i = nxl+1, nxr+1
1513                         emt(ispec)%preproc_emission_data(k,1,j,i) = dum_var_5d(k,1,j-1,i-1,1)
1514                      ENDDO
1515                   ENDDO
1516                ENDDO
1517
1518                DEALLOCATE ( dum_var_5d )
1519
1520             ENDDO  ! ispec
1521!
1522!-- UNITS
1523
1524             CALL get_attribute ( id_emis, "units", emt_att%units, .FALSE. , "emission_values" )
1525       
1526          ENDIF  ! LOD1 & LOD2 (default and pre-processed mode)
1527
1528          CALL close_input_file (id_emis)
1529
1530#endif
1531
1532       ENDIF ! LOD0 (parameterized mode)
1533
1534    END SUBROUTINE netcdf_data_input_chemistry_data
1535
1536
1537!------------------------------------------------------------------------------!
1538! Description:
1539! ------------
1540!> Reads surface classification data, such as vegetation and soil type, etc. .
1541!------------------------------------------------------------------------------!
1542    SUBROUTINE netcdf_data_input_surface_data
1543
1544       USE control_parameters,                                                 &
1545           ONLY:  land_surface, urban_surface
1546
1547       USE exchange_horiz_mod,                                                 &
1548           ONLY:  exchange_horiz_2d, exchange_horiz_2d_byte, exchange_horiz_2d_int
1549
1550       USE indices,                                                            &
1551           ONLY:  nbgp, nxl, nxr, nyn, nys
1552
1553
1554       IMPLICIT NONE
1555
1556       CHARACTER(LEN=100), DIMENSION(:), ALLOCATABLE ::  var_names  !< variable names in static input file
1557
1558       INTEGER(iwp) ::  id_surf   !< NetCDF id of input file
1559       INTEGER(iwp) ::  k         !< running index along z-direction
1560       INTEGER(iwp) ::  k2        !< running index
1561       INTEGER(iwp) ::  num_vars  !< number of variables in input file
1562       INTEGER(iwp) ::  nz_soil   !< number of soil layers in file
1563
1564!
1565!--    If not static input file is available, skip this routine
1566       IF ( .NOT. input_pids_static )  RETURN
1567!
1568!--    Measure CPU time
1569       CALL cpu_log( log_point_s(82), 'NetCDF input', 'start' )
1570!
1571!--    Skip the following if no land-surface or urban-surface module are
1572!--    applied. This case, no one of the following variables is used anyway.
1573       IF (  .NOT. land_surface  .AND.  .NOT. urban_surface )  RETURN
1574
1575#if defined ( __netcdf )
1576!
1577!--    Open file in read-only mode
1578       CALL open_read_file( TRIM( input_file_static ) //                       &
1579                            TRIM( coupling_char ) , id_surf )
1580!
1581!--    Inquire all variable names.
1582!--    This will be used to check whether an optional input variable exist
1583!--    or not.
1584       CALL inquire_num_variables( id_surf, num_vars )
1585
1586       ALLOCATE( var_names(1:num_vars) )
1587       CALL inquire_variable_names( id_surf, var_names )
1588!
1589!--    Read vegetation type and required attributes
1590       IF ( check_existence( var_names, 'vegetation_type' ) )  THEN
1591          vegetation_type_f%from_file = .TRUE.
1592          CALL get_attribute( id_surf, char_fill,                              &
1593                              vegetation_type_f%fill,                          &
1594                              .FALSE., 'vegetation_type' )
1595
1596          ALLOCATE ( vegetation_type_f%var(nys:nyn,nxl:nxr)  )
1597
1598          CALL get_variable( id_surf, 'vegetation_type',                       &
1599                             vegetation_type_f%var, nxl, nxr, nys, nyn )
1600       ELSE
1601          vegetation_type_f%from_file = .FALSE.
1602       ENDIF
1603
1604!
1605!--    Read soil type and required attributes
1606       IF ( check_existence( var_names, 'soil_type' ) )  THEN
1607          soil_type_f%from_file = .TRUE.
1608!
1609!--       Note, lod is currently not on file; skip for the moment
1610!           CALL get_attribute( id_surf, char_lod,                       &
1611!                                      soil_type_f%lod,                  &
1612!                                      .FALSE., 'soil_type' )
1613          CALL get_attribute( id_surf, char_fill,                              &
1614                              soil_type_f%fill,                                &
1615                              .FALSE., 'soil_type' )
1616
1617          IF ( soil_type_f%lod == 1 )  THEN
1618
1619             ALLOCATE ( soil_type_f%var_2d(nys:nyn,nxl:nxr)  )
1620
1621             CALL get_variable( id_surf, 'soil_type', soil_type_f%var_2d,      &
1622                                nxl, nxr, nys, nyn )
1623
1624          ELSEIF ( soil_type_f%lod == 2 )  THEN
1625!
1626!--          Obtain number of soil layers from file.
1627             CALL get_dimension_length( id_surf, nz_soil, 'zsoil' )
1628
1629             ALLOCATE ( soil_type_f%var_3d(0:nz_soil,nys:nyn,nxl:nxr) )
1630
1631             CALL get_variable( id_surf, 'soil_type', soil_type_f%var_3d,      &
1632                                nxl, nxr, nys, nyn, 0, nz_soil )
1633 
1634          ENDIF
1635       ELSE
1636          soil_type_f%from_file = .FALSE.
1637       ENDIF
1638
1639!
1640!--    Read pavement type and required attributes
1641       IF ( check_existence( var_names, 'pavement_type' ) )  THEN
1642          pavement_type_f%from_file = .TRUE.
1643          CALL get_attribute( id_surf, char_fill,                              &
1644                              pavement_type_f%fill, .FALSE.,                   &
1645                              'pavement_type' )
1646
1647          ALLOCATE ( pavement_type_f%var(nys:nyn,nxl:nxr)  )
1648
1649          CALL get_variable( id_surf, 'pavement_type', pavement_type_f%var,    &
1650                             nxl, nxr, nys, nyn )
1651       ELSE
1652          pavement_type_f%from_file = .FALSE.
1653       ENDIF
1654
1655!
1656!--    Read water type and required attributes
1657       IF ( check_existence( var_names, 'water_type' ) )  THEN
1658          water_type_f%from_file = .TRUE.
1659          CALL get_attribute( id_surf, char_fill, water_type_f%fill,           &
1660                              .FALSE., 'water_type' )
1661
1662          ALLOCATE ( water_type_f%var(nys:nyn,nxl:nxr)  )
1663
1664          CALL get_variable( id_surf, 'water_type', water_type_f%var,          &
1665                             nxl, nxr, nys, nyn )
1666
1667       ELSE
1668          water_type_f%from_file = .FALSE.
1669       ENDIF
1670!
1671!--    Read relative surface fractions of vegetation, pavement and water.
1672       IF ( check_existence( var_names, 'surface_fraction' ) )  THEN
1673          surface_fraction_f%from_file = .TRUE.
1674          CALL get_attribute( id_surf, char_fill,                              &
1675                              surface_fraction_f%fill,                         &
1676                              .FALSE., 'surface_fraction' )
1677!
1678!--       Inquire number of surface fractions
1679          CALL get_dimension_length( id_surf,                                  &
1680                                     surface_fraction_f%nf,                    &
1681                                     'nsurface_fraction' )
1682!
1683!--       Allocate dimension array and input array for surface fractions
1684          ALLOCATE( surface_fraction_f%nfracs(0:surface_fraction_f%nf-1) )
1685          ALLOCATE( surface_fraction_f%frac(0:surface_fraction_f%nf-1,         &
1686                                            nys:nyn,nxl:nxr) )
1687!
1688!--       Get dimension of surface fractions
1689          CALL get_variable( id_surf, 'nsurface_fraction',                     &
1690                             surface_fraction_f%nfracs )
1691!
1692!--       Read surface fractions
1693          CALL get_variable( id_surf, 'surface_fraction',                      &
1694                             surface_fraction_f%frac, nxl, nxr, nys, nyn,      &
1695                             0, surface_fraction_f%nf-1 )
1696       ELSE
1697          surface_fraction_f%from_file = .FALSE.
1698       ENDIF
1699!
1700!--    Read building parameters and related information
1701       IF ( check_existence( var_names, 'building_pars' ) )  THEN
1702          building_pars_f%from_file = .TRUE.
1703          CALL get_attribute( id_surf, char_fill,                              &
1704                              building_pars_f%fill,                            &
1705                              .FALSE., 'building_pars' )
1706!
1707!--       Inquire number of building parameters
1708          CALL get_dimension_length( id_surf,                                  &
1709                                      building_pars_f%np,                      &
1710                                      'nbuilding_pars' )
1711!
1712!--       Allocate dimension array and input array for building parameters
1713          ALLOCATE( building_pars_f%pars(0:building_pars_f%np-1) )
1714          ALLOCATE( building_pars_f%pars_xy(0:building_pars_f%np-1,            &
1715                                            nys:nyn,nxl:nxr) )
1716!
1717!--       Get dimension of building parameters
1718          CALL get_variable( id_surf, 'nbuilding_pars',                        &
1719                             building_pars_f%pars )
1720!
1721!--       Read building_pars
1722          CALL get_variable( id_surf, 'building_pars',                         &
1723                             building_pars_f%pars_xy, nxl, nxr, nys, nyn,      &
1724                             0, building_pars_f%np-1 )
1725       ELSE
1726          building_pars_f%from_file = .FALSE.
1727       ENDIF
1728!
1729!--    Read building surface parameters
1730       IF ( check_existence( var_names, 'building_surface_pars' ) )  THEN
1731          building_surface_pars_f%from_file = .TRUE.
1732          CALL get_attribute( id_surf, char_fill,                              &
1733                              building_surface_pars_f%fill,                    &
1734                              .FALSE., 'building_surface_pars' )
1735!
1736!--       Read building_surface_pars
1737          CALL get_variable_surf( id_surf, 'building_surface_pars', &
1738                                  building_surface_pars_f )
1739       ELSE
1740          building_surface_pars_f%from_file = .FALSE.
1741       ENDIF
1742
1743!
1744!--    Read albedo type and required attributes
1745       IF ( check_existence( var_names, 'albedo_type' ) )  THEN
1746          albedo_type_f%from_file = .TRUE.
1747          CALL get_attribute( id_surf, char_fill, albedo_type_f%fill,          &
1748                              .FALSE.,  'albedo_type' )
1749
1750          ALLOCATE ( albedo_type_f%var(nys:nyn,nxl:nxr)  )
1751         
1752          CALL get_variable( id_surf, 'albedo_type', albedo_type_f%var,        &
1753                             nxl, nxr, nys, nyn )
1754       ELSE
1755          albedo_type_f%from_file = .FALSE.
1756       ENDIF
1757!
1758!--    Read albedo parameters and related information
1759       IF ( check_existence( var_names, 'albedo_pars' ) )  THEN
1760          albedo_pars_f%from_file = .TRUE.
1761          CALL get_attribute( id_surf, char_fill, albedo_pars_f%fill,          &
1762                              .FALSE., 'albedo_pars' )
1763!
1764!--       Inquire number of albedo parameters
1765          CALL get_dimension_length( id_surf,                                  &
1766                                     albedo_pars_f%np,                         &
1767                                     'nalbedo_pars' )
1768!
1769!--       Allocate dimension array and input array for albedo parameters
1770          ALLOCATE( albedo_pars_f%pars(0:albedo_pars_f%np-1) )
1771          ALLOCATE( albedo_pars_f%pars_xy(0:albedo_pars_f%np-1,                &
1772                                          nys:nyn,nxl:nxr) )
1773!
1774!--       Get dimension of albedo parameters
1775          CALL get_variable( id_surf, 'nalbedo_pars', albedo_pars_f%pars )
1776
1777          CALL get_variable( id_surf, 'albedo_pars', albedo_pars_f%pars_xy,    &
1778                             nxl, nxr, nys, nyn,                               &
1779                             0, albedo_pars_f%np-1 )
1780       ELSE
1781          albedo_pars_f%from_file = .FALSE.
1782       ENDIF
1783
1784!
1785!--    Read pavement parameters and related information
1786       IF ( check_existence( var_names, 'pavement_pars' ) )  THEN
1787          pavement_pars_f%from_file = .TRUE.
1788          CALL get_attribute( id_surf, char_fill,                              &
1789                              pavement_pars_f%fill,                            &
1790                              .FALSE., 'pavement_pars' )
1791!
1792!--       Inquire number of pavement parameters
1793          CALL get_dimension_length( id_surf,                                  &
1794                                     pavement_pars_f%np,                       &
1795                                     'npavement_pars' )
1796!
1797!--       Allocate dimension array and input array for pavement parameters
1798          ALLOCATE( pavement_pars_f%pars(0:pavement_pars_f%np-1) )
1799          ALLOCATE( pavement_pars_f%pars_xy(0:pavement_pars_f%np-1,            &
1800                                            nys:nyn,nxl:nxr) )
1801!
1802!--       Get dimension of pavement parameters
1803          CALL get_variable( id_surf, 'npavement_pars', pavement_pars_f%pars )
1804
1805          CALL get_variable( id_surf, 'pavement_pars', pavement_pars_f%pars_xy,&
1806                             nxl, nxr, nys, nyn,                               &
1807                             0, pavement_pars_f%np-1 )
1808       ELSE
1809          pavement_pars_f%from_file = .FALSE.
1810       ENDIF
1811
1812!
1813!--    Read pavement subsurface parameters and related information
1814       IF ( check_existence( var_names, 'pavement_subsurface_pars' ) )         &
1815       THEN
1816          pavement_subsurface_pars_f%from_file = .TRUE.
1817          CALL get_attribute( id_surf, char_fill,                              &
1818                              pavement_subsurface_pars_f%fill,                 &
1819                              .FALSE., 'pavement_subsurface_pars' )
1820!
1821!--       Inquire number of parameters
1822          CALL get_dimension_length( id_surf,                                  &
1823                                     pavement_subsurface_pars_f%np,            &
1824                                     'npavement_subsurface_pars' )
1825!
1826!--       Inquire number of soil layers
1827          CALL get_dimension_length( id_surf,                                  &
1828                                     pavement_subsurface_pars_f%nz,            &
1829                                     'zsoil' )
1830!
1831!--       Allocate dimension array and input array for pavement parameters
1832          ALLOCATE( pavement_subsurface_pars_f%pars                            &
1833                            (0:pavement_subsurface_pars_f%np-1) )
1834          ALLOCATE( pavement_subsurface_pars_f%pars_xyz                        &
1835                            (0:pavement_subsurface_pars_f%np-1,                &
1836                             0:pavement_subsurface_pars_f%nz-1,                &
1837                             nys:nyn,nxl:nxr) )
1838!
1839!--       Get dimension of pavement parameters
1840          CALL get_variable( id_surf, 'npavement_subsurface_pars',             &
1841                             pavement_subsurface_pars_f%pars )
1842
1843          CALL get_variable( id_surf, 'pavement_subsurface_pars',              &
1844                             pavement_subsurface_pars_f%pars_xyz,              &
1845                             nxl, nxr, nys, nyn,                               &
1846                             0, pavement_subsurface_pars_f%nz-1,               &
1847                             0, pavement_subsurface_pars_f%np-1 )
1848       ELSE
1849          pavement_subsurface_pars_f%from_file = .FALSE.
1850       ENDIF
1851
1852
1853!
1854!--    Read vegetation parameters and related information
1855       IF ( check_existence( var_names, 'vegetation_pars' ) )  THEN
1856          vegetation_pars_f%from_file = .TRUE.
1857          CALL get_attribute( id_surf, char_fill,                              &
1858                              vegetation_pars_f%fill,                          &
1859                              .FALSE.,  'vegetation_pars' )
1860!
1861!--       Inquire number of vegetation parameters
1862          CALL get_dimension_length( id_surf,                                  &
1863                                     vegetation_pars_f%np,                     &
1864                                     'nvegetation_pars' )
1865!
1866!--       Allocate dimension array and input array for surface fractions
1867          ALLOCATE( vegetation_pars_f%pars(0:vegetation_pars_f%np-1) )
1868          ALLOCATE( vegetation_pars_f%pars_xy(0:vegetation_pars_f%np-1,        &
1869                                              nys:nyn,nxl:nxr) )
1870!
1871!--       Get dimension of the parameters
1872          CALL get_variable( id_surf, 'nvegetation_pars',                      &
1873                             vegetation_pars_f%pars )
1874
1875          CALL get_variable( id_surf, 'vegetation_pars',                       &
1876                             vegetation_pars_f%pars_xy, nxl, nxr, nys, nyn,    &
1877                             0, vegetation_pars_f%np-1 )
1878       ELSE
1879          vegetation_pars_f%from_file = .FALSE.
1880       ENDIF
1881
1882!
1883!--    Read root parameters/distribution and related information
1884       IF ( check_existence( var_names, 'soil_pars' ) )  THEN
1885          soil_pars_f%from_file = .TRUE.
1886          CALL get_attribute( id_surf, char_fill,                              &
1887                              soil_pars_f%fill,                                &
1888                              .FALSE., 'soil_pars' )
1889
1890          CALL get_attribute( id_surf, char_lod,                               &
1891                              soil_pars_f%lod,                                 &
1892                              .FALSE., 'soil_pars' )
1893
1894!
1895!--       Inquire number of soil parameters
1896          CALL get_dimension_length( id_surf,                                  &
1897                                     soil_pars_f%np,                           &
1898                                     'nsoil_pars' )
1899!
1900!--       Read parameters array
1901          ALLOCATE( soil_pars_f%pars(0:soil_pars_f%np-1) )
1902          CALL get_variable( id_surf, 'nsoil_pars', soil_pars_f%pars )
1903
1904!
1905!--       In case of level of detail 2, also inquire number of vertical
1906!--       soil layers, allocate memory and read the respective dimension
1907          IF ( soil_pars_f%lod == 2 )  THEN
1908             CALL get_dimension_length( id_surf,                               &
1909                                        soil_pars_f%nz,                        &
1910                                        'zsoil' )
1911
1912             ALLOCATE( soil_pars_f%layers(0:soil_pars_f%nz-1) )
1913             CALL get_variable( id_surf, 'zsoil', soil_pars_f%layers )
1914
1915          ENDIF
1916
1917!
1918!--       Read soil parameters, depending on level of detail
1919          IF ( soil_pars_f%lod == 1 )  THEN
1920             ALLOCATE( soil_pars_f%pars_xy(0:soil_pars_f%np-1,                 &
1921                                           nys:nyn,nxl:nxr) )
1922                 
1923             CALL get_variable( id_surf, 'soil_pars', soil_pars_f%pars_xy,     &
1924                                nxl, nxr, nys, nyn, 0, soil_pars_f%np-1 )
1925
1926          ELSEIF ( soil_pars_f%lod == 2 )  THEN
1927             ALLOCATE( soil_pars_f%pars_xyz(0:soil_pars_f%np-1,                &
1928                                            0:soil_pars_f%nz-1,                &
1929                                            nys:nyn,nxl:nxr) )
1930             CALL get_variable( id_surf, 'soil_pars',                          &
1931                                soil_pars_f%pars_xyz,                          &
1932                                nxl, nxr, nys, nyn, 0, soil_pars_f%nz-1,       &
1933                                0, soil_pars_f%np-1 )
1934
1935          ENDIF
1936       ELSE
1937          soil_pars_f%from_file = .FALSE.
1938       ENDIF
1939
1940!
1941!--    Read water parameters and related information
1942       IF ( check_existence( var_names, 'water_pars' ) )  THEN
1943          water_pars_f%from_file = .TRUE.
1944          CALL get_attribute( id_surf, char_fill,                              &
1945                              water_pars_f%fill,                               &
1946                              .FALSE., 'water_pars' )
1947!
1948!--       Inquire number of water parameters
1949          CALL get_dimension_length( id_surf,                                  &
1950                                     water_pars_f%np,                          &
1951                                     'nwater_pars' )
1952!
1953!--       Allocate dimension array and input array for water parameters
1954          ALLOCATE( water_pars_f%pars(0:water_pars_f%np-1) )
1955          ALLOCATE( water_pars_f%pars_xy(0:water_pars_f%np-1,                  &
1956                                         nys:nyn,nxl:nxr) )
1957!
1958!--       Get dimension of water parameters
1959          CALL get_variable( id_surf, 'nwater_pars', water_pars_f%pars )
1960
1961          CALL get_variable( id_surf, 'water_pars', water_pars_f%pars_xy,      &
1962                             nxl, nxr, nys, nyn, 0, water_pars_f%np-1 )
1963       ELSE
1964          water_pars_f%from_file = .FALSE.
1965       ENDIF
1966!
1967!--    Read root area density - parametrized vegetation
1968       IF ( check_existence( var_names, 'root_area_dens_s' ) )  THEN
1969          root_area_density_lsm_f%from_file = .TRUE.
1970          CALL get_attribute( id_surf, char_fill,                              &
1971                              root_area_density_lsm_f%fill,                    &
1972                              .FALSE., 'root_area_dens_s' )
1973!
1974!--       Obtain number of soil layers from file and allocate variable
1975          CALL get_dimension_length( id_surf,                                  &
1976                                     root_area_density_lsm_f%nz,               &
1977                                     'zsoil' )
1978          ALLOCATE( root_area_density_lsm_f%var                                &
1979                                        (0:root_area_density_lsm_f%nz-1,       &
1980                                         nys:nyn,nxl:nxr) )
1981
1982!
1983!--       Read root-area density
1984          CALL get_variable( id_surf, 'root_area_dens_s',                      &
1985                             root_area_density_lsm_f%var,                      &
1986                             nxl, nxr, nys, nyn,                               &
1987                             0, root_area_density_lsm_f%nz-1 )
1988
1989       ELSE
1990          root_area_density_lsm_f%from_file = .FALSE.
1991       ENDIF
1992!
1993!--    Read street type and street crossing
1994       IF ( check_existence( var_names, 'street_type' ) )  THEN
1995          street_type_f%from_file = .TRUE.
1996          CALL get_attribute( id_surf, char_fill,                              &
1997                              street_type_f%fill, .FALSE.,                     &
1998                              'street_type' )
1999
2000          ALLOCATE ( street_type_f%var(nys:nyn,nxl:nxr)  )
2001         
2002          CALL get_variable( id_surf, 'street_type', street_type_f%var,        &
2003                             nxl, nxr, nys, nyn )
2004       ELSE
2005          street_type_f%from_file = .FALSE.
2006       ENDIF
2007
2008       IF ( check_existence( var_names, 'street_crossing' ) )  THEN
2009          street_crossing_f%from_file = .TRUE.
2010          CALL get_attribute( id_surf, char_fill,                              &
2011                              street_crossing_f%fill, .FALSE.,                 &
2012                              'street_crossing' )
2013
2014          ALLOCATE ( street_crossing_f%var(nys:nyn,nxl:nxr)  )
2015
2016          CALL get_variable( id_surf, 'street_crossing',                       &
2017                             street_crossing_f%var, nxl, nxr, nys, nyn )
2018
2019       ELSE
2020          street_crossing_f%from_file = .FALSE.
2021       ENDIF
2022!
2023!--    Still missing: root_resolved and building_surface_pars.
2024!--    Will be implemented as soon as they are available.
2025
2026!
2027!--    Finally, close input file
2028       CALL close_input_file( id_surf )
2029#endif
2030!
2031!--    End of CPU measurement
2032       CALL cpu_log( log_point_s(82), 'NetCDF input', 'stop' )
2033!
2034!--    Exchange ghost points for surface variables. Therefore, resize
2035!--    variables.
2036       IF ( albedo_type_f%from_file )  THEN
2037          CALL resize_array_2d_int8( albedo_type_f%var, nys, nyn, nxl, nxr )
2038          CALL exchange_horiz_2d_byte( albedo_type_f%var, nys, nyn, nxl, nxr,  &
2039                                       nbgp )
2040       ENDIF
2041       IF ( pavement_type_f%from_file )  THEN
2042          CALL resize_array_2d_int8( pavement_type_f%var, nys, nyn, nxl, nxr )
2043          CALL exchange_horiz_2d_byte( pavement_type_f%var, nys, nyn, nxl, nxr,&
2044                                       nbgp )
2045       ENDIF
2046       IF ( soil_type_f%from_file  .AND.  ALLOCATED( soil_type_f%var_2d ) )  THEN
2047          CALL resize_array_2d_int8( soil_type_f%var_2d, nys, nyn, nxl, nxr )
2048          CALL exchange_horiz_2d_byte( soil_type_f%var_2d, nys, nyn, nxl, nxr, &
2049                                       nbgp )
2050       ENDIF
2051       IF ( vegetation_type_f%from_file )  THEN
2052          CALL resize_array_2d_int8( vegetation_type_f%var, nys, nyn, nxl, nxr )
2053          CALL exchange_horiz_2d_byte( vegetation_type_f%var, nys, nyn, nxl,   &
2054                                       nxr, nbgp )
2055       ENDIF
2056       IF ( water_type_f%from_file )  THEN
2057          CALL resize_array_2d_int8( water_type_f%var, nys, nyn, nxl, nxr )
2058          CALL exchange_horiz_2d_byte( water_type_f%var, nys, nyn, nxl, nxr,   &
2059                                       nbgp )
2060       ENDIF
2061!
2062!--    Exchange ghost points for 3/4-D variables. For the sake of simplicity,
2063!--    loop further dimensions to use 2D exchange routines. Unfortunately this
2064!--    is necessary, else new MPI-data types need to be introduced just for
2065!--    2 variables.
2066       IF ( soil_type_f%from_file  .AND.  ALLOCATED( soil_type_f%var_3d ) )    &
2067       THEN
2068          CALL resize_array_3d_int8( soil_type_f%var_3d, 0, nz_soil,           &
2069                                     nys, nyn, nxl, nxr )
2070          DO  k = 0, nz_soil
2071             CALL exchange_horiz_2d_byte(                                       &
2072                        soil_type_f%var_3d(k,:,:), nys, nyn, nxl, nxr, nbgp )
2073          ENDDO
2074       ENDIF
2075
2076       IF ( surface_fraction_f%from_file )  THEN
2077          CALL resize_array_3d_real( surface_fraction_f%frac,                  &
2078                                     0, surface_fraction_f%nf-1,               &
2079                                     nys, nyn, nxl, nxr )
2080          DO  k = 0, surface_fraction_f%nf-1
2081             CALL exchange_horiz_2d( surface_fraction_f%frac(k,:,:) )
2082          ENDDO
2083       ENDIF
2084
2085       IF ( building_pars_f%from_file )  THEN         
2086          CALL resize_array_3d_real( building_pars_f%pars_xy,                  &
2087                                     0, building_pars_f%np-1,                  &
2088                                     nys, nyn, nxl, nxr )
2089          DO  k = 0, building_pars_f%np-1
2090             CALL exchange_horiz_2d( building_pars_f%pars_xy(k,:,:) )
2091          ENDDO
2092       ENDIF
2093
2094       IF ( albedo_pars_f%from_file )  THEN         
2095          CALL resize_array_3d_real( albedo_pars_f%pars_xy,                    &
2096                                     0, albedo_pars_f%np-1,                    &
2097                                     nys, nyn, nxl, nxr )
2098          DO  k = 0, albedo_pars_f%np-1
2099             CALL exchange_horiz_2d( albedo_pars_f%pars_xy(k,:,:) )
2100          ENDDO
2101       ENDIF
2102
2103       IF ( pavement_pars_f%from_file )  THEN         
2104          CALL resize_array_3d_real( pavement_pars_f%pars_xy,                  &
2105                                     0, pavement_pars_f%np-1,                  &
2106                                     nys, nyn, nxl, nxr )
2107          DO  k = 0, pavement_pars_f%np-1
2108             CALL exchange_horiz_2d( pavement_pars_f%pars_xy(k,:,:) )
2109          ENDDO
2110       ENDIF
2111
2112       IF ( vegetation_pars_f%from_file )  THEN
2113          CALL resize_array_3d_real( vegetation_pars_f%pars_xy,                &
2114                                     0, vegetation_pars_f%np-1,                &
2115                                     nys, nyn, nxl, nxr )
2116          DO  k = 0, vegetation_pars_f%np-1
2117             CALL exchange_horiz_2d( vegetation_pars_f%pars_xy(k,:,:) )
2118          ENDDO
2119       ENDIF
2120
2121       IF ( water_pars_f%from_file )  THEN
2122          CALL resize_array_3d_real( water_pars_f%pars_xy,                     &
2123                                     0, water_pars_f%np-1,                     &
2124                                     nys, nyn, nxl, nxr )
2125          DO  k = 0, water_pars_f%np-1
2126             CALL exchange_horiz_2d( water_pars_f%pars_xy(k,:,:) )
2127          ENDDO
2128       ENDIF
2129
2130       IF ( root_area_density_lsm_f%from_file )  THEN
2131          CALL resize_array_3d_real( root_area_density_lsm_f%var,              &
2132                                     0, root_area_density_lsm_f%nz-1,          &
2133                                     nys, nyn, nxl, nxr )
2134          DO  k = 0, root_area_density_lsm_f%nz-1
2135             CALL exchange_horiz_2d( root_area_density_lsm_f%var(k,:,:) )
2136          ENDDO
2137       ENDIF
2138
2139       IF ( soil_pars_f%from_file )  THEN
2140          IF ( soil_pars_f%lod == 1 )  THEN
2141         
2142             CALL resize_array_3d_real( soil_pars_f%pars_xy,                   &
2143                                        0, soil_pars_f%np-1,                   &
2144                                        nys, nyn, nxl, nxr )
2145             DO  k = 0, soil_pars_f%np-1
2146                CALL exchange_horiz_2d( soil_pars_f%pars_xy(k,:,:) )
2147             ENDDO
2148             
2149          ELSEIF ( soil_pars_f%lod == 2 )  THEN
2150             CALL resize_array_4d_real( soil_pars_f%pars_xyz,                  &
2151                                        0, soil_pars_f%np-1,                   &
2152                                        0, soil_pars_f%nz-1,                   &
2153                                        nys, nyn, nxl, nxr )
2154
2155             DO  k2 = 0, soil_pars_f%nz-1
2156                DO  k = 0, soil_pars_f%np-1
2157                   CALL exchange_horiz_2d( soil_pars_f%pars_xyz(k,k2,:,:) )
2158                ENDDO
2159             ENDDO
2160          ENDIF
2161       ENDIF
2162
2163       IF ( pavement_subsurface_pars_f%from_file )  THEN         
2164          CALL resize_array_4d_real( pavement_subsurface_pars_f%pars_xyz,      &
2165                                     0, pavement_subsurface_pars_f%np-1,       &
2166                                     0, pavement_subsurface_pars_f%nz-1,       &
2167                                     nys, nyn, nxl, nxr )
2168
2169          DO  k2 = 0, pavement_subsurface_pars_f%nz-1
2170             DO  k = 0, pavement_subsurface_pars_f%np-1
2171                CALL exchange_horiz_2d(                                        &
2172                           pavement_subsurface_pars_f%pars_xyz(k,k2,:,:) )
2173             ENDDO
2174          ENDDO
2175       ENDIF
2176
2177    END SUBROUTINE netcdf_data_input_surface_data
2178
2179!------------------------------------------------------------------------------!
2180! Description:
2181! ------------
2182!> Reads uvem lookup table information.
2183!------------------------------------------------------------------------------!
2184    SUBROUTINE netcdf_data_input_uvem
2185       
2186       USE indices,                                                            &
2187           ONLY:  nxl, nxr, nyn, nys
2188
2189       IMPLICIT NONE
2190
2191       CHARACTER(LEN=100), DIMENSION(:), ALLOCATABLE ::  var_names  !< variable names in static input file
2192
2193
2194       INTEGER(iwp) ::  id_uvem       !< NetCDF id of uvem lookup table input file
2195       INTEGER(iwp) ::  nli = 35      !< dimension length of lookup table in x
2196       INTEGER(iwp) ::  nlj =  9      !< dimension length of lookup table in y
2197       INTEGER(iwp) ::  nlk = 90      !< dimension length of lookup table in z
2198       INTEGER(iwp) ::  num_vars      !< number of variables in netcdf input file
2199!
2200!--    Input via uv exposure model lookup table input
2201       IF ( input_pids_uvem )  THEN
2202
2203#if defined ( __netcdf )
2204!
2205!--       Open file in read-only mode
2206          CALL open_read_file( TRIM( input_file_uvem ) //                    &
2207                               TRIM( coupling_char ), id_uvem )
2208!
2209!--       At first, inquire all variable names.
2210!--       This will be used to check whether an input variable exist or not.
2211          CALL inquire_num_variables( id_uvem, num_vars )
2212!
2213!--       Allocate memory to store variable names and inquire them.
2214          ALLOCATE( var_names(1:num_vars) )
2215          CALL inquire_variable_names( id_uvem, var_names )
2216!
2217!--       uvem integration
2218          IF ( check_existence( var_names, 'int_factors' ) )  THEN
2219             uvem_integration_f%from_file = .TRUE.
2220!
2221!--          Input 2D uvem integration.
2222             ALLOCATE ( uvem_integration_f%var(0:nlj,0:nli)  )
2223             
2224             CALL get_variable( id_uvem, 'int_factors', uvem_integration_f%var, 0, nli, 0, nlj )
2225          ELSE
2226             uvem_integration_f%from_file = .FALSE.
2227          ENDIF
2228!
2229!--       uvem irradiance
2230          IF ( check_existence( var_names, 'irradiance' ) )  THEN
2231             uvem_irradiance_f%from_file = .TRUE.
2232!
2233!--          Input 2D uvem irradiance.
2234             ALLOCATE ( uvem_irradiance_f%var(0:nlk, 0:2)  )
2235             
2236             CALL get_variable( id_uvem, 'irradiance', uvem_irradiance_f%var, 0, 2, 0, nlk )
2237          ELSE
2238             uvem_irradiance_f%from_file = .FALSE.
2239          ENDIF
2240!
2241!--       uvem porjection areas
2242          IF ( check_existence( var_names, 'projarea' ) )  THEN
2243             uvem_projarea_f%from_file = .TRUE.
2244!
2245!--          Input 3D uvem projection area (human geometgry)
2246             ALLOCATE ( uvem_projarea_f%var(0:2,0:nlj,0:nli)  )
2247           
2248             CALL get_variable( id_uvem, 'projarea', uvem_projarea_f%var, 0, nli, 0, nlj, 0, 2 )
2249          ELSE
2250             uvem_projarea_f%from_file = .FALSE.
2251          ENDIF
2252!
2253!--       uvem radiance
2254          IF ( check_existence( var_names, 'radiance' ) )  THEN
2255             uvem_radiance_f%from_file = .TRUE.
2256!
2257!--          Input 3D uvem radiance
2258             ALLOCATE ( uvem_radiance_f%var(0:nlk,0:nlj,0:nli)  )
2259             
2260             CALL get_variable( id_uvem, 'radiance', uvem_radiance_f%var, 0, nli, 0, nlj, 0, nlk )
2261          ELSE
2262             uvem_radiance_f%from_file = .FALSE.
2263          ENDIF
2264!
2265!--       Read building obstruction
2266          IF ( check_existence( var_names, 'obstruction' ) )  THEN
2267             building_obstruction_full%from_file = .TRUE.
2268!--          Input 3D uvem building obstruction
2269              ALLOCATE ( building_obstruction_full%var_3d(0:44,0:2,0:2) )
2270              CALL get_variable( id_uvem, 'obstruction', building_obstruction_full%var_3d,0, 2, 0, 2, 0, 44 )       
2271          ELSE
2272             building_obstruction_full%from_file = .FALSE.
2273          ENDIF
2274!
2275          IF ( check_existence( var_names, 'obstruction' ) )  THEN
2276             building_obstruction_f%from_file = .TRUE.
2277!
2278!--          Input 3D uvem building obstruction
2279             ALLOCATE ( building_obstruction_f%var_3d(0:44,nys:nyn,nxl:nxr) )
2280!
2281             CALL get_variable( id_uvem, 'obstruction', building_obstruction_f%var_3d,      &
2282                                nxl, nxr, nys, nyn, 0, 44 )       
2283          ELSE
2284             building_obstruction_f%from_file = .FALSE.
2285          ENDIF
2286!
2287!--       Close uvem lookup table input file
2288          CALL close_input_file( id_uvem )
2289#else
2290          CONTINUE
2291#endif
2292       ENDIF
2293    END SUBROUTINE netcdf_data_input_uvem
2294
2295!------------------------------------------------------------------------------!
2296! Description:
2297! ------------
2298!> Reads orography and building information.
2299!------------------------------------------------------------------------------!
2300    SUBROUTINE netcdf_data_input_topo
2301
2302       USE control_parameters,                                                 &
2303           ONLY:  message_string, topography
2304
2305       USE exchange_horiz_mod,                                                 &
2306           ONLY:  exchange_horiz_2d_byte, exchange_horiz_2d_int
2307
2308       USE grid_variables,                                                     &
2309           ONLY:  dx, dy   
2310           
2311       USE indices,                                                            &
2312           ONLY:  nbgp, nx, nxl, nxr, ny, nyn, nys, nzb
2313
2314
2315       IMPLICIT NONE
2316
2317       CHARACTER(LEN=100), DIMENSION(:), ALLOCATABLE ::  var_names  !< variable names in static input file
2318
2319
2320       INTEGER(iwp) ::  i             !< running index along x-direction
2321       INTEGER(iwp) ::  ii            !< running index for IO blocks
2322       INTEGER(iwp) ::  id_topo       !< NetCDF id of topograhy input file
2323       INTEGER(iwp) ::  io_status     !< status after reading the ascii topo file
2324       INTEGER(iwp) ::  j             !< running index along y-direction
2325       INTEGER(iwp) ::  num_vars      !< number of variables in netcdf input file
2326       INTEGER(iwp) ::  skip_n_rows   !< counting variable to skip rows while reading topography file
2327
2328       REAL(wp) ::  dum           !< dummy variable to skip columns while reading topography file
2329!
2330!--    CPU measurement
2331       CALL cpu_log( log_point_s(83), 'NetCDF/ASCII input topo', 'start' )
2332
2333!
2334!--    Input via palm-input data standard
2335       IF ( input_pids_static )  THEN
2336#if defined ( __netcdf )
2337!
2338!--       Open file in read-only mode
2339          CALL open_read_file( TRIM( input_file_static ) //                    &
2340                               TRIM( coupling_char ), id_topo )
2341!
2342!--       At first, inquire all variable names.
2343!--       This will be used to check whether an  input variable exist
2344!--       or not.
2345          CALL inquire_num_variables( id_topo, num_vars )
2346!
2347!--       Allocate memory to store variable names and inquire them.
2348          ALLOCATE( var_names(1:num_vars) )
2349          CALL inquire_variable_names( id_topo, var_names )
2350!
2351!--       Read x, y - dimensions. Only required for consistency checks.
2352          CALL get_dimension_length( id_topo, dim_static%nx, 'x' )
2353          CALL get_dimension_length( id_topo, dim_static%ny, 'y' )
2354          ALLOCATE( dim_static%x(0:dim_static%nx-1) )
2355          ALLOCATE( dim_static%y(0:dim_static%ny-1) )
2356          CALL get_variable( id_topo, 'x', dim_static%x )
2357          CALL get_variable( id_topo, 'y', dim_static%y )
2358!
2359!--       Check whether dimension size in input file matches the model dimensions
2360          IF ( dim_static%nx-1 /= nx  .OR.  dim_static%ny-1 /= ny )  THEN
2361             message_string = 'Static input file: horizontal dimension in ' // &
2362                              'x- and/or y-direction ' //                      &
2363                              'do not match the respective model dimension'
2364             CALL message( 'netcdf_data_input_mod', 'PA0548', 1, 2, 0, 6, 0 )
2365          ENDIF
2366!
2367!--       Check if grid spacing of provided input data matches the respective
2368!--       grid spacing in the model.
2369          IF ( ABS( dim_static%x(1) - dim_static%x(0) - dx ) > 10E-6_wp  .OR.  &
2370               ABS( dim_static%y(1) - dim_static%y(0) - dy ) > 10E-6_wp )  THEN
2371             message_string = 'Static input file: horizontal grid spacing ' // &
2372                              'in x- and/or y-direction ' //                   &
2373                              'do not match the respective model grid spacing.'
2374             CALL message( 'netcdf_data_input_mod', 'PA0549', 1, 2, 0, 6, 0 )
2375          ENDIF
2376!
2377!--       Terrain height. First, get variable-related _FillValue attribute
2378          IF ( check_existence( var_names, 'zt' ) )  THEN
2379             terrain_height_f%from_file = .TRUE.
2380             CALL get_attribute( id_topo, char_fill, terrain_height_f%fill,    &
2381                                 .FALSE., 'zt' )
2382!
2383!--          Input 2D terrain height.
2384             ALLOCATE ( terrain_height_f%var(nys:nyn,nxl:nxr)  )
2385             
2386             CALL get_variable( id_topo, 'zt', terrain_height_f%var,           &
2387                                nxl, nxr, nys, nyn )
2388
2389          ELSE
2390             terrain_height_f%from_file = .FALSE.
2391          ENDIF
2392
2393!
2394!--       Read building height. First, read its _FillValue attribute,
2395!--       as well as lod attribute
2396          buildings_f%from_file = .FALSE.
2397          IF ( check_existence( var_names, 'buildings_2d' ) )  THEN
2398             buildings_f%from_file = .TRUE.
2399             CALL get_attribute( id_topo, char_lod, buildings_f%lod,           &
2400                                 .FALSE., 'buildings_2d' )
2401
2402             CALL get_attribute( id_topo, char_fill, buildings_f%fill1,        &
2403                                 .FALSE., 'buildings_2d' )
2404
2405!
2406!--          Read 2D buildings
2407             IF ( buildings_f%lod == 1 )  THEN
2408                ALLOCATE ( buildings_f%var_2d(nys:nyn,nxl:nxr) )
2409
2410                CALL get_variable( id_topo, 'buildings_2d',                    &
2411                                   buildings_f%var_2d,                         &
2412                                   nxl, nxr, nys, nyn )
2413             ELSE
2414                message_string = 'NetCDF attribute lod ' //                    &
2415                                 '(level of detail) is not set ' //            &
2416                                 'properly for buildings_2d.'
2417                CALL message( 'netcdf_data_input_mod', 'PA0540',               &
2418                               1, 2, 0, 6, 0 )
2419             ENDIF
2420          ENDIF
2421!
2422!--       If available, also read 3D building information. If both are
2423!--       available, use 3D information.
2424          IF ( check_existence( var_names, 'buildings_3d' ) )  THEN
2425             buildings_f%from_file = .TRUE.
2426             CALL get_attribute( id_topo, char_lod, buildings_f%lod,           &
2427                                 .FALSE., 'buildings_3d' )     
2428
2429             CALL get_attribute( id_topo, char_fill, buildings_f%fill2,        &
2430                                 .FALSE., 'buildings_3d' )
2431
2432             CALL get_dimension_length( id_topo, buildings_f%nz, 'z' )
2433!
2434!--          Read 3D buildings
2435             IF ( buildings_f%lod == 2 )  THEN
2436                ALLOCATE( buildings_f%z(nzb:buildings_f%nz-1) )
2437                CALL get_variable( id_topo, 'z', buildings_f%z )
2438
2439                ALLOCATE( buildings_f%var_3d(nzb:buildings_f%nz-1,             &
2440                                             nys:nyn,nxl:nxr) )
2441                buildings_f%var_3d = 0
2442               
2443                CALL get_variable( id_topo, 'buildings_3d',                    &
2444                                   buildings_f%var_3d,                         &
2445                                   nxl, nxr, nys, nyn, 0, buildings_f%nz-1 )
2446             ELSE
2447                message_string = 'NetCDF attribute lod ' //                    &
2448                                 '(level of detail) is not set ' //            &
2449                                 'properly for buildings_3d.'
2450                CALL message( 'netcdf_data_input_mod', 'PA0541',               &
2451                               1, 2, 0, 6, 0 )
2452             ENDIF
2453          ENDIF
2454!
2455!--       Read building IDs and its FillValue attribute. Further required
2456!--       for mapping buildings on top of orography.
2457          IF ( check_existence( var_names, 'building_id' ) )  THEN
2458             building_id_f%from_file = .TRUE.
2459             CALL get_attribute( id_topo, char_fill,                           &
2460                                 building_id_f%fill, .FALSE.,                  &
2461                                 'building_id' )
2462
2463             ALLOCATE ( building_id_f%var(nys:nyn,nxl:nxr) )
2464             
2465             CALL get_variable( id_topo, 'building_id', building_id_f%var,     &
2466                                nxl, nxr, nys, nyn )
2467          ELSE
2468             building_id_f%from_file = .FALSE.
2469          ENDIF
2470!
2471!--       Read building_type and required attributes.
2472          IF ( check_existence( var_names, 'building_type' ) )  THEN
2473             building_type_f%from_file = .TRUE.
2474             CALL get_attribute( id_topo, char_fill,                           &
2475                                 building_type_f%fill, .FALSE.,                &
2476                                 'building_type' )
2477
2478             ALLOCATE ( building_type_f%var(nys:nyn,nxl:nxr) )
2479
2480             CALL get_variable( id_topo, 'building_type', building_type_f%var, &
2481                                nxl, nxr, nys, nyn )
2482
2483          ELSE
2484             building_type_f%from_file = .FALSE.
2485          ENDIF
2486!
2487!--       Close topography input file
2488          CALL close_input_file( id_topo )
2489#else
2490          CONTINUE
2491#endif
2492!
2493!--    ASCII input
2494       ELSEIF ( TRIM( topography ) == 'read_from_file' )  THEN
2495             
2496          DO  ii = 0, io_blocks-1
2497             IF ( ii == io_group )  THEN
2498
2499                OPEN( 90, FILE='TOPOGRAPHY_DATA'//TRIM( coupling_char ),       &
2500                      STATUS='OLD', FORM='FORMATTED', IOSTAT=io_status )
2501
2502                IF ( io_status > 0 )  THEN
2503                   message_string = 'file TOPOGRAPHY_DATA'//                      &
2504                                    TRIM( coupling_char )// ' does not exist'
2505                   CALL message( 'netcdf_data_input_mod', 'PA0208', 1, 2, 0, 6, 0 )
2506                ENDIF
2507
2508!
2509!--             Read topography PE-wise. Rows are read from nyn to nys, columns
2510!--             are read from nxl to nxr. At first, ny-nyn rows need to be skipped.
2511                skip_n_rows = 0
2512                DO WHILE ( skip_n_rows < ny - nyn )
2513                   READ( 90, * )
2514                   skip_n_rows = skip_n_rows + 1
2515                ENDDO
2516!
2517!--             Read data from nyn to nys and nxl to nxr. Therefore, skip
2518!--             column until nxl-1 is reached
2519                ALLOCATE ( buildings_f%var_2d(nys:nyn,nxl:nxr) )
2520                DO  j = nyn, nys, -1
2521
2522                   READ( 90, *, IOSTAT=io_status )                               &
2523                                   ( dum, i = 0, nxl-1 ),                      &
2524                                   ( buildings_f%var_2d(j,i), i = nxl, nxr )
2525
2526                   IF ( io_status > 0 )  THEN
2527                      WRITE( message_string, '(A,1X,I5,1X,A)' ) 'error reading line', ny-j+1,      &
2528                                                   'of file TOPOGRAPHY_DATA'//TRIM( coupling_char )
2529                      CALL message( 'netcdf_data_input_mod', 'PA0209', 2, 2, myid, 6, 0 )
2530                   ELSEIF ( io_status < 0 )  THEN
2531                      WRITE( message_string, '(A,1X,I5)' ) 'end of line or file detected for '// &
2532                                  'file TOPOGRAPHY_DATA'//TRIM( coupling_char )//' at line', ny-j+1
2533                      CALL message( 'netcdf_data_input_mod', 'PA0704', 2, 2, myid, 6, 0 )
2534                   ENDIF
2535
2536                ENDDO
2537
2538                CLOSE( 90 )
2539                buildings_f%from_file = .TRUE.
2540
2541             ENDIF
2542#if defined( __parallel )
2543             CALL MPI_BARRIER( comm2d, ierr )
2544#endif
2545          ENDDO
2546
2547       ENDIF
2548!
2549!--    End of CPU measurement
2550       CALL cpu_log( log_point_s(83), 'NetCDF/ASCII input topo', 'stop' )
2551!
2552!--    Check for minimum requirement to setup building topography. If buildings
2553!--    are provided, also an ID and a type are required.
2554!--    Note, doing this check in check_parameters
2555!--    will be too late (data will be used for grid inititialization before).
2556       IF ( input_pids_static )  THEN
2557          IF ( buildings_f%from_file  .AND.                                    &
2558               .NOT. building_id_f%from_file )  THEN
2559             message_string = 'If building heights are prescribed in ' //      &
2560                              'static input file, also an ID is required.'
2561             CALL message( 'netcdf_data_input_mod', 'PA0542', 1, 2, 0, 6, 0 )
2562          ENDIF
2563       ENDIF
2564
2565       IF ( terrain_height_f%from_file )  THEN
2566!
2567!--       Check orography for fill-values.
2568!--       For the moment, give an error message. More advanced methods, e.g. a
2569!--       nearest neighbor algorithm as used in GIS systems might be implemented
2570!--       later.
2571!--       NOTE: This check must be placed here as terrain_height_f is altered
2572!--       within init_grid which is called before netcdf_data_input_check_static
2573          IF ( ANY( terrain_height_f%var == terrain_height_f%fill ) )  THEN
2574             message_string = 'NetCDF variable zt is not ' //                  &
2575                              'allowed to have missing data'
2576             CALL message( 'netcdf_data_input_mod', 'PA0550', 2, 2, myid, 6, 0 )
2577          ENDIF
2578       ELSE
2579!
2580!--       In case no terrain height is provided by static input file, allocate
2581!--       array nevertheless and set terrain height to 0, which simplifies
2582!--       topography initialization.
2583          ALLOCATE ( terrain_height_f%var(nys:nyn,nxl:nxr) )
2584          terrain_height_f%var = 0.0_wp
2585       ENDIF
2586!
2587!--    Finally, exchange 1 ghost point for building ID and type.
2588!--    In case of non-cyclic boundary conditions set Neumann conditions at the
2589!--    lateral boundaries.
2590       IF ( building_id_f%from_file )  THEN
2591          CALL resize_array_2d_int32( building_id_f%var, nys, nyn, nxl, nxr )
2592          CALL exchange_horiz_2d_int( building_id_f%var, nys, nyn, nxl, nxr,   &
2593                                      nbgp )
2594       ENDIF
2595
2596       IF ( building_type_f%from_file )  THEN
2597          CALL resize_array_2d_int8( building_type_f%var, nys, nyn, nxl, nxr )
2598          CALL exchange_horiz_2d_byte( building_type_f%var, nys, nyn, nxl, nxr,   &
2599                                       nbgp )
2600       ENDIF
2601
2602    END SUBROUTINE netcdf_data_input_topo
2603
2604!------------------------------------------------------------------------------!
2605! Description:
2606! ------------
2607!> Reads initialization data of u, v, w, pt, q, geostrophic wind components,
2608!> as well as soil moisture and soil temperature, derived from larger-scale
2609!> model (COSMO) by Inifor.
2610!------------------------------------------------------------------------------!
2611    SUBROUTINE netcdf_data_input_init_3d
2612
2613       USE arrays_3d,                                                          &
2614           ONLY:  q, pt, u, v, w, zu, zw
2615
2616       USE control_parameters,                                                 &
2617           ONLY:  air_chemistry, bc_lr_cyc, bc_ns_cyc, humidity,               &
2618                  message_string, neutral
2619
2620       USE indices,                                                            &
2621           ONLY:  nx, nxl, nxlu, nxr, ny, nyn, nys, nysv, nzb, nz, nzt
2622
2623       IMPLICIT NONE
2624
2625       CHARACTER(LEN=100), DIMENSION(:), ALLOCATABLE ::  var_names
2626
2627       LOGICAL      ::  dynamic_3d = .TRUE. !< flag indicating that 3D data is read from dynamic file
2628       
2629       INTEGER(iwp) ::  id_dynamic !< NetCDF id of dynamic input file
2630       INTEGER(iwp) ::  n          !< running index for chemistry variables
2631       INTEGER(iwp) ::  num_vars   !< number of variables in netcdf input file
2632
2633       LOGICAL      ::  check_passed !< flag indicating if a check passed
2634
2635!
2636!--    Skip routine if no input file with dynamic input data is available.
2637       IF ( .NOT. input_pids_dynamic )  RETURN
2638!
2639!--    Please note, Inifor is designed to provide initial data for u and v for
2640!--    the prognostic grid points in case of lateral Dirichlet conditions.
2641!--    This means that Inifor provides data from nxlu:nxr (for u) and
2642!--    from nysv:nyn (for v) at the left and south domain boundary, respectively.
2643!--    However, as work-around for the moment, PALM will run with cyclic
2644!--    conditions and will be initialized with data provided by Inifor
2645!--    boundaries in case of Dirichlet.
2646!--    Hence, simply set set nxlu/nysv to 1 (will be reset to its original value
2647!--    at the end of this routine.
2648       IF ( bc_lr_cyc  .AND.  nxl == 0 )  nxlu = 1
2649       IF ( bc_ns_cyc  .AND.  nys == 0 )  nysv = 1
2650
2651!
2652!--    CPU measurement
2653       CALL cpu_log( log_point_s(85), 'NetCDF input init', 'start' )
2654
2655#if defined ( __netcdf )
2656!
2657!--    Open file in read-only mode
2658       CALL open_read_file( TRIM( input_file_dynamic ) //                      &
2659                            TRIM( coupling_char ), id_dynamic )
2660
2661!
2662!--    At first, inquire all variable names.
2663       CALL inquire_num_variables( id_dynamic, num_vars )
2664!
2665!--    Allocate memory to store variable names.
2666       ALLOCATE( var_names(1:num_vars) )
2667       CALL inquire_variable_names( id_dynamic, var_names )
2668!
2669!--    Read vertical dimension of scalar und w grid.
2670       CALL get_dimension_length( id_dynamic, init_3d%nzu, 'z'     )
2671       CALL get_dimension_length( id_dynamic, init_3d%nzw, 'zw'    )
2672!
2673!--    Read also the horizontal dimensions. These are used just used fo
2674!--    checking the compatibility with the PALM grid before reading.
2675       CALL get_dimension_length( id_dynamic, init_3d%nx,  'x'  )
2676       CALL get_dimension_length( id_dynamic, init_3d%nxu, 'xu' )
2677       CALL get_dimension_length( id_dynamic, init_3d%ny,  'y'  )
2678       CALL get_dimension_length( id_dynamic, init_3d%nyv, 'yv' )
2679
2680!
2681!--    Check for correct horizontal and vertical dimension. Please note,
2682!--    checks are performed directly here and not called from
2683!--    check_parameters as some varialbes are still not allocated there.
2684!--    Moreover, please note, u- and v-grid has 1 grid point less on
2685!--    Inifor grid.
2686       IF ( init_3d%nx-1 /= nx  .OR.  init_3d%nxu-1 /= nx - 1  .OR.            &
2687            init_3d%ny-1 /= ny  .OR.  init_3d%nyv-1 /= ny - 1 )  THEN
2688          message_string = 'Number of horizontal grid points in '//            &
2689                           'dynamic input file does not match ' //             &
2690                           'the number of numeric grid points.'
2691          CALL message( 'netcdf_data_input_mod', 'PA0543', 1, 2, 0, 6, 0 )
2692       ENDIF
2693
2694       IF ( init_3d%nzu /= nz )  THEN
2695          message_string = 'Number of vertical grid points in '//              &
2696                           'dynamic input file does not match ' //             &
2697                           'the number of numeric grid points.'
2698          CALL message( 'netcdf_data_input_mod', 'PA0543', 1, 2, 0, 6, 0 )
2699       ENDIF
2700!
2701!--    Read vertical dimensions. Later, these are required for eventual
2702!--    inter- and extrapolations of the initialization data.
2703       IF ( check_existence( var_names, 'z' ) )  THEN
2704          ALLOCATE( init_3d%zu_atmos(1:init_3d%nzu) )
2705          CALL get_variable( id_dynamic, 'z', init_3d%zu_atmos )
2706       ENDIF
2707       IF ( check_existence( var_names, 'zw' ) )  THEN
2708          ALLOCATE( init_3d%zw_atmos(1:init_3d%nzw) )
2709          CALL get_variable( id_dynamic, 'zw', init_3d%zw_atmos )
2710       ENDIF
2711!
2712!--    Check for consistency between vertical coordinates in dynamic
2713!--    driver and numeric grid.
2714!--    Please note, depending on compiler options both may be
2715!--    equal up to a certain threshold, and differences between
2716!--    the numeric grid and vertical coordinate in the driver can built-
2717!--    up to 10E-1-10E-0 m. For this reason, the check is performed not
2718!--    for exactly matching values.
2719       IF ( ANY( ABS( zu(1:nzt)   - init_3d%zu_atmos(1:init_3d%nzu) )    &
2720                      > 10E-1 )  .OR.                                    &
2721            ANY( ABS( zw(1:nzt-1) - init_3d%zw_atmos(1:init_3d%nzw) )    &
2722                      > 10E-1 ) )  THEN
2723          message_string = 'Vertical grid in dynamic driver does not '// &
2724                           'match the numeric grid.'
2725          CALL message( 'netcdf_data_input_mod', 'PA0543', 1, 2, 0, 6, 0 )
2726       ENDIF
2727!
2728!--    Read initial geostrophic wind components at
2729!--    t = 0 (index 1 in file).
2730       IF ( check_existence( var_names, 'ls_forcing_ug' ) )  THEN
2731          ALLOCATE( init_3d%ug_init(nzb:nzt+1) )
2732          init_3d%ug_init = 0.0_wp
2733
2734          CALL get_variable_pr( id_dynamic, 'ls_forcing_ug', 1,          &
2735                                init_3d%ug_init(1:nzt) )
2736!
2737!--       Set top-boundary condition (Neumann)
2738          init_3d%ug_init(nzt+1) = init_3d%ug_init(nzt)
2739
2740          init_3d%from_file_ug = .TRUE.
2741       ELSE
2742          init_3d%from_file_ug = .FALSE.
2743       ENDIF
2744       IF ( check_existence( var_names, 'ls_forcing_vg' ) )  THEN
2745          ALLOCATE( init_3d%vg_init(nzb:nzt+1) )
2746          init_3d%vg_init = 0.0_wp
2747
2748          CALL get_variable_pr( id_dynamic, 'ls_forcing_vg', 1,          &
2749                                init_3d%vg_init(1:nzt) )
2750!
2751!--       Set top-boundary condition (Neumann)
2752          init_3d%vg_init(nzt+1) = init_3d%vg_init(nzt)
2753
2754          init_3d%from_file_vg = .TRUE.
2755       ELSE
2756          init_3d%from_file_vg = .FALSE.
2757       ENDIF
2758!
2759!--    Read inital 3D data of u, v, w, pt and q,
2760!--    derived from COSMO model. Read PE-wise yz-slices.
2761!--    Please note, the u-, v- and w-component are defined on different
2762!--    grids with one element less in the x-, y-,
2763!--    and z-direction, respectively. Hence, reading is subdivided
2764!--    into separate loops. 
2765!--    Read u-component
2766       IF ( check_existence( var_names, 'init_atmosphere_u' ) )  THEN
2767!
2768!--       Read attributes for the fill value and level-of-detail
2769          CALL get_attribute( id_dynamic, char_fill, init_3d%fill_u,           &
2770                              .FALSE., 'init_atmosphere_u' )
2771          CALL get_attribute( id_dynamic, char_lod, init_3d%lod_u,             &
2772                              .FALSE., 'init_atmosphere_u' )
2773!
2774!--       level-of-detail 1 - read initialization profile
2775          IF ( init_3d%lod_u == 1 )  THEN
2776             ALLOCATE( init_3d%u_init(nzb:nzt+1) )
2777             init_3d%u_init = 0.0_wp
2778
2779             CALL get_variable( id_dynamic, 'init_atmosphere_u',               &
2780                                init_3d%u_init(nzb+1:nzt) )
2781!
2782!--          Set top-boundary condition (Neumann)
2783             init_3d%u_init(nzt+1) = init_3d%u_init(nzt)
2784!
2785!--       level-of-detail 2 - read 3D initialization data
2786          ELSEIF ( init_3d%lod_u == 2 )  THEN
2787             CALL get_variable( id_dynamic, 'init_atmosphere_u',               &
2788                                u(nzb+1:nzt,nys:nyn,nxlu:nxr),                 &
2789                                nxlu, nys+1, nzb+1,                            &
2790                                nxr-nxlu+1, nyn-nys+1, init_3d%nzu,            &
2791                                dynamic_3d )
2792!
2793!--          Set value at leftmost model grid point nxl = 0. This is because
2794!--          Inifor provides data only from 1:nx-1 since it assumes non-cyclic
2795!--          conditions.
2796             IF ( nxl == 0 )                                                   &
2797                u(nzb+1:nzt,nys:nyn,nxl) = u(nzb+1:nzt,nys:nyn,nxlu)
2798!
2799!--          Set bottom and top-boundary
2800             u(nzb,:,:)   = u(nzb+1,:,:)
2801             u(nzt+1,:,:) = u(nzt,:,:)
2802             
2803          ENDIF
2804          init_3d%from_file_u = .TRUE.
2805       ELSE
2806          message_string = 'Missing initial data for u-component'
2807          CALL message( 'netcdf_data_input_mod', 'PA0544', 1, 2, 0, 6, 0 )
2808       ENDIF
2809!
2810!--    Read v-component
2811       IF ( check_existence( var_names, 'init_atmosphere_v' ) )  THEN
2812!
2813!--       Read attributes for the fill value and level-of-detail
2814          CALL get_attribute( id_dynamic, char_fill, init_3d%fill_v,           &
2815                              .FALSE., 'init_atmosphere_v' )
2816          CALL get_attribute( id_dynamic, char_lod, init_3d%lod_v,             &
2817                              .FALSE., 'init_atmosphere_v' )
2818!
2819!--       level-of-detail 1 - read initialization profile
2820          IF ( init_3d%lod_v == 1 )  THEN
2821             ALLOCATE( init_3d%v_init(nzb:nzt+1) )
2822             init_3d%v_init = 0.0_wp
2823
2824             CALL get_variable( id_dynamic, 'init_atmosphere_v',               &
2825                                init_3d%v_init(nzb+1:nzt) )
2826!
2827!--          Set top-boundary condition (Neumann)
2828             init_3d%v_init(nzt+1) = init_3d%v_init(nzt)
2829!
2830!--       level-of-detail 2 - read 3D initialization data
2831          ELSEIF ( init_3d%lod_v == 2 )  THEN
2832         
2833             CALL get_variable( id_dynamic, 'init_atmosphere_v',               &
2834                                v(nzb+1:nzt,nysv:nyn,nxl:nxr),                 &
2835                                nxl+1, nysv, nzb+1,                            &
2836                                nxr-nxl+1, nyn-nysv+1, init_3d%nzu,            &
2837                                dynamic_3d )
2838!
2839!--          Set value at southmost model grid point nys = 0. This is because
2840!--          Inifor provides data only from 1:ny-1 since it assumes non-cyclic
2841!--          conditions.
2842             IF ( nys == 0 )                                                   &
2843                v(nzb+1:nzt,nys,nxl:nxr) = v(nzb+1:nzt,nysv,nxl:nxr)                               
2844!
2845!--          Set bottom and top-boundary
2846             v(nzb,:,:)   = v(nzb+1,:,:)
2847             v(nzt+1,:,:) = v(nzt,:,:)
2848             
2849          ENDIF
2850          init_3d%from_file_v = .TRUE.
2851       ELSE
2852          message_string = 'Missing initial data for v-component'
2853          CALL message( 'netcdf_data_input_mod', 'PA0544', 1, 2, 0, 6, 0 )
2854       ENDIF
2855!
2856!--    Read w-component
2857       IF ( check_existence( var_names, 'init_atmosphere_w' ) )  THEN
2858!
2859!--       Read attributes for the fill value and level-of-detail
2860          CALL get_attribute( id_dynamic, char_fill, init_3d%fill_w,           &
2861                              .FALSE., 'init_atmosphere_w' )
2862          CALL get_attribute( id_dynamic, char_lod, init_3d%lod_w,             &
2863                              .FALSE., 'init_atmosphere_w' )
2864!
2865!--       level-of-detail 1 - read initialization profile
2866          IF ( init_3d%lod_w == 1 )  THEN
2867             ALLOCATE( init_3d%w_init(nzb:nzt+1) )
2868             init_3d%w_init = 0.0_wp
2869
2870             CALL get_variable( id_dynamic, 'init_atmosphere_w',               &
2871                                init_3d%w_init(nzb+1:nzt-1) )
2872!
2873!--          Set top-boundary condition (Neumann)
2874             init_3d%w_init(nzt:nzt+1) = init_3d%w_init(nzt-1)
2875!
2876!--       level-of-detail 2 - read 3D initialization data
2877          ELSEIF ( init_3d%lod_w == 2 )  THEN
2878
2879             CALL get_variable( id_dynamic, 'init_atmosphere_w',                &
2880                                w(nzb+1:nzt-1,nys:nyn,nxl:nxr),                 &
2881                                nxl+1, nys+1, nzb+1,                            &
2882                                nxr-nxl+1, nyn-nys+1, init_3d%nzw,              &
2883                                dynamic_3d )
2884!
2885!--          Set bottom and top-boundary                               
2886             w(nzb,:,:)   = 0.0_wp 
2887             w(nzt,:,:)   = w(nzt-1,:,:)
2888             w(nzt+1,:,:) = w(nzt-1,:,:)
2889
2890          ENDIF
2891          init_3d%from_file_w = .TRUE.
2892       ELSE
2893          message_string = 'Missing initial data for w-component'
2894          CALL message( 'netcdf_data_input_mod', 'PA0544', 1, 2, 0, 6, 0 )
2895       ENDIF
2896!
2897!--    Read potential temperature
2898       IF ( .NOT. neutral )  THEN
2899          IF ( check_existence( var_names, 'init_atmosphere_pt' ) )  THEN
2900!
2901!--          Read attributes for the fill value and level-of-detail
2902             CALL get_attribute( id_dynamic, char_fill, init_3d%fill_pt,       &
2903                                 .FALSE., 'init_atmosphere_pt' )
2904             CALL get_attribute( id_dynamic, char_lod, init_3d%lod_pt,         &
2905                                 .FALSE., 'init_atmosphere_pt' )
2906!
2907!--          level-of-detail 1 - read initialization profile
2908             IF ( init_3d%lod_pt == 1 )  THEN
2909                ALLOCATE( init_3d%pt_init(nzb:nzt+1) )
2910
2911                CALL get_variable( id_dynamic, 'init_atmosphere_pt',           &
2912                                   init_3d%pt_init(nzb+1:nzt) )
2913!
2914!--             Set Neumann top and surface boundary condition for initial
2915!--             profil
2916                init_3d%pt_init(nzb)   = init_3d%pt_init(nzb+1)
2917                init_3d%pt_init(nzt+1) = init_3d%pt_init(nzt)
2918!
2919!--          level-of-detail 2 - read 3D initialization data
2920             ELSEIF ( init_3d%lod_pt == 2 )  THEN
2921
2922                CALL get_variable( id_dynamic, 'init_atmosphere_pt',           &
2923                                   pt(nzb+1:nzt,nys:nyn,nxl:nxr),              &
2924                                   nxl+1, nys+1, nzb+1,                        &
2925                                   nxr-nxl+1, nyn-nys+1, init_3d%nzu,          &
2926                                   dynamic_3d )
2927                                   
2928!
2929!--             Set bottom and top-boundary
2930                pt(nzb,:,:)   = pt(nzb+1,:,:)
2931                pt(nzt+1,:,:) = pt(nzt,:,:)             
2932
2933             ENDIF
2934             init_3d%from_file_pt = .TRUE.
2935          ELSE
2936             message_string = 'Missing initial data for ' //                   &
2937                              'potential temperature'
2938             CALL message( 'netcdf_data_input_mod', 'PA0544', 1, 2, 0, 6, 0 )
2939          ENDIF
2940       ENDIF
2941!
2942!--    Read mixing ratio
2943       IF ( humidity )  THEN
2944          IF ( check_existence( var_names, 'init_atmosphere_qv' ) )  THEN
2945!
2946!--          Read attributes for the fill value and level-of-detail
2947             CALL get_attribute( id_dynamic, char_fill, init_3d%fill_q,        &
2948                                 .FALSE., 'init_atmosphere_qv' )
2949             CALL get_attribute( id_dynamic, char_lod, init_3d%lod_q,          &
2950                                 .FALSE., 'init_atmosphere_qv' )
2951!
2952!--          level-of-detail 1 - read initialization profile
2953             IF ( init_3d%lod_q == 1 )  THEN
2954                ALLOCATE( init_3d%q_init(nzb:nzt+1) )
2955
2956                CALL get_variable( id_dynamic, 'init_atmosphere_qv',           &
2957                                    init_3d%q_init(nzb+1:nzt) )
2958!
2959!--             Set bottom and top boundary condition (Neumann)
2960                init_3d%q_init(nzb)   = init_3d%q_init(nzb+1)
2961                init_3d%q_init(nzt+1) = init_3d%q_init(nzt)
2962!
2963!--          level-of-detail 2 - read 3D initialization data
2964             ELSEIF ( init_3d%lod_q == 2 )  THEN
2965             
2966                CALL get_variable( id_dynamic, 'init_atmosphere_qv',           &
2967                                   q(nzb+1:nzt,nys:nyn,nxl:nxr),               &
2968                                   nxl+1, nys+1, nzb+1,                        &
2969                                   nxr-nxl+1, nyn-nys+1, init_3d%nzu,          &
2970                                   dynamic_3d )
2971                                   
2972!
2973!--             Set bottom and top-boundary
2974                q(nzb,:,:)   = q(nzb+1,:,:)
2975                q(nzt+1,:,:) = q(nzt,:,:)
2976               
2977             ENDIF
2978             init_3d%from_file_q = .TRUE.
2979          ELSE
2980             message_string = 'Missing initial data for ' //                   &
2981                              'mixing ratio'
2982             CALL message( 'netcdf_data_input_mod', 'PA0544', 1, 2, 0, 6, 0 )
2983          ENDIF
2984       ENDIF       
2985!
2986!--    Read chemistry variables.
2987!--    Please note, for the moment, only LOD=1 is allowed
2988       IF ( air_chemistry )  THEN
2989!
2990!--       Allocate chemistry input profiles, as well as arrays for fill values
2991!--       and LOD's.
2992          ALLOCATE( init_3d%chem_init(nzb:nzt+1,                               &
2993                                      1:UBOUND(init_3d%var_names_chem, 1 )) )
2994          ALLOCATE( init_3d%fill_chem(1:UBOUND(init_3d%var_names_chem, 1)) )   
2995          ALLOCATE( init_3d%lod_chem(1:UBOUND(init_3d%var_names_chem, 1))  ) 
2996         
2997          DO  n = 1, UBOUND(init_3d%var_names_chem, 1)
2998             IF ( check_existence( var_names,                                  &
2999                                   TRIM( init_3d%var_names_chem(n) ) ) )  THEN
3000!
3001!--             Read attributes for the fill value and level-of-detail
3002                CALL get_attribute( id_dynamic, char_fill,                     &
3003                                    init_3d%fill_chem(n),                      &
3004                                    .FALSE.,                                   &
3005                                    TRIM( init_3d%var_names_chem(n) ) )
3006                CALL get_attribute( id_dynamic, char_lod,                      &
3007                                    init_3d%lod_chem(n),                       &
3008                                    .FALSE.,                                   &
3009                                    TRIM( init_3d%var_names_chem(n) ) )
3010!
3011!--             Give message that only LOD=1 is allowed.
3012                IF ( init_3d%lod_chem(n) /= 1 )  THEN               
3013                   message_string = 'For chemistry variables only LOD=1 is ' //&
3014                                    'allowed.'
3015                   CALL message( 'netcdf_data_input_mod', 'PA0586',            &
3016                                 1, 2, 0, 6, 0 )
3017                ENDIF
3018!
3019!--             level-of-detail 1 - read initialization profile
3020                CALL get_variable( id_dynamic,                                 &
3021                                   TRIM( init_3d%var_names_chem(n) ),          &
3022                                   init_3d%chem_init(nzb+1:nzt,n) )
3023!
3024!--             Set bottom and top boundary condition (Neumann)
3025                init_3d%chem_init(nzb,n)   = init_3d%chem_init(nzb+1,n)
3026                init_3d%chem_init(nzt+1,n) = init_3d%chem_init(nzt,n)
3027               
3028                init_3d%from_file_chem(n) = .TRUE.
3029             ENDIF
3030          ENDDO
3031       ENDIF
3032!
3033!--    Close input file
3034       CALL close_input_file( id_dynamic )
3035#endif
3036!
3037!--    End of CPU measurement
3038       CALL cpu_log( log_point_s(85), 'NetCDF input init', 'stop' )
3039!
3040!--    Finally, check if the input data has any fill values. Please note,
3041!--    checks depend on the LOD of the input data.
3042       IF ( init_3d%from_file_u )  THEN
3043          check_passed = .TRUE.
3044          IF ( init_3d%lod_u == 1 )  THEN
3045             IF ( ANY( init_3d%u_init(nzb+1:nzt+1) == init_3d%fill_u ) )       &
3046                check_passed = .FALSE.
3047          ELSEIF ( init_3d%lod_u == 2 )  THEN
3048             IF ( ANY( u(nzb+1:nzt+1,nys:nyn,nxlu:nxr) == init_3d%fill_u ) )   &
3049                check_passed = .FALSE.
3050          ENDIF
3051          IF ( .NOT. check_passed )  THEN
3052             message_string = 'NetCDF input for init_atmosphere_u must ' //    &
3053                              'not contain any _FillValues'
3054             CALL message( 'netcdf_data_input_mod', 'PA0545', 2, 2, 0, 6, 0 )
3055          ENDIF
3056       ENDIF
3057
3058       IF ( init_3d%from_file_v )  THEN
3059          check_passed = .TRUE.
3060          IF ( init_3d%lod_v == 1 )  THEN
3061             IF ( ANY( init_3d%v_init(nzb+1:nzt+1) == init_3d%fill_v ) )       &
3062                check_passed = .FALSE.
3063          ELSEIF ( init_3d%lod_v == 2 )  THEN
3064             IF ( ANY( v(nzb+1:nzt+1,nysv:nyn,nxl:nxr) == init_3d%fill_v ) )   &
3065                check_passed = .FALSE.
3066          ENDIF
3067          IF ( .NOT. check_passed )  THEN
3068             message_string = 'NetCDF input for init_atmosphere_v must ' //    &
3069                              'not contain any _FillValues'
3070             CALL message( 'netcdf_data_input_mod', 'PA0545', 2, 2, 0, 6, 0 )
3071          ENDIF
3072       ENDIF
3073
3074       IF ( init_3d%from_file_w )  THEN
3075          check_passed = .TRUE.
3076          IF ( init_3d%lod_w == 1 )  THEN
3077             IF ( ANY( init_3d%w_init(nzb+1:nzt) == init_3d%fill_w ) )         &
3078                check_passed = .FALSE.
3079          ELSEIF ( init_3d%lod_w == 2 )  THEN
3080             IF ( ANY( w(nzb+1:nzt,nys:nyn,nxl:nxr) == init_3d%fill_w ) )      &
3081                check_passed = .FALSE.
3082          ENDIF
3083          IF ( .NOT. check_passed )  THEN
3084             message_string = 'NetCDF input for init_atmosphere_w must ' //    &
3085                              'not contain any _FillValues'
3086             CALL message( 'netcdf_data_input_mod', 'PA0545', 2, 2, 0, 6, 0 )
3087          ENDIF
3088       ENDIF
3089
3090       IF ( init_3d%from_file_pt )  THEN
3091          check_passed = .TRUE.
3092          IF ( init_3d%lod_pt == 1 )  THEN
3093             IF ( ANY( init_3d%pt_init(nzb+1:nzt+1) == init_3d%fill_pt ) )     &
3094                check_passed = .FALSE.
3095          ELSEIF ( init_3d%lod_pt == 2 )  THEN
3096             IF ( ANY( pt(nzb+1:nzt+1,nys:nyn,nxl:nxr) == init_3d%fill_pt ) )  &
3097                check_passed = .FALSE.
3098          ENDIF
3099          IF ( .NOT. check_passed )  THEN
3100             message_string = 'NetCDF input for init_atmosphere_pt must ' //   &
3101                              'not contain any _FillValues'
3102             CALL message( 'netcdf_data_input_mod', 'PA0545', 2, 2, 0, 6, 0 )
3103          ENDIF
3104       ENDIF
3105
3106       IF ( init_3d%from_file_q )  THEN
3107          check_passed = .TRUE.
3108          IF ( init_3d%lod_q == 1 )  THEN
3109             IF ( ANY( init_3d%q_init(nzb+1:nzt+1) == init_3d%fill_q ) )       &
3110                check_passed = .FALSE.
3111          ELSEIF ( init_3d%lod_q == 2 )  THEN
3112             IF ( ANY( q(nzb+1:nzt+1,nys:nyn,nxl:nxr) == init_3d%fill_q ) )    &
3113                check_passed = .FALSE.
3114          ENDIF
3115          IF ( .NOT. check_passed )  THEN
3116             message_string = 'NetCDF input for init_atmosphere_q must ' //    &
3117                              'not contain any _FillValues'
3118             CALL message( 'netcdf_data_input_mod', 'PA0545', 2, 2, 0, 6, 0 )
3119          ENDIF
3120       ENDIF
3121!
3122!--    Workaround for cyclic conditions. Please see above for further explanation.
3123       IF ( bc_lr_cyc  .AND.  nxl == 0 )  nxlu = nxl
3124       IF ( bc_ns_cyc  .AND.  nys == 0 )  nysv = nys
3125
3126    END SUBROUTINE netcdf_data_input_init_3d
3127
3128!------------------------------------------------------------------------------!
3129! Description:
3130! ------------
3131!> Checks input file for consistency and minimum requirements.
3132!------------------------------------------------------------------------------!
3133    SUBROUTINE netcdf_data_input_check_dynamic
3134
3135       USE control_parameters,                                                 &
3136           ONLY:  initializing_actions, message_string
3137
3138       IMPLICIT NONE
3139!
3140!--    Dynamic input file must also be present if initialization via inifor is
3141!--    prescribed.
3142       IF ( .NOT. input_pids_dynamic  .AND.                                    &
3143            INDEX( initializing_actions, 'inifor' ) /= 0 )  THEN
3144          message_string = 'initializing_actions = inifor requires dynamic ' //&
3145                           'input file ' // TRIM( input_file_dynamic ) //      &
3146                           TRIM( coupling_char )
3147          CALL message( 'netcdf_data_input_mod', 'PA0547', 1, 2, 0, 6, 0 )
3148       ENDIF
3149
3150    END SUBROUTINE netcdf_data_input_check_dynamic
3151
3152!------------------------------------------------------------------------------!
3153! Description:
3154! ------------
3155!> Checks input file for consistency and minimum requirements.
3156!------------------------------------------------------------------------------!
3157    SUBROUTINE netcdf_data_input_check_static
3158
3159       USE arrays_3d,                                                          &
3160           ONLY:  zu
3161
3162       USE control_parameters,                                                 &
3163           ONLY:  land_surface, message_string, urban_surface
3164
3165       USE indices,                                                            &
3166           ONLY:  nxl, nxr, nyn, nys, wall_flags_total_0
3167
3168       IMPLICIT NONE
3169
3170       INTEGER(iwp) ::  i      !< loop index along x-direction
3171       INTEGER(iwp) ::  j      !< loop index along y-direction
3172       INTEGER(iwp) ::  n_surf !< number of different surface types at given location
3173
3174       LOGICAL      ::  check_passed !< flag indicating if a check passed
3175
3176!
3177!--    Return if no static input file is available
3178       IF ( .NOT. input_pids_static )  RETURN
3179!
3180!--    Check for correct dimension of surface_fractions, should run from 0-2.
3181       IF ( surface_fraction_f%from_file )  THEN
3182          IF ( surface_fraction_f%nf-1 > 2 )  THEN
3183             message_string = 'nsurface_fraction must not be larger than 3.' 
3184             CALL message( 'netcdf_data_input_mod', 'PA0580', 1, 2, 0, 6, 0 )
3185          ENDIF
3186       ENDIF
3187!
3188!--    If 3D buildings are read, check if building information is consistent
3189!--    to numeric grid.
3190       IF ( buildings_f%from_file )  THEN
3191          IF ( buildings_f%lod == 2 )  THEN
3192             IF ( buildings_f%nz > SIZE( zu ) )  THEN
3193                message_string = 'Reading 3D building data - too much ' //     &
3194                                 'data points along the vertical coordinate.'
3195                CALL message( 'netcdf_data_input_mod', 'PA0552', 2, 2, 0, 6, 0 )
3196             ENDIF
3197
3198             IF ( ANY( ABS( buildings_f%z(0:buildings_f%nz-1) -                &
3199                       zu(0:buildings_f%nz-1) ) > 1E-6_wp ) )  THEN
3200                message_string = 'Reading 3D building data - vertical ' //     &
3201                                 'coordinate do not match numeric grid.'
3202                CALL message( 'netcdf_data_input_mod', 'PA0553', 2, 2, 0, 6, 0 )
3203             ENDIF
3204          ENDIF
3205       ENDIF
3206
3207!
3208!--    Skip further checks concerning buildings and natural surface properties
3209!--    if no urban surface and land surface model are applied.
3210       IF (  .NOT. land_surface  .AND.  .NOT. urban_surface )  RETURN
3211!
3212!--    Check for minimum requirement of surface-classification data in case
3213!--    static input file is used.
3214       IF ( ( .NOT. vegetation_type_f%from_file  .OR.                          &
3215              .NOT. pavement_type_f%from_file    .OR.                          &
3216              .NOT. water_type_f%from_file       .OR.                          &
3217              .NOT. soil_type_f%from_file             ) .OR.                   &
3218             ( urban_surface  .AND.  .NOT. building_type_f%from_file ) )  THEN
3219          message_string = 'Minimum requirement for surface classification ' //&
3220                           'is not fulfilled. At least ' //                    &
3221                           'vegetation_type, pavement_type, ' //               &
3222                           'soil_type and water_type are '//                   &
3223                           'required. If urban-surface model is applied, ' //  &
3224                           'also building_type is required'
3225          CALL message( 'netcdf_data_input_mod', 'PA0554', 1, 2, 0, 6, 0 )
3226       ENDIF
3227!
3228!--    Check for general availability of input variables.
3229!--    If vegetation_type is 0 at any location, vegetation_pars as well as
3230!--    root_area_dens_s are required.
3231       IF ( vegetation_type_f%from_file )  THEN
3232          IF ( ANY( vegetation_type_f%var == 0 ) )  THEN
3233             IF ( .NOT. vegetation_pars_f%from_file )  THEN
3234                message_string = 'If vegetation_type = 0 at any location, ' // &
3235                                 'vegetation_pars is required'
3236                CALL message( 'netcdf_data_input_mod', 'PA0555', 2, 2, myid, 6, 0 )
3237             ENDIF
3238             IF ( .NOT. root_area_density_lsm_f%from_file )  THEN
3239                message_string = 'If vegetation_type = 0 at any location, ' // &
3240                                 'root_area_dens_s is required'
3241                CALL message( 'netcdf_data_input_mod', 'PA0556', 2, 2, myid, 6, 0 )
3242             ENDIF
3243          ENDIF
3244       ENDIF
3245!
3246!--    If soil_type is zero at any location, soil_pars is required.
3247       IF ( soil_type_f%from_file )  THEN
3248          check_passed = .TRUE.
3249          IF ( ALLOCATED( soil_type_f%var_2d ) )  THEN
3250             IF ( ANY( soil_type_f%var_2d == 0 ) )  THEN
3251                IF ( .NOT. soil_pars_f%from_file )  check_passed = .FALSE.
3252             ENDIF
3253          ELSE
3254             IF ( ANY( soil_type_f%var_3d == 0 ) )  THEN
3255                IF ( .NOT. soil_pars_f%from_file )  check_passed = .FALSE.
3256             ENDIF
3257          ENDIF
3258          IF ( .NOT. check_passed )  THEN
3259             message_string = 'If soil_type = 0 at any location, ' //          &
3260                              'soil_pars is required'
3261             CALL message( 'netcdf_data_input_mod', 'PA0557', 2, 2, myid, 6, 0 )
3262          ENDIF
3263       ENDIF
3264!
3265!--    Buildings require a type in case of urban-surface model.
3266       IF ( buildings_f%from_file  .AND.  .NOT. building_type_f%from_file  )  THEN
3267          message_string = 'If buildings are provided, also building_type ' // &
3268                           'is required'
3269          CALL message( 'netcdf_data_input_mod', 'PA0581', 2, 2, 0, 6, 0 )
3270       ENDIF
3271!
3272!--    Buildings require an ID.
3273       IF ( buildings_f%from_file  .AND.  .NOT. building_id_f%from_file  )  THEN
3274          message_string = 'If buildings are provided, also building_id ' //   &
3275                           'is required'
3276          CALL message( 'netcdf_data_input_mod', 'PA0582', 2, 2, 0, 6, 0 )
3277       ENDIF
3278!
3279!--    If building_type is zero at any location, building_pars is required.
3280       IF ( building_type_f%from_file )  THEN
3281          IF ( ANY( building_type_f%var == 0 ) )  THEN
3282             IF ( .NOT. building_pars_f%from_file )  THEN
3283                message_string = 'If building_type = 0 at any location, ' //   &
3284                                 'building_pars is required'
3285                CALL message( 'netcdf_data_input_mod', 'PA0558', 2, 2, myid, 6, 0 )
3286             ENDIF
3287          ENDIF
3288       ENDIF
3289!
3290!--    If building_type is provided, also building_id is needed (due to the
3291!--    filtering algorithm).
3292       IF ( building_type_f%from_file  .AND.  .NOT. building_id_f%from_file )  &
3293       THEN
3294          message_string = 'If building_type is provided, also building_id '// &
3295                           'is required'
3296          CALL message( 'netcdf_data_input_mod', 'PA0519', 2, 2, 0, 6, 0 )
3297       ENDIF       
3298!
3299!--    If albedo_type is zero at any location, albedo_pars is required.
3300       IF ( albedo_type_f%from_file )  THEN
3301          IF ( ANY( albedo_type_f%var == 0 ) )  THEN
3302             IF ( .NOT. albedo_pars_f%from_file )  THEN
3303                message_string = 'If albedo_type = 0 at any location, ' //     &
3304                                 'albedo_pars is required'
3305                CALL message( 'netcdf_data_input_mod', 'PA0559', 2, 2, myid, 6, 0 )
3306             ENDIF
3307          ENDIF
3308       ENDIF
3309!
3310!--    If pavement_type is zero at any location, pavement_pars is required.
3311       IF ( pavement_type_f%from_file )  THEN
3312          IF ( ANY( pavement_type_f%var == 0 ) )  THEN
3313             IF ( .NOT. pavement_pars_f%from_file )  THEN
3314                message_string = 'If pavement_type = 0 at any location, ' //   &
3315                                 'pavement_pars is required'
3316                CALL message( 'netcdf_data_input_mod', 'PA0560', 2, 2, myid, 6, 0 )
3317             ENDIF
3318          ENDIF
3319       ENDIF
3320!
3321!--    If pavement_type is zero at any location, also pavement_subsurface_pars
3322!--    is required.
3323       IF ( pavement_type_f%from_file )  THEN
3324          IF ( ANY( pavement_type_f%var == 0 ) )  THEN
3325             IF ( .NOT. pavement_subsurface_pars_f%from_file )  THEN
3326                message_string = 'If pavement_type = 0 at any location, ' //   &
3327                                 'pavement_subsurface_pars is required'
3328                CALL message( 'netcdf_data_input_mod', 'PA0561', 2, 2, myid, 6, 0 )
3329             ENDIF
3330          ENDIF
3331       ENDIF
3332!
3333!--    If water_type is zero at any location, water_pars is required.
3334       IF ( water_type_f%from_file )  THEN
3335          IF ( ANY( water_type_f%var == 0 ) )  THEN
3336             IF ( .NOT. water_pars_f%from_file )  THEN
3337                message_string = 'If water_type = 0 at any location, ' //      &
3338                                 'water_pars is required'
3339                CALL message( 'netcdf_data_input_mod', 'PA0562', 2, 2,myid, 6, 0 )
3340             ENDIF
3341          ENDIF
3342       ENDIF
3343!
3344!--    Check for local consistency of the input data.
3345       DO  i = nxl, nxr
3346          DO  j = nys, nyn
3347!
3348!--          For each (y,x)-location at least one of the parameters
3349!--          vegetation_type, pavement_type, building_type, or water_type
3350!--          must be set to a non­missing value.
3351             IF ( land_surface  .AND.  .NOT. urban_surface )  THEN
3352                IF ( vegetation_type_f%var(j,i) == vegetation_type_f%fill  .AND.&
3353                     pavement_type_f%var(j,i)   == pavement_type_f%fill    .AND.&
3354                     water_type_f%var(j,i)      == water_type_f%fill )  THEN
3355                   WRITE( message_string, * )                                  &
3356                                    'At least one of the parameters '//        &
3357                                    'vegetation_type, pavement_type, '     //  &
3358                                    'or water_type must be set '//             &
3359                                    'to a non-missing value. Grid point: ', j, i
3360                   CALL message( 'netcdf_data_input_mod', 'PA0563', 2, 2, myid, 6, 0 )
3361                ENDIF
3362             ELSEIF ( land_surface  .AND.  urban_surface )  THEN
3363                IF ( vegetation_type_f%var(j,i) == vegetation_type_f%fill  .AND.&
3364                     pavement_type_f%var(j,i)   == pavement_type_f%fill    .AND.&
3365                     building_type_f%var(j,i)   == building_type_f%fill    .AND.&
3366                     water_type_f%var(j,i)      == water_type_f%fill )  THEN
3367                   WRITE( message_string, * )                                  &
3368                                 'At least one of the parameters '//           &
3369                                 'vegetation_type, pavement_type, '  //        &
3370                                 'building_type, or water_type must be set '// &
3371                                 'to a non-missing value. Grid point: ', j, i
3372                   CALL message( 'netcdf_data_input_mod', 'PA0563', 2, 2, myid, 6, 0 )
3373                ENDIF
3374             ENDIF
3375               
3376!
3377!--          Note that a soil_type is required for each location (y,x) where
3378!--          either vegetation_type or pavement_type is a non­missing value.
3379             IF ( ( vegetation_type_f%var(j,i) /= vegetation_type_f%fill  .OR. &
3380                    pavement_type_f%var(j,i)   /= pavement_type_f%fill ) )  THEN
3381                check_passed = .TRUE.
3382                IF ( ALLOCATED( soil_type_f%var_2d ) )  THEN
3383                   IF ( soil_type_f%var_2d(j,i) == soil_type_f%fill )          &
3384                      check_passed = .FALSE.
3385                ELSE
3386                   IF ( ANY( soil_type_f%var_3d(:,j,i) == soil_type_f%fill) )  &
3387                      check_passed = .FALSE.
3388                ENDIF
3389
3390                IF ( .NOT. check_passed )  THEN
3391                   message_string = 'soil_type is required for each '//        &
3392                                 'location (y,x) where vegetation_type or ' // &
3393                                 'pavement_type is a non-missing value.'
3394                   CALL message( 'netcdf_data_input_mod', 'PA0564',            &
3395                                  2, 2, myid, 6, 0 )
3396                ENDIF
3397             ENDIF
3398!
3399!--          Check for consistency of given types. At the moment, only one
3400!--          of vegetation, pavement, or water-type can be set. This is
3401!--          because no tile approach is yet implemented in the land-surface
3402!--          model. Later, when this is possible, surface fraction need to be
3403!--          given and the sum must not  be larger than 1. Please note, in case
3404!--          more than one type is given at a pixel, an error message will be
3405!--          given.
3406             n_surf = 0
3407             IF ( vegetation_type_f%var(j,i) /= vegetation_type_f%fill )       &
3408                n_surf = n_surf + 1
3409             IF ( water_type_f%var(j,i)      /= water_type_f%fill )            &
3410                n_surf = n_surf + 1
3411             IF ( pavement_type_f%var(j,i)   /= pavement_type_f%fill )         &
3412                n_surf = n_surf + 1
3413
3414             IF ( n_surf > 1 )  THEN
3415                WRITE( message_string, * )                                     &
3416                                 'More than one surface type (vegetation, '//  &
3417                                 'pavement, water) is given at a location. '// &
3418                                 'Please note, this is not possible at ' //    &
3419                                 'the moment as no tile approach has been ' // &
3420                                 'yet implemented. (i,j) = ', i, j
3421                CALL message( 'netcdf_data_input_mod', 'PA0565',               &
3422                               2, 2, myid, 6, 0 )
3423
3424!                 IF ( .NOT. surface_fraction_f%from_file )  THEN
3425!                    message_string = 'More than one surface type (vegetation '//&
3426!                                  'pavement, water) is given at a location. '// &
3427!                                  'Please note, this is not possible at ' //    &
3428!                                  'the moment as no tile approach is yet ' //   &
3429!                                  'implemented.'
3430!                    message_string = 'If more than one surface type is ' //     &
3431!                                  'given at a location, surface_fraction ' //   &
3432!                                  'must be provided.'
3433!                    CALL message( 'netcdf_data_input_mod', 'PA0565',            &
3434!                                   2, 2, myid, 6, 0 )
3435!                 ELSEIF ( ANY ( surface_fraction_f%frac(:,j,i) ==               &
3436!                                surface_fraction_f%fill ) )  THEN
3437!                    message_string = 'If more than one surface type is ' //     &
3438!                                  'given at a location, surface_fraction ' //   &
3439!                                  'must be provided.'
3440!                    CALL message( 'netcdf_data_input_mod', 'PA0565',            &
3441!                                   2, 2, myid, 6, 0 )
3442!                 ENDIF
3443             ENDIF
3444!
3445!--          Check for further mismatches. e.g. relative fractions exceed 1 or
3446!--          vegetation_type is set but surface vegetation fraction is zero,
3447!--          etc..
3448             IF ( surface_fraction_f%from_file )  THEN
3449!
3450!--             If surface fractions is given, also check that only one type
3451!--             is given.
3452                IF ( SUM( MERGE( 1, 0, surface_fraction_f%frac(:,j,i) /= 0.0_wp&
3453                                .AND.  surface_fraction_f%frac(:,j,i) /=       &
3454                                       surface_fraction_f%fill  ) ) > 1 )  THEN
3455                   WRITE( message_string, * )                                  &
3456                                    'surface_fraction is given for more ' //   &
3457                                    'than one type. ' //                       &
3458                                    'Please note, this is not possible at ' // &
3459                                    'the moment as no tile approach has '//    &
3460                                    'yet been implemented. (i, j) = ', i, j
3461                   CALL message( 'netcdf_data_input_mod', 'PA0676',            &
3462                                  2, 2, myid, 6, 0 )
3463                ENDIF
3464!
3465!--             Sum of relative fractions must be 1. Note, attributed to type
3466!--             conversions due to reading, the sum of surface fractions
3467!--             might be not exactly 1. Hence, the sum is check with a
3468!--             tolerance. Later, in the land-surface model, the relative
3469!--             fractions are normalized to one. Actually, surface fractions
3470!--             shall be _FillValue at building grid points, however, in order
3471!--             to relax this requirement and allow that surface-fraction can
3472!--             also be zero at these grid points, only perform this check
3473!--             at locations where some vegetation, pavement or water is defined.
3474                IF ( vegetation_type_f%var(j,i) /= vegetation_type_f%fill  .OR.&
3475                     pavement_type_f%var(j,i)   /= pavement_type_f%fill    .OR.&
3476                     water_type_f%var(j,i)      /= water_type_f%fill )  THEN
3477                   IF ( SUM ( surface_fraction_f%frac(0:2,j,i) ) >             &
3478                        1.0_wp + 1E-8_wp  .OR.                                 &
3479                        SUM ( surface_fraction_f%frac(0:2,j,i) ) <             &
3480                        1.0_wp - 1E-8_wp )  THEN
3481                      WRITE( message_string, * )                               &
3482                                    'The sum of all land-surface fractions ' //&
3483                                    'must equal 1. (i, j) = ', i, j
3484                      CALL message( 'netcdf_data_input_mod', 'PA0566',         &
3485                                     2, 2, myid, 6, 0 )
3486                   ENDIF
3487                ENDIF
3488!
3489!--             Relative fraction for a type must not be zero at locations where
3490!--             this type is set.
3491                IF (                                                           &
3492                  ( vegetation_type_f%var(j,i) /= vegetation_type_f%fill  .AND.&
3493                 ( surface_fraction_f%frac(ind_veg_wall,j,i) == 0.0_wp .OR.    &
3494                   surface_fraction_f%frac(ind_veg_wall,j,i) ==                &
3495                                                     surface_fraction_f%fill ) &
3496                  )  .OR.                                                      &
3497                  ( pavement_type_f%var(j,i) /= pavement_type_f%fill     .AND. &
3498                 ( surface_fraction_f%frac(ind_pav_green,j,i) == 0.0_wp .OR.   &
3499                   surface_fraction_f%frac(ind_pav_green,j,i) ==               &
3500                                                     surface_fraction_f%fill ) &
3501                  )  .OR.                                                      &
3502                  ( water_type_f%var(j,i) /= water_type_f%fill           .AND. &
3503                 ( surface_fraction_f%frac(ind_wat_win,j,i) == 0.0_wp .OR.     &
3504                   surface_fraction_f%frac(ind_wat_win,j,i) ==                 &
3505                                                     surface_fraction_f%fill ) &
3506                  ) )  THEN
3507                   WRITE( message_string, * ) 'Mismatch in setting of '     // &
3508                             'surface_fraction. Vegetation-, pavement-, or '// &
3509                             'water surface is given at (i,j) = ( ', i, j,     &
3510                             ' ), but surface fraction is 0 for the given type.'
3511                   CALL message( 'netcdf_data_input_mod', 'PA0567',            &
3512                                  2, 2, myid, 6, 0 )
3513                ENDIF
3514!
3515!--             Relative fraction for a type must not contain non-zero values
3516!--             if this type is not set.
3517                IF (                                                           &
3518                  ( vegetation_type_f%var(j,i) == vegetation_type_f%fill  .AND.&
3519                 ( surface_fraction_f%frac(ind_veg_wall,j,i) /= 0.0_wp .AND.   &
3520                   surface_fraction_f%frac(ind_veg_wall,j,i) /=                &
3521                                                     surface_fraction_f%fill ) &
3522                  )  .OR.                                                      &
3523                  ( pavement_type_f%var(j,i) == pavement_type_f%fill     .AND. &
3524                 ( surface_fraction_f%frac(ind_pav_green,j,i) /= 0.0_wp .AND.  &
3525                   surface_fraction_f%frac(ind_pav_green,j,i) /=               &
3526                                                     surface_fraction_f%fill ) &
3527                  )  .OR.                                                      &
3528                  ( water_type_f%var(j,i) == water_type_f%fill           .AND. &
3529                 ( surface_fraction_f%frac(ind_wat_win,j,i) /= 0.0_wp .AND.    &
3530                   surface_fraction_f%frac(ind_wat_win,j,i) /=                 &
3531                                                     surface_fraction_f%fill ) &
3532                  ) )  THEN
3533                   WRITE( message_string, * ) 'Mismatch in setting of '     // &
3534                             'surface_fraction. Vegetation-, pavement-, or '// &
3535                             'water surface is not given at (i,j) = ( ', i, j, &
3536                             ' ), but surface fraction is not 0 for the ' //   &
3537                             'given type.'
3538                   CALL message( 'netcdf_data_input_mod', 'PA0568',            &
3539                                  2, 2, myid, 6, 0 )
3540                ENDIF
3541             ENDIF
3542!
3543!--          Check vegetation_pars. If vegetation_type is 0, all parameters
3544!--          need to be set, otherwise, single parameters set by
3545!--          vegetation_type can be overwritten.
3546             IF ( vegetation_type_f%from_file )  THEN
3547                IF ( vegetation_type_f%var(j,i) == 0 )  THEN
3548                   IF ( ANY( vegetation_pars_f%pars_xy(:,j,i) ==               &
3549                             vegetation_pars_f%fill ) )  THEN
3550                      message_string = 'If vegetation_type(y,x) = 0, all '  // &
3551                                       'parameters of vegetation_pars at '//   &
3552                                       'this location must be set.'
3553                      CALL message( 'netcdf_data_input_mod', 'PA0569',         &
3554                                     2, 2, myid, 6, 0 )
3555                   ENDIF
3556                ENDIF
3557             ENDIF
3558!
3559!--          Check root distribution. If vegetation_type is 0, all levels must
3560!--          be set.
3561             IF ( vegetation_type_f%from_file )  THEN
3562                IF ( vegetation_type_f%var(j,i) == 0 )  THEN
3563                   IF ( ANY( root_area_density_lsm_f%var(:,j,i) ==             &
3564                             root_area_density_lsm_f%fill ) )  THEN
3565                      message_string = 'If vegetation_type(y,x) = 0, all ' //  &
3566                                       'levels of root_area_dens_s ' //        &
3567                                       'must be set at this location.'
3568                      CALL message( 'netcdf_data_input_mod', 'PA0570',         &
3569                                     2, 2, myid, 6, 0 )
3570                   ENDIF
3571                ENDIF
3572             ENDIF
3573!
3574!--          Check soil parameters. If soil_type is 0, all parameters
3575!--          must be set.
3576             IF ( soil_type_f%from_file )  THEN
3577                check_passed = .TRUE.
3578                IF ( ALLOCATED( soil_type_f%var_2d ) )  THEN
3579                   IF ( soil_type_f%var_2d(j,i) == 0 )  THEN
3580                      IF ( ANY( soil_pars_f%pars_xy(:,j,i) ==                  &
3581                                soil_pars_f%fill ) )  check_passed = .FALSE.
3582                   ENDIF
3583                ELSE
3584                   IF ( ANY( soil_type_f%var_3d(:,j,i) == 0 ) )  THEN
3585                      IF ( ANY( soil_pars_f%pars_xy(:,j,i) ==                  &
3586                                soil_pars_f%fill ) )  check_passed = .FALSE.
3587                   ENDIF
3588                ENDIF
3589                IF ( .NOT. check_passed )  THEN
3590                   message_string = 'If soil_type(y,x) = 0, all levels of '  //&
3591                                    'soil_pars at this location must be set.'
3592                   CALL message( 'netcdf_data_input_mod', 'PA0571',            &
3593                                  2, 2, myid, 6, 0 )
3594                ENDIF
3595             ENDIF
3596
3597!
3598!--          Check building parameters. If building_type is 0, all parameters
3599!--          must be set.
3600             IF ( building_type_f%from_file )  THEN
3601                IF ( building_type_f%var(j,i) == 0 )  THEN
3602                   IF ( ANY( building_pars_f%pars_xy(:,j,i) ==                 &
3603                             building_pars_f%fill ) )  THEN
3604                      message_string = 'If building_type(y,x) = 0, all ' //    &
3605                                       'parameters of building_pars at this '//&
3606                                       'location must be set.'
3607                      CALL message( 'netcdf_data_input_mod', 'PA0572',         &
3608                                     2, 2, myid, 6, 0 )
3609                   ENDIF
3610                ENDIF
3611             ENDIF
3612!
3613!--          Check if building_type is set at each building and vice versa.
3614!--          Please note, buildings are already processed and filtered.
3615!--          For this reason, consistency checks are based on wall_flags_total_0
3616!--          rather than buildings_f (buildings are represented by bit 6 in
3617!--          wall_flags_total_0).
3618             IF ( building_type_f%from_file  .AND.  buildings_f%from_file )  THEN
3619                IF ( ANY( BTEST ( wall_flags_total_0(:,j,i), 6 ) )  .AND.      &
3620                     building_type_f%var(j,i) == building_type_f%fill  .OR.    &
3621               .NOT. ANY( BTEST ( wall_flags_total_0(:,j,i), 6 ) )  .AND.      &
3622                     building_type_f%var(j,i) /= building_type_f%fill )  THEN
3623                   WRITE( message_string, * ) 'Each location where a ' //      &
3624                                   'building is set requires a type ' //       &
3625                                   '( and vice versa ) in case the ' //        &
3626                                   'urban-surface model is applied. ' //       &
3627                                   'i, j = ', i, j
3628                   CALL message( 'netcdf_data_input_mod', 'PA0573',            &
3629                                  2, 2, myid, 6, 0 )
3630                ENDIF
3631             ENDIF
3632!
3633!--          Check if at each location where a building is present also an ID
3634!--          is set and vice versa.
3635             IF ( buildings_f%from_file )  THEN
3636                IF ( ANY( BTEST ( wall_flags_total_0(:,j,i), 6 ) )  .AND.     &
3637                     building_id_f%var(j,i) == building_id_f%fill  .OR.       &
3638               .NOT. ANY( BTEST ( wall_flags_total_0(:,j,i), 6 ) )  .AND.     &
3639                     building_id_f%var(j,i) /= building_id_f%fill )  THEN
3640                   WRITE( message_string, * ) 'Each location where a ' //     &
3641                                   'building is set requires an ID ' //       &
3642                                   '( and vice versa ). i, j = ', i, j
3643                   CALL message( 'netcdf_data_input_mod', 'PA0574',           &
3644                                  2, 2, myid, 6, 0 )
3645                ENDIF
3646             ENDIF
3647!
3648!--          Check if building ID is set where a bulding is defined.
3649             IF ( buildings_f%from_file )  THEN
3650                IF ( ANY( BTEST ( wall_flags_total_0(:,j,i), 6 ) )  .AND.     &
3651                     building_id_f%var(j,i) == building_id_f%fill )  THEN
3652                   WRITE( message_string, * ) 'Each building grid point '//   &
3653                                              'requires an ID.', i, j
3654                   CALL message( 'netcdf_data_input_mod', 'PA0575',           &
3655                                  2, 2, myid, 6, 0 )
3656                ENDIF
3657             ENDIF
3658!
3659!--          Check albedo parameters. If albedo_type is 0, all parameters
3660!--          must be set.
3661             IF ( albedo_type_f%from_file )  THEN
3662                IF ( albedo_type_f%var(j,i) == 0 )  THEN
3663                   IF ( ANY( albedo_pars_f%pars_xy(:,j,i) ==                   &
3664                             albedo_pars_f%fill ) )  THEN
3665                      message_string = 'If albedo_type(y,x) = 0, all ' //      &
3666                                       'parameters of albedo_pars at this ' // &
3667                                       'location must be set.'
3668                      CALL message( 'netcdf_data_input_mod', 'PA0576',         &
3669                                     2, 2, myid, 6, 0 )
3670                   ENDIF
3671                ENDIF
3672             ENDIF
3673
3674!
3675!--          Check pavement parameters. If pavement_type is 0, all parameters
3676!--          of pavement_pars must be set at this location.
3677             IF ( pavement_type_f%from_file )  THEN
3678                IF ( pavement_type_f%var(j,i) == 0 )  THEN
3679                   IF ( ANY( pavement_pars_f%pars_xy(:,j,i) ==                 &
3680                             pavement_pars_f%fill ) )  THEN
3681                      message_string = 'If pavement_type(y,x) = 0, all ' //    &
3682                                       'parameters of pavement_pars at this '//&
3683                                       'location must be set.'
3684                      CALL message( 'netcdf_data_input_mod', 'PA0577',         &
3685                                     2, 2, myid, 6, 0 )
3686                   ENDIF
3687                ENDIF
3688             ENDIF
3689!
3690!--          Check pavement-subsurface parameters. If pavement_type is 0,
3691!--          all parameters of pavement_subsurface_pars must be set  at this
3692!--          location.
3693             IF ( pavement_type_f%from_file )  THEN
3694                IF ( pavement_type_f%var(j,i) == 0 )  THEN
3695                   IF ( ANY( pavement_subsurface_pars_f%pars_xyz(:,:,j,i) ==   &
3696                             pavement_subsurface_pars_f%fill ) )  THEN
3697                      message_string = 'If pavement_type(y,x) = 0, all ' //    &
3698                                       'parameters of '                  //    &
3699                                       'pavement_subsurface_pars at this '//   &
3700                                       'location must be set.'
3701                      CALL message( 'netcdf_data_input_mod', 'PA0578',         &
3702                                     2, 2, myid, 6, 0 )
3703                   ENDIF
3704                ENDIF
3705             ENDIF
3706
3707!
3708!--          Check water parameters. If water_type is 0, all parameters
3709!--          must be set  at this location.
3710             IF ( water_type_f%from_file )  THEN
3711                IF ( water_type_f%var(j,i) == 0 )  THEN
3712                   IF ( ANY( water_pars_f%pars_xy(:,j,i) ==                    &
3713                             water_pars_f%fill ) )  THEN
3714                      message_string = 'If water_type(y,x) = 0, all ' //       &
3715                                       'parameters of water_pars at this ' //  &
3716                                       'location must be set.'
3717                      CALL message( 'netcdf_data_input_mod', 'PA0579',         &
3718                                     2, 2, myid, 6, 0 )
3719                   ENDIF
3720                ENDIF
3721             ENDIF
3722
3723          ENDDO
3724       ENDDO
3725
3726    END SUBROUTINE netcdf_data_input_check_static
3727
3728!------------------------------------------------------------------------------!
3729! Description:
3730! ------------
3731!> Resize 8-bit 2D Integer array: (nys:nyn,nxl:nxr) -> (nysg:nyng,nxlg:nxrg)
3732!------------------------------------------------------------------------------!
3733    SUBROUTINE resize_array_2d_int8( var, js, je, is, ie )
3734   
3735       IMPLICIT NONE
3736
3737       INTEGER(iwp) ::  je !< upper index bound along y direction
3738       INTEGER(iwp) ::  js !< lower index bound along y direction
3739       INTEGER(iwp) ::  ie !< upper index bound along x direction
3740       INTEGER(iwp) ::  is !< lower index bound along x direction
3741       
3742       INTEGER(KIND=1), DIMENSION(:,:), ALLOCATABLE ::  var     !< treated variable
3743       INTEGER(KIND=1), DIMENSION(:,:), ALLOCATABLE ::  var_tmp !< temporary copy
3744!
3745!--    Allocate temporary variable
3746       ALLOCATE( var_tmp(js-nbgp:je+nbgp,is-nbgp:ie+nbgp) )
3747!
3748!--    Temporary copy of the variable
3749       var_tmp(js:je,is:ie) = var(js:je,is:ie)
3750!
3751!--    Resize the array
3752       DEALLOCATE( var )
3753       ALLOCATE( var(js-nbgp:je+nbgp,is-nbgp:ie+nbgp) )
3754!
3755!--    Transfer temporary copy back to original array
3756       var(js:je,is:ie) = var_tmp(js:je,is:ie)
3757
3758    END SUBROUTINE resize_array_2d_int8
3759   
3760!------------------------------------------------------------------------------!
3761! Description:
3762! ------------
3763!> Resize 32-bit 2D Integer array: (nys:nyn,nxl:nxr) -> (nysg:nyng,nxlg:nxrg)
3764!------------------------------------------------------------------------------!
3765    SUBROUTINE resize_array_2d_int32( var, js, je, is, ie )
3766
3767       IMPLICIT NONE
3768       
3769       INTEGER(iwp) ::  je !< upper index bound along y direction
3770       INTEGER(iwp) ::  js !< lower index bound along y direction
3771       INTEGER(iwp) ::  ie !< upper index bound along x direction
3772       INTEGER(iwp) ::  is !< lower index bound along x direction
3773
3774       INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE ::  var     !< treated variable
3775       INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE ::  var_tmp !< temporary copy
3776!
3777!--    Allocate temporary variable
3778       ALLOCATE( var_tmp(js-nbgp:je+nbgp,is-nbgp:ie+nbgp) )
3779!
3780!--    Temporary copy of the variable
3781       var_tmp(js:je,is:ie) = var(js:je,is:ie)
3782!
3783!--    Resize the array
3784       DEALLOCATE( var )
3785       ALLOCATE( var(js-nbgp:je+nbgp,is-nbgp:ie+nbgp) )
3786!
3787!--    Transfer temporary copy back to original array
3788       var(js:je,is:ie) = var_tmp(js:je,is:ie)
3789
3790    END SUBROUTINE resize_array_2d_int32
3791   
3792!------------------------------------------------------------------------------!
3793! Description:
3794! ------------
3795!> Resize 8-bit 3D Integer array: (:,nys:nyn,nxl:nxr) -> (:,nysg:nyng,nxlg:nxrg)
3796!------------------------------------------------------------------------------!
3797    SUBROUTINE resize_array_3d_int8( var, ks, ke, js, je, is, ie )
3798
3799       IMPLICIT NONE
3800
3801       INTEGER(iwp) ::  je !< upper index bound along y direction
3802       INTEGER(iwp) ::  js !< lower index bound along y direction
3803       INTEGER(iwp) ::  ie !< upper index bound along x direction
3804       INTEGER(iwp) ::  is !< lower index bound along x direction
3805       INTEGER(iwp) ::  ke !< upper bound of treated array in z-direction 
3806       INTEGER(iwp) ::  ks !< lower bound of treated array in z-direction 
3807       
3808       INTEGER(KIND=1), DIMENSION(:,:,:), ALLOCATABLE ::  var     !< treated variable
3809       INTEGER(KIND=1), DIMENSION(:,:,:), ALLOCATABLE ::  var_tmp !< temporary copy
3810!
3811!--    Allocate temporary variable
3812       ALLOCATE( var_tmp(ks:ke,js-nbgp:je+nbgp,is-nbgp:ie+nbgp) )
3813!
3814!--    Temporary copy of the variable
3815       var_tmp(ks:ke,js:je,is:ie) = var(ks:ke,js:je,is:ie)
3816!
3817!--    Resize the array
3818       DEALLOCATE( var )
3819       ALLOCATE( var(ks:ke,js-nbgp:je+nbgp,is-nbgp:ie+nbgp) )
3820!
3821!--    Transfer temporary copy back to original array
3822       var(ks:ke,js:je,is:ie) = var_tmp(ks:ke,js:je,is:ie)
3823
3824    END SUBROUTINE resize_array_3d_int8
3825   
3826!------------------------------------------------------------------------------!
3827! Description:
3828! ------------
3829!> Resize 3D Real array: (:,nys:nyn,nxl:nxr) -> (:,nysg:nyng,nxlg:nxrg)
3830!------------------------------------------------------------------------------!
3831    SUBROUTINE resize_array_3d_real( var, ks, ke, js, je, is, ie )
3832
3833       IMPLICIT NONE
3834
3835       INTEGER(iwp) ::  je !< upper index bound along y direction
3836       INTEGER(iwp) ::  js !< lower index bound along y direction
3837       INTEGER(iwp) ::  ie !< upper index bound along x direction
3838       INTEGER(iwp) ::  is !< lower index bound along x direction
3839       INTEGER(iwp) ::  ke !< upper bound of treated array in z-direction 
3840       INTEGER(iwp) ::  ks !< lower bound of treated array in z-direction 
3841       
3842       REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  var     !< treated variable
3843       REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  var_tmp !< temporary copy
3844!
3845!--    Allocate temporary variable
3846       ALLOCATE( var_tmp(ks:ke,js-nbgp:je+nbgp,is-nbgp:ie+nbgp) )
3847!
3848!--    Temporary copy of the variable
3849       var_tmp(ks:ke,js:je,is:ie) = var(ks:ke,js:je,is:ie)
3850!
3851!--    Resize the array
3852       DEALLOCATE( var )
3853       ALLOCATE( var(ks:ke,js-nbgp:je+nbgp,is-nbgp:ie+nbgp) )
3854!
3855!--    Transfer temporary copy back to original array
3856       var(ks:ke,js:je,is:ie) = var_tmp(ks:ke,js:je,is:ie)
3857
3858    END SUBROUTINE resize_array_3d_real
3859   
3860!------------------------------------------------------------------------------!
3861! Description:
3862! ------------
3863!> Resize 4D Real array: (:,:,nys:nyn,nxl:nxr) -> (:,nysg:nyng,nxlg:nxrg)
3864!------------------------------------------------------------------------------!
3865    SUBROUTINE resize_array_4d_real( var, k1s, k1e, k2s, k2e, js, je, is, ie )
3866
3867       IMPLICIT NONE
3868       
3869       INTEGER(iwp) ::  je  !< upper index bound along y direction
3870       INTEGER(iwp) ::  js  !< lower index bound along y direction
3871       INTEGER(iwp) ::  ie  !< upper index bound along x direction
3872       INTEGER(iwp) ::  is  !< lower index bound along x direction
3873       INTEGER(iwp) ::  k1e !< upper bound of treated array in z-direction 
3874       INTEGER(iwp) ::  k1s !< lower bound of treated array in z-direction
3875       INTEGER(iwp) ::  k2e !< upper bound of treated array along parameter space 
3876       INTEGER(iwp) ::  k2s !< lower bound of treated array along parameter space 
3877       
3878       REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE ::  var     !< treated variable
3879       REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE ::  var_tmp !< temporary copy
3880!
3881!--    Allocate temporary variable
3882       ALLOCATE( var_tmp(k1s:k1e,k2s:k2e,js-nbgp:je+nbgp,is-nbgp:ie+nbgp) )
3883!
3884!--    Temporary copy of the variable
3885       var_tmp(k1s:k1e,k2s:k2e,js:je,is:ie) = var(k1s:k1e,k2s:k2e,js:je,is:ie)
3886!
3887!--    Resize the array
3888       DEALLOCATE( var )
3889       ALLOCATE( var(k1s:k1e,k2s:k2e,js-nbgp:je+nbgp,is-nbgp:ie+nbgp) )
3890!
3891!--    Transfer temporary copy back to original array
3892       var(k1s:k1e,k2s:k2e,js:je,is:ie) = var_tmp(k1s:k1e,k2s:k2e,js:je,is:ie)
3893
3894    END SUBROUTINE resize_array_4d_real
3895
3896!------------------------------------------------------------------------------!
3897! Description:
3898! ------------
3899!> Checks if a given variables is on file
3900!------------------------------------------------------------------------------!
3901    FUNCTION check_existence( vars_in_file, var_name )
3902
3903       IMPLICIT NONE
3904
3905       CHARACTER(LEN=*) ::  var_name                   !< variable to be checked
3906       CHARACTER(LEN=*), DIMENSION(:) ::  vars_in_file !< list of variables in file
3907
3908       INTEGER(iwp) ::  i                              !< loop variable
3909
3910       LOGICAL ::  check_existence                     !< flag indicating whether a variable exist or not - actual return value
3911
3912       i = 1
3913       check_existence = .FALSE.
3914       DO  WHILE ( i <= SIZE( vars_in_file ) )
3915          check_existence = TRIM( vars_in_file(i) ) == TRIM( var_name )  .OR.  &
3916                            check_existence
3917          i = i + 1
3918       ENDDO
3919
3920       RETURN
3921
3922    END FUNCTION check_existence
3923
3924
3925!------------------------------------------------------------------------------!
3926! Description:
3927! ------------
3928!> Closes an existing netCDF file.
3929!------------------------------------------------------------------------------!
3930    SUBROUTINE close_input_file( id )
3931
3932       USE pegrid
3933
3934       IMPLICIT NONE
3935
3936       INTEGER(iwp), INTENT(INOUT)        ::  id        !< file id
3937
3938#if defined( __netcdf )
3939       nc_stat = NF90_CLOSE( id )
3940       CALL handle_error( 'close', 540 )
3941#endif
3942    END SUBROUTINE close_input_file
3943
3944!------------------------------------------------------------------------------!
3945! Description:
3946! ------------
3947!> Opens an existing netCDF file for reading only and returns its id.
3948!------------------------------------------------------------------------------!
3949    SUBROUTINE open_read_file( filename, id )
3950
3951       USE pegrid
3952
3953       IMPLICIT NONE
3954
3955       CHARACTER (LEN=*), INTENT(IN) ::  filename  !< filename
3956       INTEGER(iwp), INTENT(INOUT)   ::  id        !< file id
3957
3958#if defined( __netcdf )
3959
3960#if defined( __netcdf4_parallel )
3961!
3962!--    If __netcdf4_parallel is defined, parrallel NetCDF will be used.
3963       nc_stat = NF90_OPEN( filename, IOR( NF90_NOWRITE, NF90_MPIIO ), id,     &
3964                            COMM = comm2d, INFO = MPI_INFO_NULL )
3965!
3966!--    In case the previous open call fails, check for possible Netcdf 3 file,
3967!--    and open it. However, this case, disable parallel access.
3968       IF( nc_stat /= NF90_NOERR )  THEN
3969          nc_stat = NF90_OPEN( filename, NF90_NOWRITE, id )
3970          collective_read = .FALSE.
3971       ELSE
3972#if defined( __nec )
3973          collective_read = .FALSE.   ! collective read causes hang situations on NEC Aurora
3974#else
3975          collective_read = .TRUE.
3976#endif
3977       ENDIF
3978#else
3979!
3980!--    All MPI processes open the file and read it (but not in parallel).
3981       nc_stat = NF90_OPEN( filename, NF90_NOWRITE, id )
3982#endif
3983
3984       CALL handle_error( 'open_read_file', 539 )
3985
3986#endif
3987    END SUBROUTINE open_read_file
3988
3989!------------------------------------------------------------------------------!
3990! Description:
3991! ------------
3992!> Reads global or variable-related attributes of type INTEGER (32-bit)
3993!------------------------------------------------------------------------------!
3994     SUBROUTINE get_attribute_int32( id, attribute_name, value, global,        &
3995                                     variable_name )
3996
3997       USE pegrid
3998
3999       IMPLICIT NONE
4000
4001       CHARACTER(LEN=*)            ::  attribute_name   !< attribute name
4002       CHARACTER(LEN=*), OPTIONAL  ::  variable_name    !< variable name
4003
4004       INTEGER(iwp), INTENT(IN)    ::  id               !< file id
4005       INTEGER(iwp)                ::  id_var           !< variable id
4006       INTEGER(iwp), INTENT(INOUT) ::  value            !< read value
4007
4008       LOGICAL, INTENT(IN) ::  global                   !< flag indicating global attribute
4009#if defined( __netcdf )
4010
4011!
4012!--    Read global attribute
4013       IF ( global )  THEN
4014          nc_stat = NF90_GET_ATT( id, NF90_GLOBAL, TRIM( attribute_name ), value )
4015          CALL handle_error( 'get_attribute_int32 global', 522, attribute_name )
4016!
4017!--    Read attributes referring to a single variable. Therefore, first inquire
4018!--    variable id
4019       ELSE
4020          nc_stat = NF90_INQ_VARID( id, TRIM( variable_name ), id_var )
4021          CALL handle_error( 'get_attribute_int32', 522, attribute_name )
4022          nc_stat = NF90_GET_ATT( id, id_var, TRIM( attribute_name ), value )
4023          CALL handle_error( 'get_attribute_int32', 522, attribute_name )
4024       ENDIF
4025#endif
4026    END SUBROUTINE get_attribute_int32
4027
4028!------------------------------------------------------------------------------!
4029! Description:
4030! ------------
4031!> Reads global or variable-related attributes of type INTEGER (8-bit)
4032!------------------------------------------------------------------------------!
4033     SUBROUTINE get_attribute_int8( id, attribute_name, value, global,         &
4034                                    variable_name )
4035
4036       USE pegrid
4037
4038       IMPLICIT NONE
4039
4040       CHARACTER(LEN=*)            ::  attribute_name   !< attribute name
4041       CHARACTER(LEN=*), OPTIONAL  ::  variable_name    !< variable name
4042
4043       INTEGER(iwp), INTENT(IN)    ::  id               !< file id
4044       INTEGER(iwp)                ::  id_var           !< variable id
4045       INTEGER(KIND=1), INTENT(INOUT) ::  value         !< read value
4046
4047       LOGICAL, INTENT(IN) ::  global                   !< flag indicating global attribute
4048#if defined( __netcdf )
4049
4050!
4051!--    Read global attribute
4052       IF ( global )  THEN
4053          nc_stat = NF90_GET_ATT( id, NF90_GLOBAL, TRIM( attribute_name ), value )
4054          CALL handle_error( 'get_attribute_int8 global', 523, attribute_name )
4055!
4056!--    Read attributes referring to a single variable. Therefore, first inquire
4057!--    variable id
4058       ELSE
4059          nc_stat = NF90_INQ_VARID( id, TRIM( variable_name ), id_var )
4060          CALL handle_error( 'get_attribute_int8', 523, attribute_name )
4061          nc_stat = NF90_GET_ATT( id, id_var, TRIM( attribute_name ), value )
4062          CALL handle_error( 'get_attribute_int8', 523, attribute_name )
4063       ENDIF
4064#endif
4065    END SUBROUTINE get_attribute_int8
4066
4067!------------------------------------------------------------------------------!
4068! Description:
4069! ------------
4070!> Reads global or variable-related attributes of type REAL
4071!------------------------------------------------------------------------------!
4072     SUBROUTINE get_attribute_real( id, attribute_name, value, global,         &
4073                                    variable_name )
4074
4075       USE pegrid
4076
4077       IMPLICIT NONE
4078
4079       CHARACTER(LEN=*)            ::  attribute_name   !< attribute name
4080       CHARACTER(LEN=*), OPTIONAL  ::  variable_name    !< variable name
4081
4082       INTEGER(iwp), INTENT(IN)    ::  id               !< file id
4083       INTEGER(iwp)                ::  id_var           !< variable id
4084
4085       LOGICAL, INTENT(IN) ::  global                   !< flag indicating global attribute
4086
4087       REAL(wp), INTENT(INOUT)     ::  value            !< read value
4088#if defined( __netcdf )
4089
4090
4091!
4092!-- Read global attribute
4093       IF ( global )  THEN
4094          nc_stat = NF90_GET_ATT( id, NF90_GLOBAL, TRIM( attribute_name ), value )
4095          CALL handle_error( 'get_attribute_real global', 524, attribute_name )
4096!
4097!-- Read attributes referring to a single variable. Therefore, first inquire
4098!-- variable id
4099       ELSE
4100          nc_stat = NF90_INQ_VARID( id, TRIM( variable_name ), id_var )
4101          CALL handle_error( 'get_attribute_real', 524, attribute_name )
4102          nc_stat = NF90_GET_ATT( id, id_var, TRIM( attribute_name ), value )
4103          CALL handle_error( 'get_attribute_real', 524, attribute_name )
4104       ENDIF
4105#endif
4106    END SUBROUTINE get_attribute_real
4107
4108!------------------------------------------------------------------------------!
4109! Description:
4110! ------------
4111!> Reads global or variable-related attributes of type CHARACTER
4112!> Remark: reading attributes of type NF_STRING return an error code 56 -
4113!> Attempt to convert between text & numbers.
4114!------------------------------------------------------------------------------!
4115     SUBROUTINE get_attribute_string( id, attribute_name, value, global,       &
4116                                      variable_name, no_abort )
4117
4118       USE pegrid
4119
4120       IMPLICIT NONE
4121
4122       CHARACTER(LEN=*)                ::  attribute_name   !< attribute name
4123       CHARACTER(LEN=*), OPTIONAL      ::  variable_name    !< variable name
4124       CHARACTER(LEN=*), INTENT(INOUT) ::  value            !< read value
4125
4126       INTEGER(iwp), INTENT(IN)    ::  id               !< file id
4127       INTEGER(iwp)                ::  id_var           !< variable id
4128
4129       LOGICAL ::  check_error                          !< flag indicating if handle_error shall be checked
4130       LOGICAL, INTENT(IN) ::  global                   !< flag indicating global attribute
4131       LOGICAL, INTENT(IN), OPTIONAL ::  no_abort       !< flag indicating if errors should be checked
4132#if defined( __netcdf )
4133
4134       IF ( PRESENT( no_abort ) )  THEN
4135          check_error = no_abort
4136       ELSE
4137          check_error = .TRUE.
4138       ENDIF
4139!
4140!--    Read global attribute
4141       IF ( global )  THEN
4142          nc_stat = NF90_GET_ATT( id, NF90_GLOBAL, TRIM( attribute_name ), value )
4143          IF ( check_error)  CALL handle_error( 'get_attribute_string global', 525, attribute_name )
4144!
4145!--    Read attributes referring to a single variable. Therefore, first inquire
4146!--    variable id
4147       ELSE
4148          nc_stat = NF90_INQ_VARID( id, TRIM( variable_name ), id_var )
4149          IF ( check_error)  CALL handle_error( 'get_attribute_string', 525, attribute_name )
4150
4151          nc_stat = NF90_GET_ATT( id, id_var, TRIM( attribute_name ), value )
4152          IF ( check_error)  CALL handle_error( 'get_attribute_string',525, attribute_name )
4153
4154       ENDIF
4155#endif
4156    END SUBROUTINE get_attribute_string
4157
4158
4159
4160!------------------------------------------------------------------------------!
4161! Description:
4162! ------------
4163!> Get dimension array for a given dimension
4164!------------------------------------------------------------------------------!
4165     SUBROUTINE get_dimension_length( id, dim_len, variable_name )
4166       USE pegrid
4167
4168       IMPLICIT NONE
4169
4170       CHARACTER(LEN=*)            ::  variable_name    !< dimension name
4171       CHARACTER(LEN=100)          ::  dum              !< dummy variable to receive return character
4172
4173       INTEGER(iwp)                ::  dim_len          !< dimension size
4174       INTEGER(iwp), INTENT(IN)    ::  id               !< file id
4175       INTEGER(iwp)                ::  id_dim           !< dimension id
4176
4177#if defined( __netcdf )
4178!
4179!--    First, inquire dimension ID
4180       nc_stat = NF90_INQ_DIMID( id, TRIM( variable_name ), id_dim )
4181       CALL handle_error( 'get_dimension_length', 526, variable_name )
4182!
4183!--    Inquire dimension length
4184       nc_stat = NF90_INQUIRE_DIMENSION( id, id_dim, dum, LEN = dim_len )
4185       CALL handle_error( 'get_dimension_length', 526, variable_name )
4186
4187#endif
4188    END SUBROUTINE get_dimension_length
4189
4190!------------------------------------------------------------------------------!
4191! Description:
4192! ------------
4193!> Routine for reading-in a character string from the chem emissions netcdf
4194!> input file. 
4195!------------------------------------------------------------------------------!
4196
4197    SUBROUTINE get_variable_string( id, variable_name, var_string, names_number )
4198
4199       USE indices
4200       USE pegrid
4201
4202       IMPLICIT NONE
4203
4204       CHARACTER (LEN=25), ALLOCATABLE, DIMENSION(:), INTENT(INOUT)  :: var_string
4205
4206       CHARACTER(LEN=*)                                              :: variable_name          !> variable name
4207
4208       CHARACTER (LEN=1), ALLOCATABLE, DIMENSION(:,:)                :: tmp_var_string         !> variable to be read
4209
4210
4211       INTEGER(iwp), INTENT(IN)                                      :: id                     !> file id
4212
4213       INTEGER(iwp), INTENT(IN)                                      :: names_number           !> number of names
4214
4215       INTEGER(iwp)                                                  :: id_var                 !> variable id
4216
4217       INTEGER(iwp)                                                  :: i,j                    !> index to go through the length of the dimensions
4218
4219       INTEGER(iwp)                                                  :: max_string_length=25   !> this is both the maximum length of a name, but also 
4220                                                                                            ! the number of the components of the first dimensions
4221                                                                                            ! (rows)
4222#if defined( __netcdf )
4223
4224       ALLOCATE(tmp_var_string(max_string_length,names_number))
4225
4226       ALLOCATE(var_string(names_number))
4227
4228    !-- Inquire variable id
4229       nc_stat = NF90_INQ_VARID( id, TRIM( variable_name ), id_var )
4230
4231
4232    !-- Get variable
4233    !-- Start cycle over the emission species
4234       DO i = 1, names_number
4235       !-- read the first letter of each component
4236          nc_stat = NF90_GET_VAR( id, id_var, var_string(i), start = (/ 1,i /), &
4237                                 count = (/ 1,1 /) )
4238          CALL handle_error( 'get_variable_string', 701 )
4239
4240       !-- Start cycle over charachters
4241          DO j = 1, max_string_length
4242                       
4243          !-- read the rest of the components of the name
4244             nc_stat = NF90_GET_VAR( id, id_var, tmp_var_string(j,i), start = (/ j,i /),&
4245                                     count = (/ 1,1 /) )
4246             CALL handle_error( 'get_variable_string', 702 )
4247
4248             IF ( iachar(tmp_var_string(j,i) ) == 0 ) THEN
4249                  tmp_var_string(j,i)=''
4250             ENDIF
4251
4252             IF ( j>1 ) THEN
4253             !-- Concatenate first letter of the name and the others
4254                var_string(i)=TRIM(var_string(i)) // TRIM(tmp_var_string(j,i))
4255
4256             ENDIF
4257          ENDDO
4258       ENDDO
4259
4260#endif
4261    END SUBROUTINE get_variable_string
4262
4263
4264!
4265!------------------------------------------------------------------------------!
4266! Description:
4267! ------------
4268!> Generalized routine for reading strings from a netCDF variable
4269!> to replace existing get_variable_string ( )
4270!>
4271!> Improvements:
4272!>   - Expanded string size limit from 25 to 512
4273!>   - No longer removes spaces between text magically (this seems to have
4274!>     been aimed at a very specific application, but I don't know what)
4275!>   - Simplified implementation
4276!>
4277!> Caveats:
4278!>   - Somehow I could not get the subroutine to work with str_array(:,:)
4279!>     so I reverted to a hard-coded str_array(:,512), hopefully large enough
4280!>     for most general applications.  This also means the character variable
4281!>     used for str_array must be of size (:,512)
4282!>     (ecc 20200128)   
4283!------------------------------------------------------------------------------!
4284
4285 SUBROUTINE get_variable_string_generic ( id, var_name, str_array, num_str, str_len )
4286
4287    IMPLICIT NONE
4288
4289    CHARACTER(LEN=*),                INTENT(IN)    :: var_name       !> netCDF variable name
4290    CHARACTER(LEN=512), ALLOCATABLE, INTENT(INOUT) :: str_array(:)   !> output string array
4291
4292    INTEGER(iwp)              :: buf_len   !> string buffer size
4293    INTEGER(iwp)              :: id_var    !> netCDF variable ID
4294    INTEGER(iwp)              :: k         !> generic counter
4295
4296    INTEGER(iwp), INTENT(IN)  :: id        !> netCDF file ID
4297    INTEGER(iwp), INTENT(IN)  :: num_str   !> number of string elements in array
4298    INTEGER(iwp), INTENT(IN)  :: str_len   !> size of each string element
4299
4300#if defined( __netcdf )
4301
4302!
4303!-- set buffer length to up to hard-coded string size
4304
4305    buf_len = MIN( ABS(str_len), 512 )
4306
4307!
4308!-- allocate necessary memories for string variables
4309
4310    ALLOCATE(str_array(num_str))
4311!
4312!-- get variable id
4313
4314    nc_stat = NF90_INQ_VARID( id, TRIM(var_name), id_var )
4315!
4316!-- extract string variables
4317
4318    DO k = 1, num_str
4319       str_array(k) = ''
4320       nc_stat = NF90_GET_VAR( id, id_var, str_array(k),  &
4321                      start = (/ 1, k /), count = (/ buf_len, 1 /)  )
4322       CALL handle_error ( 'get_variable_string_generic', 702 )
4323    ENDDO
4324
4325#endif
4326
4327 END SUBROUTINE get_variable_string_generic
4328
4329
4330!------------------------------------------------------------------------------!
4331! Description:
4332! ------------
4333!> Reads a character variable in a 1D array
4334!------------------------------------------------------------------------------!
4335     SUBROUTINE get_variable_1d_char( id, variable_name, var )
4336
4337       USE pegrid
4338
4339       IMPLICIT NONE
4340
4341       CHARACTER(LEN=*)            ::  variable_name          !< variable name
4342       CHARACTER(LEN=*), DIMENSION(:), INTENT(INOUT) ::  var  !< variable to be read
4343
4344       INTEGER(iwp)                ::  i                !< running index over variable dimension
4345       INTEGER(iwp), INTENT(IN)    ::  id               !< file id
4346       INTEGER(iwp)                ::  id_var           !< dimension id
4347       
4348       INTEGER(iwp), DIMENSION(2)  ::  dimid            !< dimension IDs
4349       INTEGER(iwp), DIMENSION(2)  ::  dimsize          !< dimension size
4350
4351#if defined( __netcdf )
4352
4353!
4354!--    First, inquire variable ID
4355       nc_stat = NF90_INQ_VARID( id, TRIM( variable_name ), id_var )
4356       CALL handle_error( 'get_variable_1d_int', 527, variable_name )
4357!
4358!--    Inquire dimension IDs
4359       nc_stat = NF90_INQUIRE_VARIABLE( id, id_var, dimids = dimid(1:2) )
4360       CALL handle_error( 'get_variable_1d_char', 527, variable_name )
4361!
4362!--    Read dimesnion length
4363       nc_stat = NF90_INQUIRE_DIMENSION( id, dimid(1), LEN = dimsize(1) )
4364       nc_stat = NF90_INQUIRE_DIMENSION( id, dimid(2), LEN = dimsize(2) )
4365       
4366!
4367!--    Read character array. Note, each element is read individually, in order
4368!--    to better separate single strings.
4369       DO  i = 1, dimsize(2)
4370          nc_stat = NF90_GET_VAR( id, id_var, var(i),                          &
4371                                  start = (/ 1, i /),                          &
4372                                  count = (/ dimsize(1), 1 /) )
4373          CALL handle_error( 'get_variable_1d_char', 527, variable_name )
4374       ENDDO     
4375                         
4376#endif
4377    END SUBROUTINE get_variable_1d_char
4378
4379   
4380!------------------------------------------------------------------------------!
4381! Description:
4382! ------------
4383!> Reads a 1D integer variable from file.
4384!------------------------------------------------------------------------------!
4385     SUBROUTINE get_variable_1d_int( id, variable_name, var )
4386
4387       USE pegrid
4388
4389       IMPLICIT NONE
4390
4391       CHARACTER(LEN=*)            ::  variable_name    !< variable name
4392
4393       INTEGER(iwp), INTENT(IN)    ::  id               !< file id
4394       INTEGER(iwp)                ::  id_var           !< dimension id
4395
4396       INTEGER(iwp), DIMENSION(:), INTENT(INOUT) ::  var  !< variable to be read
4397#if defined( __netcdf )
4398
4399!
4400!--    First, inquire variable ID
4401       nc_stat = NF90_INQ_VARID( id, TRIM( variable_name ), id_var )
4402       CALL handle_error( 'get_variable_1d_int', 527, variable_name )
4403!
4404!--    Inquire dimension length
4405       nc_stat = NF90_GET_VAR( id, id_var, var )
4406       CALL handle_error( 'get_variable_1d_int', 527, variable_name )
4407
4408#endif
4409    END SUBROUTINE get_variable_1d_int
4410
4411!------------------------------------------------------------------------------!
4412! Description:
4413! ------------
4414!> Reads a 1D float variable from file.
4415!------------------------------------------------------------------------------!
4416     SUBROUTINE get_variable_1d_real( id, variable_name, var )
4417
4418       USE pegrid
4419
4420       IMPLICIT NONE
4421
4422       CHARACTER(LEN=*)            ::  variable_name    !< variable name
4423
4424       INTEGER(iwp), INTENT(IN)    ::  id               !< file id
4425       INTEGER(iwp)                ::  id_var           !< dimension id
4426
4427       REAL(wp), DIMENSION(:), INTENT(INOUT) ::  var    !< variable to be read
4428#if defined( __netcdf )
4429
4430!
4431!--    First, inquire variable ID
4432       nc_stat = NF90_INQ_VARID( id, TRIM( variable_name ), id_var )
4433       CALL handle_error( 'get_variable_1d_real', 528, variable_name )
4434!
4435!--    Inquire dimension length
4436       nc_stat = NF90_GET_VAR( id, id_var, var )
4437       CALL handle_error( 'get_variable_1d_real', 528, variable_name )
4438
4439#endif
4440    END SUBROUTINE get_variable_1d_real
4441
4442
4443!------------------------------------------------------------------------------!
4444! Description:
4445! ------------
4446!> Reads a time-dependent 1D float variable from file.
4447!------------------------------------------------------------------------------!
4448    SUBROUTINE get_variable_pr( id, variable_name, t, var )
4449
4450       USE pegrid
4451
4452       IMPLICIT NONE
4453
4454       CHARACTER(LEN=*)                      ::  variable_name    !< variable name
4455
4456       INTEGER(iwp), INTENT(IN)              ::  id               !< file id
4457       INTEGER(iwp), DIMENSION(1:2)          ::  id_dim           !< dimension ids
4458       INTEGER(iwp)                          ::  id_var           !< dimension id
4459       INTEGER(iwp)                          ::  n_file           !< number of data-points in file along z dimension
4460       INTEGER(iwp), INTENT(IN)              ::  t                !< timestep number
4461
4462       REAL(wp), DIMENSION(:), INTENT(INOUT) ::  var  !< variable to be read
4463
4464#if defined( __netcdf )
4465!
4466!--    First, inquire variable ID
4467       nc_stat = NF90_INQ_VARID( id, TRIM( variable_name ), id_var )
4468!
4469!--    Inquire dimension size of vertical dimension
4470       nc_stat = NF90_INQUIRE_VARIABLE( id, id_var, DIMIDS = id_dim )
4471       nc_stat = NF90_INQUIRE_DIMENSION( id, id_dim(1), LEN = n_file )
4472!
4473!--    Read variable.
4474       nc_stat = NF90_GET_VAR( id, id_var, var,                                &
4475                               start = (/ 1,      t     /),                    &
4476                               count = (/ n_file, 1     /) )
4477       CALL handle_error( 'get_variable_pr', 529, variable_name )
4478
4479#endif
4480    END SUBROUTINE get_variable_pr
4481
4482
4483!------------------------------------------------------------------------------!
4484! Description:
4485! ------------
4486!> Reads a per-surface pars variable from file. Because all surfaces are stored
4487!> as flat 1-D array, each PE has to scan the data and find the surface indices
4488!> belonging to its subdomain. During this scan, it also builds a necessary
4489!> (j,i) index.
4490!------------------------------------------------------------------------------!
4491    SUBROUTINE get_variable_surf( id, variable_name, surf )
4492
4493       USE pegrid
4494
4495       USE indices,                                            &
4496           ONLY:  nxl, nxr, nys, nyn
4497
4498       USE control_parameters,                                 &
4499           ONLY: dz, message_string
4500
4501       USE grid_variables,                                     &
4502           ONLY: dx, dy
4503
4504       USE basic_constants_and_equations_mod,                  &
4505           ONLY: pi
4506
4507       IMPLICIT NONE
4508
4509       INTEGER(iwp), PARAMETER                   ::  nsurf_pars_read = 2**15 !< read buffer size (value > 10^15 makes problem with ifort)
4510
4511       CHARACTER(LEN=*)                          ::  variable_name !< variable name
4512
4513       INTEGER(iwp), DIMENSION(6)                ::  coords        !< integer coordinates of surface
4514       INTEGER(iwp)                              ::  i, j
4515       INTEGER(iwp)                              ::  isurf         !< netcdf surface index
4516       INTEGER(iwp)                              ::  is            !< local surface index
4517       INTEGER(iwp), INTENT(IN)                  ::  id            !< file id
4518       INTEGER(iwp), DIMENSION(2)                ::  id_dim        !< dimension ids
4519       INTEGER(iwp)                              ::  id_var        !< variable id
4520       INTEGER(iwp)                              ::  id_zs         !< zs variable id
4521       INTEGER(iwp)                              ::  id_ys         !< ys variable id
4522       INTEGER(iwp)                              ::  id_xs         !< xs variable id
4523       INTEGER(iwp)                              ::  id_zenith     !< zeith variable id
4524       INTEGER(iwp)                              ::  id_azimuth    !< azimuth variable id
4525       INTEGER(iwp)                              ::  is0, isc      !< read surface start and count
4526       INTEGER(iwp)                              ::  nsurf         !< total number of surfaces in file
4527       INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE ::  nsurf_ji      !< numbers of surfaces by coords
4528
4529       TYPE(pars_surf)                           ::  surf          !< parameters variable to be loaded
4530       REAL(wp), DIMENSION(:,:), ALLOCATABLE     ::  pars_read     !< read buffer
4531       REAL(wp), DIMENSION(:), ALLOCATABLE       ::  zs, ys, xs    !< read buffer for zs(s), ys, xs
4532       REAL(wp), DIMENSION(:), ALLOCATABLE       ::  zenith        !< read buffer for zenith(s)
4533       REAL(wp), DIMENSION(:), ALLOCATABLE       ::  azimuth       !< read buffer for azimuth(s)
4534       REAL(wp)                                  ::  oro_max_l     !< maximum terrain height under building
4535
4536#if defined( __netcdf )
4537!
4538!--    First, inquire variable ID
4539       nc_stat = NF90_INQ_VARID( id, TRIM( variable_name ), id_var )
4540       nc_stat = NF90_INQ_VARID( id, 'zs',                  id_zs )
4541       nc_stat = NF90_INQ_VARID( id, 'ys',                  id_ys )
4542       nc_stat = NF90_INQ_VARID( id, 'xs',                  id_xs )
4543       nc_stat = NF90_INQ_VARID( id, 'zenith',              id_zenith )
4544       nc_stat = NF90_INQ_VARID( id, 'azimuth',             id_azimuth )
4545!
4546!--    Inquire dimension sizes
4547       nc_stat = NF90_INQUIRE_VARIABLE( id, id_var, DIMIDS = id_dim )
4548       nc_stat = NF90_INQUIRE_DIMENSION( id, id_dim(1), LEN = nsurf )
4549       nc_stat = NF90_INQUIRE_DIMENSION( id, id_dim(2), LEN = surf%np )
4550
4551       ALLOCATE ( pars_read( nsurf_pars_read, surf%np ),        &
4552                  zs(nsurf_pars_read), ys(nsurf_pars_read),     &
4553                  xs(nsurf_pars_read), zenith(nsurf_pars_read), &
4554                  azimuth(nsurf_pars_read),                     &
4555                  nsurf_ji(nys:nyn, nxl:nxr) )
4556
4557       nsurf_ji(:,:) = 0
4558!
4559!--    Scan surface coordinates, count local
4560       is0 = 1
4561       DO
4562          isc = MIN(nsurf_pars_read, nsurf - is0 + 1)
4563          IF ( isc <= 0 )  EXIT
4564
4565          nc_stat = NF90_GET_VAR( id, id_ys, ys,     &
4566                                  start = (/ is0 /), &
4567                                  count = (/ isc /) )
4568          nc_stat = NF90_GET_VAR( id, id_xs, xs,     &
4569                                  start = (/ is0 /), &
4570                                  count = (/ isc /) )
4571          nc_stat = NF90_GET_VAR( id, id_zenith, zenith,      &
4572                                  start = (/ is0 /), &
4573                                  count = (/ isc /) )
4574          nc_stat = NF90_GET_VAR( id, id_azimuth, azimuth,    &
4575                                  start = (/ is0 /), &
4576                                  count = (/ isc /) )
4577          CALL handle_error( 'get_variable_surf', 682, 'azimuth' )
4578         
4579          DO  isurf = 1, isc
4580!
4581!--          Parse coordinates, detect if belongs to subdomain
4582             coords = transform_coords( xs(isurf), ys(isurf),         &
4583                                        zenith(isurf), azimuth(isurf) )
4584             IF ( coords(2) < nys  .OR.  coords(2) > nyn  .OR.  &
4585                  coords(3) < nxl  .OR.  coords(3) > nxr )  CYCLE
4586
4587             nsurf_ji(coords(2), coords(3)) = nsurf_ji(coords(2), coords(3)) + 1
4588          ENDDO
4589
4590          is0 = is0 + isc
4591       ENDDO
4592!
4593!--    Populate reverse index from surface counts
4594       ALLOCATE ( surf%index_ji( 2, nys:nyn, nxl:nxr ) )
4595
4596       isurf = 1
4597       DO  j = nys, nyn
4598          DO  i = nxl, nxr
4599             surf%index_ji(:,j,i) = (/ isurf, isurf + nsurf_ji(j,i) - 1 /)
4600             isurf = isurf + nsurf_ji(j,i)
4601          ENDDO
4602       ENDDO
4603
4604       surf%nsurf = isurf - 1
4605       ALLOCATE( surf%pars( 0:surf%np-1, surf%nsurf ), &
4606                 surf%coords( 6, surf%nsurf ) )
4607!
4608!--    Scan surfaces again, saving pars into allocated structures
4609       nsurf_ji(:,:) = 0
4610       is0 = 1
4611       DO
4612          isc = MIN(nsurf_pars_read, nsurf - is0 + 1)
4613          IF ( isc <= 0 )  EXIT
4614
4615          nc_stat = NF90_GET_VAR( id, id_var, pars_read(1:isc, 1:surf%np), &
4616                                  start = (/ is0, 1       /),              &
4617                                  count = (/ isc, surf%np /) )
4618          CALL handle_error( 'get_variable_surf', 683, variable_name )
4619
4620          nc_stat = NF90_GET_VAR( id, id_zs, zs,                           &
4621                                  start = (/ is0 /),                       &
4622                                  count = (/ isc /) )
4623          nc_stat = NF90_GET_VAR( id, id_ys, ys,                           &
4624                                  start = (/ is0 /),                       &
4625                                  count = (/ isc /) )
4626          nc_stat = NF90_GET_VAR( id, id_xs, xs,                           &
4627                                  start = (/ is0 /),                       &
4628                                  count = (/ isc /) )
4629          nc_stat = NF90_GET_VAR( id, id_zenith, zenith,                   &
4630                                  start = (/ is0 /),                       &
4631                                  count = (/ isc /) )
4632          nc_stat = NF90_GET_VAR( id, id_azimuth, azimuth,                 &
4633                                  start = (/ is0 /),                       &
4634                                  count = (/ isc /) )
4635         
4636          DO  isurf = 1, isc
4637!
4638!--          Parse coordinates, detect if belongs to subdomain
4639             coords = transform_coords( xs(isurf), ys(isurf),         &
4640                                        zenith(isurf), azimuth(isurf) )
4641             IF ( coords(2) < nys  .OR.  coords(2) > nyn  .OR.  &
4642                  coords(3) < nxl  .OR.  coords(3) > nxr )  CYCLE
4643!
4644!--          Determine maximum terrain under building (base z-coordinate). Using
4645!--          normal vector to locate building inner coordinates.
4646             oro_max_l = buildings_f%oro_max(coords(2)-coords(5), coords(3)-coords(6))
4647             IF  ( oro_max_l == buildings_f%fill1 )  THEN
4648                WRITE( message_string, * ) 'Found building surface on '   // &
4649                   'non-building coordinates (xs, ys, zenith, azimuth): ',   &
4650                   xs(isurf), ys(isurf), zenith(isurf), azimuth(isurf)
4651                CALL message( 'get_variable_surf', 'PA0684', 2, 2, myid, 6, 0 ) 
4652             ENDIF
4653!
4654!--          Urban layer has no stretching, therefore using dz(1) instead of linear
4655!--          searching through zu/zw
4656             coords(1) = NINT((zs(isurf) + oro_max_l) / dz(1) +     &
4657                              0.5_wp + 0.5_wp * coords(4), KIND=iwp)
4658!
4659!--          Save surface entry
4660             is = surf%index_ji(1, coords(2), coords(3)) + nsurf_ji(coords(2), coords(3))
4661             surf%pars(:,is) = pars_read(isurf,:)
4662             surf%coords(:,is) = coords(:)
4663
4664             nsurf_ji(coords(2), coords(3)) = nsurf_ji(coords(2), coords(3)) + 1
4665          ENDDO
4666
4667          is0 = is0 + isc
4668       ENDDO
4669
4670       DEALLOCATE( pars_read, zs, ys, xs, zenith, azimuth, nsurf_ji )
4671
4672    CONTAINS
4673
4674       PURE FUNCTION transform_coords( x, y, zenith, azimuth )
4675
4676          REAL(wp), INTENT(in)       ::  x, y    !< surface centre coordinates in metres from origin
4677          REAL(wp), INTENT(in)       ::  zenith  !< surface normal zenith angle in degrees
4678          REAL(wp), INTENT(in)       ::  azimuth !< surface normal azimuth angle in degrees
4679
4680          INTEGER(iwp), DIMENSION(6) ::  transform_coords !< (k,j,i,norm_z,norm_y,norm_x)
4681
4682          transform_coords(4) = NINT(COS(zenith*pi/180._wp), KIND=iwp)
4683          IF ( transform_coords(4) == 0 )  THEN
4684             transform_coords(5) = NINT(COS(azimuth*pi/180._wp), KIND=iwp)
4685             transform_coords(6) = NINT(SIN(azimuth*pi/180._wp), KIND=iwp)
4686          ELSE
4687             transform_coords(5) = 0._wp
4688             transform_coords(6) = 0._wp
4689          ENDIF
4690
4691          transform_coords(1) = -999._wp ! not calculated here
4692          transform_coords(2) = NINT(y/dy - 0.5_wp + 0.5_wp*transform_coords(5), KIND=iwp)
4693          transform_coords(3) = NINT(x/dx - 0.5_wp + 0.5_wp*transform_coords(6), KIND=iwp)
4694
4695       END FUNCTION transform_coords
4696
4697#endif
4698    END SUBROUTINE get_variable_surf
4699
4700
4701!------------------------------------------------------------------------------!
4702! Description:
4703! ------------
4704!> Reads a 2D REAL variable from a file. Reading is done processor-wise,
4705!> i.e. each core reads its own domain in slices along x.
4706!------------------------------------------------------------------------------!
4707    SUBROUTINE get_variable_2d_real( id, variable_name, var, is, ie, js, je )
4708
4709       USE indices
4710       USE pegrid
4711
4712       IMPLICIT NONE
4713
4714       CHARACTER(LEN=*)              ::  variable_name   !< variable name
4715
4716       INTEGER(iwp)                  ::  i               !< running index along x direction
4717       INTEGER(iwp)                  ::  ie              !< start index for subdomain input along x direction
4718       INTEGER(iwp)                  ::  is              !< end index for subdomain input along x direction
4719       INTEGER(iwp), INTENT(IN)      ::  id              !< file id
4720       INTEGER(iwp)                  ::  id_var          !< variable id
4721       INTEGER(iwp)                  ::  j               !< running index along y direction
4722       INTEGER(iwp)                  ::  je              !< start index for subdomain input along y direction
4723       INTEGER(iwp)                  ::  js              !< end index for subdomain input along y direction
4724       
4725       REAL(wp), DIMENSION(:,:), ALLOCATABLE   ::  tmp   !< temporary variable to read data from file according
4726                                                         !< to its reverse memory access
4727       REAL(wp), DIMENSION(:,:), INTENT(INOUT) ::  var   !< variable to be read
4728#if defined( __netcdf )
4729!
4730!--    Inquire variable id
4731       nc_stat = NF90_INQ_VARID( id, TRIM( variable_name ), id_var )
4732!
4733!--    Check for collective read-operation and set respective NetCDF flags if
4734!--    required.
4735       IF ( collective_read )  THEN
4736#if defined( __netcdf4_parallel )
4737          nc_stat = NF90_VAR_PAR_ACCESS (id, id_var, NF90_COLLECTIVE)
4738#endif
4739       ENDIF
4740
4741!
4742!-- Allocate temporary variable according to memory access on file.
4743       ALLOCATE( tmp(is:ie,js:je) )
4744!
4745!-- Get variable
4746       nc_stat = NF90_GET_VAR( id, id_var, tmp,            &
4747                      start = (/ is+1,      js+1 /),       &
4748                      count = (/ ie-is + 1, je-js+1 /) )   
4749          CALL handle_error( 'get_variable_2d_real', 530, variable_name )
4750!
4751!-- Resort data. Please note, dimension subscripts of var all start at 1.
4752          DO  i = is, ie
4753             DO  j = js, je
4754                var(j-js+1,i-is+1) = tmp(i,j)
4755             ENDDO
4756          ENDDO
4757       
4758          DEALLOCATE( tmp )
4759
4760#endif
4761    END SUBROUTINE get_variable_2d_real
4762
4763!------------------------------------------------------------------------------!
4764! Description:
4765! ------------
4766!> Reads a 2D 32-bit INTEGER variable from file. Reading is done processor-wise,
4767!> i.e. each core reads its own domain in slices along x.
4768!------------------------------------------------------------------------------!
4769    SUBROUTINE get_variable_2d_int32( id, variable_name, var, is, ie, js, je )
4770
4771       USE indices
4772       USE pegrid
4773
4774       IMPLICIT NONE
4775
4776       CHARACTER(LEN=*)              ::  variable_name   !< variable name
4777
4778       INTEGER(iwp)                  ::  i               !< running index along x direction
4779       INTEGER(iwp)                  ::  ie              !< start index for subdomain input along x direction
4780       INTEGER(iwp)                  ::  is              !< end index for subdomain input along x direction
4781       INTEGER(iwp), INTENT(IN)      ::  id              !< file id
4782       INTEGER(iwp)                  ::  id_var          !< variable id
4783       INTEGER(iwp)                  ::  j               !< running index along y direction
4784       INTEGER(iwp)                  ::  je              !< start index for subdomain input along y direction
4785       INTEGER(iwp)                  ::  js              !< end index for subdomain input along y direction
4786       
4787       INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE   ::  tmp  !< temporary variable to read data from file according
4788                                                            !< to its reverse memory access
4789       INTEGER(iwp), DIMENSION(:,:), INTENT(INOUT) ::  var  !< variable to be read
4790#if defined( __netcdf )
4791!
4792!--    Inquire variable id
4793       nc_stat = NF90_INQ_VARID( id, TRIM( variable_name ), id_var )
4794!
4795!--    Check for collective read-operation and set respective NetCDF flags if
4796!--    required.
4797       IF ( collective_read )  THEN
4798#if defined( __netcdf4_parallel )       
4799          nc_stat = NF90_VAR_PAR_ACCESS (id, id_var, NF90_COLLECTIVE)
4800#endif
4801       ENDIF
4802!
4803!--    Allocate temporary variable according to memory access on file.
4804       ALLOCATE( tmp(is:ie,js:je) )
4805!
4806!--    Get variable
4807       nc_stat = NF90_GET_VAR( id, id_var, tmp,                                &
4808                               start = (/ is+1,      js+1 /),                  &
4809                               count = (/ ie-is + 1, je-js+1 /) )   
4810                               
4811       CALL handle_error( 'get_variable_2d_int32', 531, variable_name )                             
4812!
4813!--    Resort data. Please note, dimension subscripts of var all start at 1.
4814       DO  i = is, ie
4815          DO  j = js, je
4816             var(j-js+1,i-is+1) = tmp(i,j)
4817          ENDDO
4818       ENDDO
4819       
4820       DEALLOCATE( tmp )
4821
4822#endif
4823    END SUBROUTINE get_variable_2d_int32
4824
4825!------------------------------------------------------------------------------!
4826! Description:
4827! ------------
4828!> Reads a 2D 8-bit INTEGER variable from file. Reading is done processor-wise,
4829!> i.e. each core reads its own domain in slices along x.
4830!------------------------------------------------------------------------------!
4831    SUBROUTINE get_variable_2d_int8( id, variable_name, var, is, ie, js, je )
4832
4833       USE indices
4834       USE pegrid
4835
4836       IMPLICIT NONE
4837
4838       CHARACTER(LEN=*)              ::  variable_name   !< variable name
4839
4840       INTEGER(iwp)                  ::  i               !< running index along x direction
4841       INTEGER(iwp)                  ::  ie              !< start index for subdomain input along x direction
4842       INTEGER(iwp)                  ::  is              !< end index for subdomain input along x direction
4843       INTEGER(iwp), INTENT(IN)      ::  id              !< file id
4844       INTEGER(iwp)                  ::  id_var          !< variable id
4845       INTEGER(iwp)                  ::  j               !< running index along y direction
4846       INTEGER(iwp)                  ::  je              !< start index for subdomain input along y direction
4847       INTEGER(iwp)                  ::  js              !< end index for subdomain input along y direction
4848       
4849       INTEGER(KIND=1), DIMENSION(:,:), ALLOCATABLE   ::  tmp  !< temporary variable to read data from file according
4850                                                               !< to its reverse memory access
4851       INTEGER(KIND=1), DIMENSION(:,:), INTENT(INOUT) ::  var  !< variable to be read
4852#if defined( __netcdf )
4853!
4854!--    Inquire variable id
4855       nc_stat = NF90_INQ_VARID( id, TRIM( variable_name ), id_var )
4856!
4857!--    Check for collective read-operation and set respective NetCDF flags if
4858!--    required.
4859       IF ( collective_read )  THEN
4860#if defined( __netcdf4_parallel )       
4861          nc_stat = NF90_VAR_PAR_ACCESS (id, id_var, NF90_COLLECTIVE)
4862#endif         
4863       ENDIF
4864!
4865!--    Allocate temporary variable according to memory access on file.
4866       ALLOCATE( tmp(is:ie,js:je) )
4867!
4868!--    Get variable
4869       nc_stat = NF90_GET_VAR( id, id_var, tmp,                                &
4870                               start = (/ is+1,      js+1 /),                  &
4871                               count = (/ ie-is + 1, je-js+1 /) )   
4872                               
4873       CALL handle_error( 'get_variable_2d_int8', 532, variable_name )
4874!
4875!--    Resort data. Please note, dimension subscripts of var all start at 1.
4876       DO  i = is, ie
4877          DO  j = js, je
4878             var(j-js+1,i-is+1) = tmp(i,j)
4879          ENDDO
4880       ENDDO
4881       
4882       DEALLOCATE( tmp )
4883
4884#endif
4885    END SUBROUTINE get_variable_2d_int8
4886
4887
4888!------------------------------------------------------------------------------!
4889! Description:
4890! ------------
4891!> Reads a 3D 8-bit INTEGER variable from file.
4892!------------------------------------------------------------------------------!
4893    SUBROUTINE get_variable_3d_int8( id, variable_name, var, is, ie, js, je,   &
4894                                     ks, ke )
4895
4896       USE indices
4897       USE pegrid
4898
4899       IMPLICIT NONE
4900
4901       CHARACTER(LEN=*)              ::  variable_name   !< variable name
4902
4903       INTEGER(iwp)                  ::  i               !< index along x direction
4904       INTEGER(iwp)                  ::  ie              !< start index for subdomain input along x direction
4905       INTEGER(iwp)                  ::  is              !< end index for subdomain input along x direction
4906       INTEGER(iwp), INTENT(IN)      ::  id              !< file id
4907       INTEGER(iwp)                  ::  id_var          !< variable id
4908       INTEGER(iwp)                  ::  j               !< index along y direction
4909       INTEGER(iwp)                  ::  je              !< start index for subdomain input along y direction
4910       INTEGER(iwp)                  ::  js              !< end index for subdomain input along y direction
4911       INTEGER(iwp)                  ::  k               !< index along any 3rd dimension
4912       INTEGER(iwp)                  ::  ke              !< start index of 3rd dimension
4913       INTEGER(iwp)                  ::  ks              !< end index of 3rd dimension
4914
4915       INTEGER(KIND=1), DIMENSION(:,:,:), ALLOCATABLE   ::  tmp  !< temporary variable to read data from file according
4916                                                                 !< to its reverse memory access
4917
4918       INTEGER(KIND=1), DIMENSION(:,:,:), INTENT(INOUT) ::  var  !< variable to be read
4919#if defined( __netcdf )
4920
4921!
4922!--    Inquire variable id
4923       nc_stat = NF90_INQ_VARID( id, TRIM( variable_name ), id_var )
4924!
4925!--    Check for collective read-operation and set respective NetCDF flags if
4926!--    required.
4927       IF ( collective_read )  THEN
4928#if defined( __netcdf4_parallel )
4929          nc_stat = NF90_VAR_PAR_ACCESS (id, id_var, NF90_COLLECTIVE)
4930#endif
4931       ENDIF
4932!
4933!--    Allocate temporary variable according to memory access on file.
4934       ALLOCATE( tmp(is:ie,js:je,ks:ke) )
4935!
4936!--    Get variable
4937       nc_stat = NF90_GET_VAR( id, id_var, tmp,                                &
4938                               start = (/ is+1,    js+1,    ks+1 /),           &
4939                               count = (/ ie-is+1, je-js+1, ke-ks+1 /) )
4940
4941       CALL handle_error( 'get_variable_3d_int8', 533, variable_name )
4942!
4943!--    Resort data. Please note, dimension subscripts of var all start at 1.
4944       DO  i = is, ie
4945          DO  j = js, je
4946             DO  k = ks, ke
4947                var(k-ks+1,j-js+1,i-is+1) = tmp(i,j,k)
4948             ENDDO
4949          ENDDO
4950       ENDDO
4951
4952       DEALLOCATE( tmp )
4953
4954#endif
4955    END SUBROUTINE get_variable_3d_int8
4956
4957
4958!------------------------------------------------------------------------------!
4959! Description:
4960! ------------
4961!> Reads a 3D float variable from file.
4962!------------------------------------------------------------------------------!
4963    SUBROUTINE get_variable_3d_real( id, variable_name, var, is, ie, js, je,   &
4964                                     ks, ke )
4965
4966       USE indices
4967       USE pegrid
4968
4969       IMPLICIT NONE
4970
4971       CHARACTER(LEN=*)              ::  variable_name   !< variable name
4972
4973       INTEGER(iwp)                  ::  i               !< index along x direction
4974       INTEGER(iwp)                  ::  ie              !< start index for subdomain input along x direction
4975       INTEGER(iwp)                  ::  is              !< end index for subdomain input along x direction
4976       INTEGER(iwp), INTENT(IN)      ::  id              !< file id
4977       INTEGER(iwp)                  ::  id_var          !< variable id
4978       INTEGER(iwp)                  ::  j               !< index along y direction
4979       INTEGER(iwp)                  ::  je              !< start index for subdomain input along y direction
4980       INTEGER(iwp)                  ::  js              !< end index for subdomain input along y direction
4981       INTEGER(iwp)                  ::  k               !< index along any 3rd dimension
4982       INTEGER(iwp)                  ::  ke              !< start index of 3rd dimension
4983       INTEGER(iwp)                  ::  ks              !< end index of 3rd dimension
4984
4985       REAL(wp), DIMENSION(:,:,:), ALLOCATABLE   ::  tmp !< temporary variable to read data from file according
4986                                                         !< to its reverse memory access
4987
4988       REAL(wp), DIMENSION(:,:,:), INTENT(INOUT) ::  var !< variable to be read
4989#if defined( __netcdf )
4990
4991!
4992!--    Inquire variable id
4993       nc_stat = NF90_INQ_VARID( id, TRIM( variable_name ), id_var ) 
4994!
4995!--    Check for collective read-operation and set respective NetCDF flags if
4996!--    required.
4997       IF ( collective_read )  THEN
4998#if defined( __netcdf4_parallel )
4999          nc_stat = NF90_VAR_PAR_ACCESS (id, id_var, NF90_COLLECTIVE)
5000#endif
5001       ENDIF
5002!
5003!--    Allocate temporary variable according to memory access on file.
5004       ALLOCATE( tmp(is:ie,js:je,ks:ke) )
5005!
5006!--    Get variable
5007       nc_stat = NF90_GET_VAR( id, id_var, tmp,                                &
5008                               start = (/ is+1,    js+1,    ks+1 /),           &
5009                               count = (/ ie-is+1, je-js+1, ke-ks+1 /) )   
5010
5011       CALL handle_error( 'get_variable_3d_real', 534, variable_name )
5012!
5013!--    Resort data. Please note, dimension subscripts of var all start at 1.
5014       DO  i = is, ie
5015          DO  j = js, je
5016             DO  k = ks, ke
5017                var(k-ks+1,j-js+1,i-is+1) = tmp(i,j,k)
5018             ENDDO
5019          ENDDO
5020       ENDDO
5021
5022       DEALLOCATE( tmp )
5023
5024#endif
5025    END SUBROUTINE get_variable_3d_real
5026
5027!------------------------------------------------------------------------------!
5028! Description:
5029! ------------
5030!> Reads a 4D float variable from file.
5031!------------------------------------------------------------------------------!
5032    SUBROUTINE get_variable_4d_real( id, variable_name, var, is, ie, js, je,   &
5033                                     k1s, k1e, k2s, k2e )
5034
5035       USE indices
5036       USE pegrid
5037
5038       IMPLICIT NONE
5039
5040       CHARACTER(LEN=*)              ::  variable_name   !< variable name
5041
5042       INTEGER(iwp)                  ::  i               !< index along x direction
5043       INTEGER(iwp)                  ::  ie              !< start index for subdomain input along x direction
5044       INTEGER(iwp)                  ::  is              !< end index for subdomain input along x direction
5045       INTEGER(iwp), INTENT(IN)      ::  id              !< file id
5046       INTEGER(iwp)                  ::  id_var          !< variable id
5047       INTEGER(iwp)                  ::  j               !< index along y direction
5048       INTEGER(iwp)                  ::  je              !< start index for subdomain input along y direction
5049       INTEGER(iwp)                  ::  js              !< end index for subdomain input along y direction
5050       INTEGER(iwp)                  ::  k1              !< index along 3rd direction
5051       INTEGER(iwp)                  ::  k1e             !< start index for 3rd dimension
5052       INTEGER(iwp)                  ::  k1s             !< end index for 3rd dimension
5053       INTEGER(iwp)                  ::  k2              !< index along 4th direction
5054       INTEGER(iwp)                  ::  k2e             !< start index for 4th dimension
5055       INTEGER(iwp)                  ::  k2s             !< end index for 4th dimension
5056
5057       REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE   ::  tmp  !< temporary variable to read data from file according
5058                                                            !< to its reverse memory access
5059       REAL(wp), DIMENSION(:,:,:,:), INTENT(INOUT) ::  var  !< variable to be read
5060#if defined( __netcdf )
5061
5062!
5063!--    Inquire variable id
5064       nc_stat = NF90_INQ_VARID( id, TRIM( variable_name ), id_var )
5065!
5066!--    Check for collective read-operation and set respective NetCDF flags if
5067!--    required.
5068       IF ( collective_read )  THEN
5069#if defined( __netcdf4_parallel )       
5070          nc_stat = NF90_VAR_PAR_ACCESS (id, id_var, NF90_COLLECTIVE)
5071#endif
5072       ENDIF
5073
5074!
5075!-- Allocate temporary variable according to memory access on file.
5076       ALLOCATE( tmp(is:ie,js:je,k1s:k1e,k2s:k2e) )
5077!
5078!-- Get variable
5079       nc_stat = NF90_GET_VAR( id, id_var, tmp,                                &
5080                      start = (/ is+1,    js+1,    k1s+1,     k2s+1 /),        &
5081                      count = (/ ie-is+1, je-js+1, k1e-k1s+1, k2e-k2s+1 /) )
5082
5083          CALL handle_error( 'get_variable_4d_real', 535, variable_name )
5084!
5085!-- Resort data. Please note, dimension subscripts of var all start at 1.
5086       DO  i = is, ie
5087          DO  j = js, je
5088             DO  k1 = k1s, k1e
5089                DO  k2 = k2s, k2e
5090                   var(k2-k2s+1,k1-k1s+1,j-js+1,i-is+1) = tmp(i,j,k1,k2)
5091                ENDDO
5092             ENDDO
5093          ENDDO
5094       ENDDO
5095
5096       DEALLOCATE( tmp )
5097
5098#endif
5099
5100    END SUBROUTINE get_variable_4d_real
5101
5102!------------------------------------------------------------------------------!
5103! Description:
5104! ------------
5105!> Reads a 4D float variable from file and store it to a 3-d variable.
5106!------------------------------------------------------------------------------!
5107    SUBROUTINE get_variable_4d_to_3d_real( id, variable_name, var, ns, is, ie, js, je,   &
5108                                           ks, ke )
5109
5110       USE indices
5111       USE pegrid
5112
5113       IMPLICIT NONE
5114
5115       CHARACTER(LEN=*)              ::  variable_name   !< variable name
5116
5117       INTEGER(iwp)                  ::  i               !< index along x direction
5118       INTEGER(iwp)                  ::  ie              !< end index for subdomain input along x direction
5119       INTEGER(iwp)                  ::  is              !< start index for subdomain input along x direction
5120       INTEGER(iwp), INTENT(IN)      ::  id              !< file id
5121       INTEGER(iwp)                  ::  id_var          !< variable id
5122       INTEGER(iwp)                  ::  j               !< index along y direction
5123       INTEGER(iwp)                  ::  je              !< end index for subdomain input along y direction
5124       INTEGER(iwp)                  ::  js              !< start index for subdomain input along y direction
5125       INTEGER(iwp)                  ::  k               !< index along any 4th dimension
5126       INTEGER(iwp)                  ::  ke              !< end index of 4th dimension
5127       INTEGER(iwp)                  ::  ks              !< start index of 4th dimension
5128       INTEGER(iwp)                  ::  ns              !< start index for subdomain input along n dimension
5129
5130       REAL</