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

Last change on this file since 4435 was 4435, checked in by raasch, 3 years ago

bugfix for message that reports about files that are read from in case that the virtual PE grid has chenged (in case of large number of files format was exceeded), detailed messages about the files are now output to the debug file, temporary bugfix to avoid compile problems with older NetCDFD libraries on IMUK machines

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