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

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

Define a default list of coordinate reference system variables used when no static driver input is available

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