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

Last change on this file since 4404 was 4404, checked in by suehring, 20 months ago

Fix misplaced preprocessor directives in netcdf data input

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