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

Last change on this file since 4400 was 4400, checked in by suehring, 16 months ago

Revision of the virtual-measurement module: data input from NetCDF file; removed binary output - instead parallel NetCDF output using the new data-output module; variable attributes added; further variables added that can be sampled, file connections added; Functions for coordinate transformation moved to basic_constants_and_equations; netcdf_data_input_mod: unused routines netcdf_data_input_att and netcdf_data_input_var removed; new routines to inquire fill values added; Preprocessing script (palm_cvd) to setup virtual-measurement input files provided; postprocessor combine_virtual_measurements deleted

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