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

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

ghost point exchange modularized, bugfix for wrong 2d-exchange

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