Home

Context Navigation

netcdf_data_input_mod.f90 @ 4346

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