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

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

Mesoscale offline nesting: enable LOD 1 (homogeneous) input of lateral and top boundary conditions; add new generic subroutines to read time-dependent profile data from dynamic input file; minor bugfix - add missing initialization of the top boundary

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