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

Last change on this file since 4641 was 4641, checked in by suehring, 5 years ago

Virtual measurement output: change content of attribute data_content, minor changes in long-name attributes and global attributes to be in agreement with (UC)2 data standard; driver creation config files: remove dummy input strings and change default values (to be in agreement with (UC)2 data standard); palm_csd_netcdf_interface: correct wrong false_easting attribute in crs variable

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