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

Last change on this file since 2897 was 2897, checked in by suehring, 4 years ago

Relax restrictions for topography input via static input file, terrain and building heights, as well as building IDs can be input separately and are not mandatory any more.

  • Property svn:keywords set to Id
File size: 189.1 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-2018 Leibniz Universitaet Hannover
18!------------------------------------------------------------------------------!
19!
20! Current revisions:
21! -----------------
22!
23!
24! Former revisions:
25! -----------------
26! $Id: netcdf_data_input_mod.f90 2897 2018-03-15 11:47:16Z suehring $
27! Relax restrictions for topography input, terrain and building heights can be
28! input separately and are not mandatory any more.
29!
30! 2874 2018-03-13 10:55:42Z knoop
31! Bugfix: wrong placement of netcdf cpp-macros fixed
32!
33! 2794 2018-02-07 14:09:43Z knoop
34! Check if 3D building input is consistent to numeric grid.
35!
36! 2773 2018-01-30 14:12:54Z suehring
37! - Enable initialization with 3D topography.
38! - Move check for correct initialization in nesting mode to check_parameters.
39!
40! 2772 2018-01-29 13:10:35Z suehring
41! Initialization of simulation independent on land-surface model.
42!
43! 2746 2018-01-15 12:06:04Z suehring
44! Read plant-canopy variables independently on land-surface model usage
45!
46! 2718 2018-01-02 08:49:38Z maronga
47! Corrected "Former revisions" section
48!
49! 2711 2017-12-20 17:04:49Z suehring
50! Rename subroutine close_file to avoid double-naming. 
51!
52! 2700 2017-12-15 14:12:35Z suehring
53!
54! 2696 2017-12-14 17:12:51Z kanani
55! Initial revision (suehring)
56!
57!
58!
59!
60! Authors:
61! --------
62! @author Matthias Suehring
63!
64! Description:
65! ------------
66!> Modulue contains routines to input data according to Palm input data       
67!> standart using dynamic and static input files.
68!>
69!> @todo - Order input alphabetically
70!> @todo - Revise error messages and error numbers
71!> @todo - Input of missing quantities (chemical species, emission rates)
72!> @todo - Defninition and input of still missing variable attributes
73!> @todo - Input of initial geostrophic wind profiles with cyclic conditions.
74!------------------------------------------------------------------------------!
75 MODULE netcdf_data_input_mod
76 
77    USE control_parameters,                                                    &
78        ONLY:  coupling_char, io_blocks, io_group
79
80    USE kinds
81
82#if defined ( __netcdf )
83    USE NETCDF
84#endif
85
86    USE pegrid
87
88!
89!-- Define data type for nesting in larger-scale models like COSMO.
90!-- Data type comprises u, v, w, pt, and q at lateral and top boundaries.
91    TYPE force_type
92
93       CHARACTER(LEN=100), DIMENSION(:), ALLOCATABLE ::  var_names
94
95       INTEGER(iwp) ::  nt     !< number of time levels in dynamic input file
96       INTEGER(iwp) ::  nzu    !< number of vertical levels on scalar grid in dynamic input file
97       INTEGER(iwp) ::  nzw    !< number of vertical levels on w grid in dynamic input file
98       INTEGER(iwp) ::  tind   !< time index for reference time in large-scale forcing data
99       INTEGER(iwp) ::  tind_p !< time index for following time in large-scale forcing data
100
101       LOGICAL      ::  init         = .FALSE.
102       LOGICAL      ::  interpolated = .FALSE.
103       LOGICAL      ::  from_file    = .FALSE.
104
105       REAL(wp), DIMENSION(:), ALLOCATABLE ::  surface_pressure !< time dependent surface pressure
106       REAL(wp), DIMENSION(:), ALLOCATABLE ::  time             !< time levels in dynamic input file
107       REAL(wp), DIMENSION(:), ALLOCATABLE ::  zu_atmos         !< vertical levels at scalar grid in dynamic input file
108       REAL(wp), DIMENSION(:), ALLOCATABLE ::  zw_atmos         !< vertical levels at w grid in dynamic input file
109
110       REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  u_left   !< u-component at left boundary
111       REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  v_left   !< v-component at left boundary
112       REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  w_left   !< w-component at left boundary
113       REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  q_left   !< mixing ratio at left boundary
114       REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  pt_left  !< potentital temperautre at left boundary
115
116       REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  u_north  !< u-component at north boundary
117       REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  v_north  !< v-component at north boundary
118       REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  w_north  !< w-component at north boundary
119       REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  q_north  !< mixing ratio at north boundary
120       REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  pt_north !< potentital temperautre at north boundary
121
122       REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  u_right  !< u-component at right boundary
123       REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  v_right  !< v-component at right boundary
124       REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  w_right  !< w-component at right boundary
125       REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  q_right  !< mixing ratio at right boundary
126       REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  pt_right !< potentital temperautre at right boundary
127
128       REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  u_south  !< u-component at south boundary
129       REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  v_south  !< v-component at south boundary
130       REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  w_south  !< w-component at south boundary
131       REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  q_south  !< mixing ratio at south boundary
132       REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  pt_south !< potentital temperautre at south boundary
133
134       REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  u_top    !< u-component at top boundary
135       REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  v_top    !< v-component at top boundary
136       REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  w_top    !< w-component at top boundary
137       REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  q_top    !< mixing ratio at top boundary
138       REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  pt_top   !< potentital temperautre at top boundary
139
140    END TYPE force_type
141
142    TYPE init_type
143
144       INTEGER(iwp) ::  lod_msoil !< level of detail - soil moisture
145       INTEGER(iwp) ::  lod_pt    !< level of detail - pt
146       INTEGER(iwp) ::  lod_q     !< level of detail - q
147       INTEGER(iwp) ::  lod_tsoil !< level of detail - soil temperature
148       INTEGER(iwp) ::  lod_u     !< level of detail - u-component
149       INTEGER(iwp) ::  lod_v     !< level of detail - v-component
150       INTEGER(iwp) ::  lod_w     !< level of detail - w-component
151       INTEGER(iwp) ::  nx        !< number of scalar grid points along x in dynamic input file
152       INTEGER(iwp) ::  nxu       !< number of u grid points along x in dynamic input file
153       INTEGER(iwp) ::  ny        !< number of scalar grid points along y in dynamic input file
154       INTEGER(iwp) ::  nyv       !< number of v grid points along y in dynamic input file
155       INTEGER(iwp) ::  nzs       !< number of vertical soil levels in dynamic input file
156       INTEGER(iwp) ::  nzu       !< number of vertical levels on scalar grid in dynamic input file
157       INTEGER(iwp) ::  nzw       !< number of vertical levels on w grid in dynamic input file
158
159       LOGICAL ::  from_file_msoil  = .FALSE. !< flag indicating whether soil moisture is already initialized from file
160       LOGICAL ::  from_file_pt     = .FALSE. !< flag indicating whether pt is already initialized from file
161       LOGICAL ::  from_file_q      = .FALSE. !< flag indicating whether q is already initialized from file 
162       LOGICAL ::  from_file_tsoil  = .FALSE. !< flag indicating whether soil temperature is already initialized from file
163       LOGICAL ::  from_file_u      = .FALSE. !< flag indicating whether u is already initialized from file
164       LOGICAL ::  from_file_v      = .FALSE. !< flag indicating whether v is already initialized from file
165       LOGICAL ::  from_file_w      = .FALSE. !< flag indicating whether w is already initialized from file
166
167       REAL(wp) ::  fill_msoil       !< fill value for soil moisture
168       REAL(wp) ::  fill_pt          !< fill value for pt
169       REAL(wp) ::  fill_q           !< fill value for q
170       REAL(wp) ::  fill_tsoil       !< fill value for soil temperature
171       REAL(wp) ::  fill_u           !< fill value for u
172       REAL(wp) ::  fill_v           !< fill value for v
173       REAL(wp) ::  fill_w           !< fill value for w
174       REAL(wp) ::  latitude         !< latitude of the southern model boundary
175       REAL(wp) ::  longitude        !< longitude of the western model boundary
176
177       REAL(wp), DIMENSION(:), ALLOCATABLE ::  msoil_init   !< initial vertical profile of soil moisture
178       REAL(wp), DIMENSION(:), ALLOCATABLE ::  pt_init      !< initial vertical profile of pt
179       REAL(wp), DIMENSION(:), ALLOCATABLE ::  q_init       !< initial vertical profile of q
180       REAL(wp), DIMENSION(:), ALLOCATABLE ::  tsoil_init   !< initial vertical profile of soil temperature
181       REAL(wp), DIMENSION(:), ALLOCATABLE ::  u_init       !< initial vertical profile of u
182       REAL(wp), DIMENSION(:), ALLOCATABLE ::  ug_init      !< initial vertical profile of ug
183       REAL(wp), DIMENSION(:), ALLOCATABLE ::  v_init       !< initial vertical profile of v
184       REAL(wp), DIMENSION(:), ALLOCATABLE ::  vg_init      !< initial vertical profile of ug
185       REAL(wp), DIMENSION(:), ALLOCATABLE ::  w_init       !< initial vertical profile of w
186       REAL(wp), DIMENSION(:), ALLOCATABLE ::  z_soil       !< vertical levels in soil in dynamic input file, used for interpolation
187       REAL(wp), DIMENSION(:), ALLOCATABLE ::  zu_atmos     !< vertical levels at scalar grid in dynamic input file, used for interpolation
188       REAL(wp), DIMENSION(:), ALLOCATABLE ::  zw_atmos     !< vertical levels at w grid in dynamic input file, used for interpolation
189
190
191       REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  msoil        !< initial 3d soil moisture provide by Inifor and interpolated onto soil grid
192       REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  tsoil        !< initial 3d soil temperature provide by Inifor and interpolated onto soil grid
193
194    END TYPE init_type
195
196!
197!-- Define data structures for different input data types.
198!-- 8-bit Integer 2D
199    TYPE int_2d_8bit
200       INTEGER(KIND=1) ::  fill = -127                      !< fill value
201       INTEGER(KIND=1), DIMENSION(:,:), ALLOCATABLE ::  var !< respective variable
202       
203       LOGICAL ::  from_file = .FALSE.  !< flag indicating whether an input variable is available and read from file or default values are used 
204    END TYPE int_2d_8bit
205!
206!-- 32-bit Integer 2D
207    TYPE int_2d_32bit
208       INTEGER(iwp) ::  fill = -9999                      !< fill value
209       INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE ::  var  !< respective variable
210
211       LOGICAL ::  from_file = .FALSE. !< flag indicating whether an input variable is available and read from file or default values are used 
212    END TYPE int_2d_32bit
213
214!
215!-- Define data type to read 2D real variables
216    TYPE real_2d
217       LOGICAL ::  from_file = .FALSE.  !< flag indicating whether an input variable is available and read from file or default values are used 
218
219       REAL(wp) ::  fill = -9999.9_wp                !< fill value
220       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  var !< respective variable
221    END TYPE real_2d
222
223!
224!-- Define data type to read 2D real variables
225    TYPE real_3d
226       LOGICAL ::  from_file = .FALSE.  !< flag indicating whether an input variable is available and read from file or default values are used 
227
228       INTEGER(iwp) ::  nz   !< number of grid points along vertical dimension                     
229
230       REAL(wp) ::  fill = -9999.9_wp                  !< fill value
231       REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  var !< respective variable
232    END TYPE real_3d
233!
234!-- Define data structure where the dimension and type of the input depends
235!-- on the given level of detail.
236!-- For buildings, the input is either 2D float, or 3d byte.
237    TYPE build_in
238       INTEGER(iwp)    ::  lod = 1                               !< level of detail                 
239       INTEGER(KIND=1) ::  fill2 = -127                          !< fill value for lod = 2
240       INTEGER(iwp)    ::  nz                                    !< number of vertical layers in file
241       INTEGER(KIND=1), DIMENSION(:,:,:), ALLOCATABLE ::  var_3d !< 3d variable (lod = 2)
242
243       REAL(wp), DIMENSION(:), ALLOCATABLE ::  z                 !< vertical coordinate for 3D building, used for consistency check
244
245       LOGICAL ::  from_file = .FALSE.  !< flag indicating whether an input variable is available and read from file or default values are used 
246
247       REAL(wp)                              ::  fill1 = -9999.9_wp !< fill values for lod = 1
248       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  var_2d             !< 2d variable (lod = 1)
249    END TYPE build_in
250
251!
252!-- For soil_type, the input is either 2D or 3D one-byte integer.
253    TYPE soil_in
254       INTEGER(iwp)                                   ::  lod = 1      !< level of detail                 
255       INTEGER(KIND=1)                                ::  fill = -127  !< fill value for lod = 2
256       INTEGER(KIND=1), DIMENSION(:,:), ALLOCATABLE   ::  var_2d       !< 2d variable (lod = 1)
257       INTEGER(KIND=1), DIMENSION(:,:,:), ALLOCATABLE ::  var_3d       !< 3d variable (lod = 2)
258
259       LOGICAL ::  from_file = .FALSE.  !< flag indicating whether an input variable is available and read from file or default values are used 
260    END TYPE soil_in
261
262!
263!-- Define data type for fractions between surface types
264    TYPE fracs
265       INTEGER(iwp)                            ::  nf             !< total number of fractions
266       INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  nfracs         !< dimension array for fraction
267
268       LOGICAL ::  from_file = .FALSE. !< flag indicating whether an input variable is available and read from file or default values are used 
269
270       REAL(wp)                                ::  fill = -9999.9_wp !< fill value
271       REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  frac              !< respective fraction between different surface types
272    END TYPE fracs
273!
274!-- Data type for parameter lists, Depending on the given level of detail,
275!-- the input is 3D or 4D
276    TYPE pars
277       INTEGER(iwp)                            ::  lod = 1         !< level of detail
278       INTEGER(iwp)                            ::  np              !< total number of parameters
279       INTEGER(iwp)                            ::  nz              !< vertical dimension - number of soil layers
280       INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  layers          !< dimension array for soil layers
281       INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  pars            !< dimension array for parameters
282
283       LOGICAL ::  from_file = .FALSE.  !< flag indicating whether an input variable is available and read from file or default values are used 
284
285       REAL(wp)                                  ::  fill = -9999.9_wp !< fill value
286       REAL(wp), DIMENSION(:,:,:), ALLOCATABLE   ::  pars_xy           !< respective parameters, level of detail = 1
287       REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE ::  pars_xyz          !< respective parameters, level of detail = 2
288    END TYPE pars
289
290    TYPE(force_type) ::  force !< input variable for lateral and top boundaries derived from large-scale model 
291
292    TYPE(init_type) ::  init_3d    !< data structure for the initialization of the 3D flow and soil fields
293    TYPE(init_type) ::  init_model !< data structure for the initialization of the model
294
295!
296!-- Define 2D variables of type NC_BYTE
297    TYPE(int_2d_8bit)  ::  albedo_type_f     !< input variable for albedo type
298    TYPE(int_2d_8bit)  ::  building_type_f   !< input variable for building type
299    TYPE(int_2d_8bit)  ::  pavement_type_f   !< input variable for pavenment type
300    TYPE(int_2d_8bit)  ::  street_crossing_f !< input variable for water type
301    TYPE(int_2d_8bit)  ::  street_type_f     !< input variable for water type
302    TYPE(int_2d_8bit)  ::  vegetation_type_f !< input variable for vegetation type
303    TYPE(int_2d_8bit)  ::  water_type_f      !< input variable for water type
304
305!
306!-- Define 2D variables of type NC_INT
307    TYPE(int_2d_32bit) ::  building_id_f     !< input variable for building ID
308!
309!-- Define 2D variables of type NC_FLOAT
310    TYPE(real_2d) ::  terrain_height_f       !< input variable for terrain height
311!
312!-- Define 3D variables of type NC_FLOAT
313    TYPE(real_3d) ::  basal_area_density_f    !< input variable for basal area density - resolved vegetation
314    TYPE(real_3d) ::  leaf_area_density_f     !< input variable for leaf area density - resolved vegetation
315    TYPE(real_3d) ::  root_area_density_lad_f !< input variable for root area density - resolved vegetation
316    TYPE(real_3d) ::  root_area_density_lsm_f !< input variable for root area density - parametrized vegetation
317
318!
319!-- Define input variable for buildings
320    TYPE(build_in) ::  buildings_f           !< input variable for buildings
321!
322!-- Define input variables for soil_type
323    TYPE(soil_in)  ::  soil_type_f           !< input variable for soil type
324
325    TYPE(fracs) ::  surface_fraction_f       !< input variable for surface fraction
326
327    TYPE(pars)  ::  albedo_pars_f              !< input variable for albedo parameters
328    TYPE(pars)  ::  building_pars_f            !< input variable for building parameters
329    TYPE(pars)  ::  pavement_pars_f            !< input variable for pavement parameters
330    TYPE(pars)  ::  pavement_subsurface_pars_f !< input variable for pavement parameters
331    TYPE(pars)  ::  soil_pars_f                !< input variable for soil parameters
332    TYPE(pars)  ::  vegetation_pars_f          !< input variable for vegetation parameters
333    TYPE(pars)  ::  water_pars_f               !< input variable for water parameters
334
335
336    CHARACTER(LEN=3)  ::  char_lod  = 'lod'         !< name of level-of-detail attribute in NetCDF file
337
338    CHARACTER(LEN=10) ::  char_fill = '_FillValue'  !< name of fill value attribute in NetCDF file
339    CHARACTER(LEN=10) ::  char_lon  = 'origin_lon'  !< name of global attribute for longitude in NetCDF file
340    CHARACTER(LEN=10) ::  char_lat  = 'origin_lat'  !< name of global attribute for latitude in NetCDF file
341
342    CHARACTER(LEN=100) ::  input_file_static  = 'PIDS_STATIC'  !< Name of file which comprises static input data
343    CHARACTER(LEN=100) ::  input_file_dynamic = 'PIDS_DYNAMIC' !< Name of file which comprises dynamic input data
344
345    INTEGER(iwp) ::  nc_stat         !< return value of nf90 function call
346
347    LOGICAL ::  input_pids_static  = .FALSE.   !< Flag indicating whether Palm-input-data-standard file containing static information exists
348    LOGICAL ::  input_pids_dynamic = .FALSE.   !< Flag indicating whether Palm-input-data-standard file containing dynamic information exists
349
350    SAVE
351
352    PRIVATE
353
354    INTERFACE netcdf_data_input_interpolate
355       MODULE PROCEDURE netcdf_data_input_interpolate_1d
356       MODULE PROCEDURE netcdf_data_input_interpolate_1d_soil
357       MODULE PROCEDURE netcdf_data_input_interpolate_2d
358       MODULE PROCEDURE netcdf_data_input_interpolate_3d
359    END INTERFACE netcdf_data_input_interpolate
360
361    INTERFACE netcdf_data_input_check_dynamic
362       MODULE PROCEDURE netcdf_data_input_check_dynamic
363    END INTERFACE netcdf_data_input_check_dynamic
364
365    INTERFACE netcdf_data_input_check_static
366       MODULE PROCEDURE netcdf_data_input_check_static
367    END INTERFACE netcdf_data_input_check_static
368
369    INTERFACE netcdf_data_input_inquire_file
370       MODULE PROCEDURE netcdf_data_input_inquire_file
371    END INTERFACE netcdf_data_input_inquire_file
372
373    INTERFACE netcdf_data_input_init
374       MODULE PROCEDURE netcdf_data_input_init
375    END INTERFACE netcdf_data_input_init
376
377    INTERFACE netcdf_data_input_init_3d
378       MODULE PROCEDURE netcdf_data_input_init_3d
379    END INTERFACE netcdf_data_input_init_3d
380
381    INTERFACE netcdf_data_input_lsf
382       MODULE PROCEDURE netcdf_data_input_lsf
383    END INTERFACE netcdf_data_input_lsf
384
385    INTERFACE netcdf_data_input_surface_data
386       MODULE PROCEDURE netcdf_data_input_surface_data
387    END INTERFACE netcdf_data_input_surface_data
388
389    INTERFACE netcdf_data_input_topo
390       MODULE PROCEDURE netcdf_data_input_topo
391    END INTERFACE netcdf_data_input_topo
392
393    INTERFACE get_variable
394       MODULE PROCEDURE get_variable_1d_int
395       MODULE PROCEDURE get_variable_1d_real
396       MODULE PROCEDURE get_variable_2d_int8
397       MODULE PROCEDURE get_variable_2d_int32
398       MODULE PROCEDURE get_variable_2d_real
399       MODULE PROCEDURE get_variable_3d_int8
400       MODULE PROCEDURE get_variable_3d_real
401       MODULE PROCEDURE get_variable_4d_real
402    END INTERFACE get_variable
403
404    INTERFACE get_attribute
405       MODULE PROCEDURE get_attribute_real
406       MODULE PROCEDURE get_attribute_int8
407       MODULE PROCEDURE get_attribute_int32
408       MODULE PROCEDURE get_attribute_string
409    END INTERFACE get_attribute
410
411!
412!-- Public variables
413    PUBLIC albedo_pars_f, albedo_type_f, basal_area_density_f, buildings_f,    &
414           building_id_f, building_pars_f, building_type_f, force, init_3d,    &
415           init_model, input_file_static, input_pids_static,                   &
416           input_pids_dynamic, leaf_area_density_f,                            &
417           pavement_pars_f, pavement_subsurface_pars_f, pavement_type_f,       &
418           root_area_density_lad_f, root_area_density_lsm_f, soil_pars_f,      &
419           soil_type_f, street_crossing_f, street_type_f, surface_fraction_f,  &
420           terrain_height_f, vegetation_pars_f, vegetation_type_f,             &
421           water_pars_f, water_type_f
422
423!
424!-- Public subroutines
425    PUBLIC netcdf_data_input_check_dynamic, netcdf_data_input_check_static,    &
426           netcdf_data_input_inquire_file,                                     &
427           netcdf_data_input_init, netcdf_data_input_init_3d,                  &
428           netcdf_data_input_interpolate, netcdf_data_input_lsf,               &
429           netcdf_data_input_surface_data, netcdf_data_input_topo
430
431 CONTAINS
432
433!------------------------------------------------------------------------------!
434! Description:
435! ------------
436!> Inquires whether NetCDF input files according to Palm-input-data standard
437!> exist. Moreover, basic checks are performed.
438!------------------------------------------------------------------------------!
439    SUBROUTINE netcdf_data_input_inquire_file
440
441       USE control_parameters,                                                 &
442           ONLY:  land_surface, message_string, topo_no_distinct, urban_surface
443
444       IMPLICIT NONE
445
446       LOGICAL ::  check_nest  !< flag indicating whether a check passed or not
447
448#if defined ( __netcdf )
449       INQUIRE( FILE = TRIM( input_file_static ) // TRIM( coupling_char ),     &
450                EXIST = input_pids_static  )
451       INQUIRE( FILE = TRIM( input_file_dynamic ) // TRIM( coupling_char ),    &
452                EXIST = input_pids_dynamic )
453#endif
454
455!
456!--    As long as topography can be input via ASCII format, no distinction
457!--    between building and terrain can be made. This case, classify all
458!--    surfaces as default type. Same in case land-surface and urban-surface
459!--    model are not applied.
460       IF ( .NOT. input_pids_static )  THEN
461          topo_no_distinct = .TRUE. 
462       ENDIF
463
464    END SUBROUTINE netcdf_data_input_inquire_file
465
466!------------------------------------------------------------------------------!
467! Description:
468! ------------
469!> Reads global attributes required for initialization of the model.
470!------------------------------------------------------------------------------!
471    SUBROUTINE netcdf_data_input_init
472
473       IMPLICIT NONE
474
475       INTEGER(iwp) ::  id_mod   !< NetCDF id of input file
476       INTEGER(iwp) ::  ii       !< running index for IO blocks
477
478       IF ( .NOT. input_pids_static )  RETURN 
479
480       DO  ii = 0, io_blocks-1
481          IF ( ii == io_group )  THEN
482#if defined ( __netcdf )
483!
484!--          Open file in read-only mode
485             CALL open_read_file( TRIM( input_file_static ) //                 &
486                                  TRIM( coupling_char ), id_mod )
487!
488!--          Read global attribute for latitude and longitude
489             CALL get_attribute( id_mod, char_lat,                             &
490                                 init_model%latitude, .TRUE. )
491
492             CALL get_attribute( id_mod, char_lon,                             &
493                                 init_model%longitude, .TRUE. )
494!
495!--          Finally, close input file
496             CALL close_input_file( id_mod )
497#endif
498          ENDIF
499#if defined( __parallel )
500          CALL MPI_BARRIER( comm2d, ierr )
501#endif
502       ENDDO
503
504    END SUBROUTINE netcdf_data_input_init
505
506!------------------------------------------------------------------------------!
507! Description:
508! ------------
509!> Reads surface classification data, such as vegetation and soil type, etc. .
510!------------------------------------------------------------------------------!
511    SUBROUTINE netcdf_data_input_surface_data
512
513       USE control_parameters,                                                 &
514           ONLY:  bc_lr_cyc, bc_ns_cyc, land_surface, message_string,          &
515                  plant_canopy, urban_surface
516
517       USE indices,                                                            &
518           ONLY:  nbgp, nx, nxl, nxlg, nxr, nxrg, ny, nyn, nyng, nys, nysg
519
520
521       IMPLICIT NONE
522
523       CHARACTER(LEN=100), DIMENSION(:), ALLOCATABLE ::  var_names  !< variable names in static input file
524
525       INTEGER(iwp) ::  i         !< running index along x-direction
526       INTEGER(iwp) ::  ii        !< running index for IO blocks
527       INTEGER(iwp) ::  id_surf   !< NetCDF id of input file
528       INTEGER(iwp) ::  j         !< running index along y-direction
529       INTEGER(iwp) ::  k         !< running index along z-direction
530       INTEGER(iwp) ::  k2        !< running index
531       INTEGER(iwp) ::  num_vars  !< number of variables in input file
532       INTEGER(iwp) ::  nz_soil   !< number of soil layers in file
533
534       INTEGER(iwp), DIMENSION(nysg:nyng,nxlg:nxrg) ::  var_exchange_int !< dummy variables used to exchange 32-bit Integer arrays
535       INTEGER(iwp), DIMENSION(:,:,:), ALLOCATABLE  ::  var_dum_int_3d !< dummy variables used to exchange real arrays
536
537       REAL(wp), DIMENSION(nysg:nyng,nxlg:nxrg) ::  var_exchange_real !< dummy variables used to exchange real arrays
538
539       REAL(wp), DIMENSION(:,:,:),   ALLOCATABLE ::  var_dum_real_3d !< dummy variables used to exchange real arrays
540       REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE ::  var_dum_real_4d !< dummy variables used to exchange real arrays
541
542!
543!--    If not static input file is available, skip this routine
544       IF ( .NOT. input_pids_static )  RETURN
545!
546!--    Read plant canopy variables.
547       IF ( plant_canopy )  THEN
548          DO  ii = 0, io_blocks-1
549             IF ( ii == io_group )  THEN
550#if defined ( __netcdf )
551!
552!--             Open file in read-only mode
553                CALL open_read_file( TRIM( input_file_static ) //              &
554                                     TRIM( coupling_char ) , id_surf )
555!
556!--             At first, inquire all variable names.
557!--             This will be used to check whether an optional input variable
558!--             exist or not.
559                CALL inquire_num_variables( id_surf, num_vars )
560
561                ALLOCATE( var_names(1:num_vars) )
562                CALL inquire_variable_names( id_surf, var_names )
563
564!
565!--             Read leaf area density - resolved vegetation
566                IF ( check_existence( var_names, 'leaf_area_density' ) )  THEN
567                   leaf_area_density_f%from_file = .TRUE. 
568                   CALL get_attribute( id_surf, char_fill,                     &
569                                       leaf_area_density_f%fill,               &
570                                       .FALSE., 'leaf_area_density' ) 
571!
572!--                Inquire number of vertical vegetation layer
573                   CALL get_dimension_length( id_surf, leaf_area_density_f%nz, &
574                                              'zlad' )
575!           
576!--                Allocate variable for leaf-area density
577                   ALLOCATE( leaf_area_density_f%var(                          &
578                                                   0:leaf_area_density_f%nz-1, &
579                                                   nys:nyn,nxl:nxr) )
580
581                   DO  i = nxl, nxr
582                      DO  j = nys, nyn
583                         CALL get_variable( id_surf, 'leaf_area_density',      &
584                                            i, j,                              &
585                                            leaf_area_density_f%var(:,j,i) ) 
586                      ENDDO
587                   ENDDO
588                ELSE
589                   leaf_area_density_f%from_file = .FALSE. 
590                ENDIF
591
592!
593!--             Read basal area density - resolved vegetation
594                IF ( check_existence( var_names, 'basal_area_density' ) )  THEN
595                   basal_area_density_f%from_file = .TRUE. 
596                   CALL get_attribute( id_surf, char_fill,                     &
597                                       basal_area_density_f%fill,              &
598                                       .FALSE., 'basal_area_density' ) 
599!
600!--                Inquire number of vertical vegetation layer
601                   CALL get_dimension_length( id_surf,                         &
602                                              basal_area_density_f%nz,         &
603                                              'zlad' )
604!           
605!--                Allocate variable
606                   ALLOCATE( basal_area_density_f%var(                         &
607                                                  0:basal_area_density_f%nz-1, &
608                                                  nys:nyn,nxl:nxr) )
609
610                   DO  i = nxl, nxr
611                      DO  j = nys, nyn
612                         CALL get_variable( id_surf, 'basal_area_density',     &
613                                            i, j,                              &
614                                            basal_area_density_f%var(:,j,i) ) 
615                      ENDDO
616                   ENDDO
617                ELSE
618                   basal_area_density_f%from_file = .FALSE. 
619                ENDIF
620
621!
622!--             Read root area density - resolved vegetation
623                IF ( check_existence( var_names, 'root_area_density_lad' ) )  THEN
624                   root_area_density_lad_f%from_file = .TRUE. 
625                   CALL get_attribute( id_surf, char_fill,                     &
626                                       root_area_density_lad_f%fill,           &
627                                       .FALSE., 'root_area_density_lad' ) 
628!
629!--                Inquire number of vertical soil layers
630                   CALL get_dimension_length( id_surf,                         &
631                                              root_area_density_lad_f%nz,      &
632                                              'zsoil' )
633!           
634!--                Allocate variable
635                   ALLOCATE( root_area_density_lad_f%var                       &
636                                               (0:root_area_density_lad_f%nz-1,&
637                                                nys:nyn,nxl:nxr) )
638
639                   DO  i = nxl, nxr
640                      DO  j = nys, nyn
641                         CALL get_variable( id_surf, 'root_area_density_lad',  &
642                                            i, j,                              &
643                                            root_area_density_lad_f%var(:,j,i) ) 
644                      ENDDO
645                   ENDDO
646                ELSE
647                   root_area_density_lad_f%from_file = .FALSE. 
648                ENDIF
649!
650!--             Finally, close input file
651                CALL close_input_file( id_surf )
652#endif
653             ENDIF
654#if defined( __parallel )
655             CALL MPI_BARRIER( comm2d, ierr )
656#endif
657          ENDDO
658!
659!--       Deallocate variable list. Will be re-allocated in case further
660!--       variables are read from file.         
661          IF ( ALLOCATED( var_names ) )  DEALLOCATE( var_names )
662
663       ENDIF
664!
665!--    Skip the following if no land-surface or urban-surface module are
666!--    applied. This case, no one of the following variables is used anyway. 
667       IF (  .NOT. land_surface  .OR.  .NOT. urban_surface )  RETURN
668!
669!--    Initialize dummy arrays used for ghost-point exchange
670       var_exchange_int  = 0
671       var_exchange_real = 0.0_wp
672
673       DO  ii = 0, io_blocks-1
674          IF ( ii == io_group )  THEN
675#if defined ( __netcdf )
676!
677!--          Open file in read-only mode
678             CALL open_read_file( TRIM( input_file_static ) //                 &
679                                  TRIM( coupling_char ) , id_surf )
680
681!
682!--          Inquire all variable names.
683!--          This will be used to check whether an optional input variable exist
684!--          or not.
685             CALL inquire_num_variables( id_surf, num_vars )
686
687             ALLOCATE( var_names(1:num_vars) )
688             CALL inquire_variable_names( id_surf, var_names )
689
690!
691!--          Read vegetation type and required attributes
692             IF ( check_existence( var_names, 'vegetation_type' ) )  THEN
693                vegetation_type_f%from_file = .TRUE.
694                CALL get_attribute( id_surf, char_fill,                        &
695                                    vegetation_type_f%fill,                    &
696                                    .FALSE., 'vegetation_type' )
697!
698!--             PE-wise reading of 2D vegetation type.
699                ALLOCATE ( vegetation_type_f%var(nys:nyn,nxl:nxr)  )
700
701                DO  i = nxl, nxr
702                   CALL get_variable( id_surf, 'vegetation_type',              &
703                                      i, vegetation_type_f%var(:,i) )
704                ENDDO
705             ELSE
706                vegetation_type_f%from_file = .FALSE.
707             ENDIF
708
709!
710!--          Read soil type and required attributes
711             IF ( check_existence( var_names, 'soil_type' ) )  THEN
712                   soil_type_f%from_file = .TRUE. 
713!
714!--             Note, lod is currently not on file; skip for the moment
715!                 CALL get_attribute( id_surf, char_lod,                       &
716!                                            soil_type_f%lod,                  &
717!                                            .FALSE., 'soil_type' )
718                CALL get_attribute( id_surf, char_fill,                        &
719                                    soil_type_f%fill,                          &
720                                    .FALSE., 'soil_type' ) 
721
722                IF ( soil_type_f%lod == 1 )  THEN
723!
724!--                PE-wise reading of 2D soil type.
725                   ALLOCATE ( soil_type_f%var_2d(nys:nyn,nxl:nxr)  )
726                   DO  i = nxl, nxr
727                      CALL get_variable( id_surf, 'soil_type',                 &
728                                         i, soil_type_f%var_2d(:,i) ) 
729                   ENDDO
730                ELSEIF ( soil_type_f%lod == 2 )  THEN
731!
732!--                Obtain number of soil layers from file.
733                   CALL get_dimension_length( id_surf, nz_soil, 'zsoil' )
734!
735!--                PE-wise reading of 3D soil type.
736                   ALLOCATE ( soil_type_f%var_3d(0:nz_soil,nys:nyn,nxl:nxr) )
737                   DO  i = nxl, nxr
738                      DO  j = nys, nyn
739                         CALL get_variable( id_surf, 'soil_type', i, j,        &
740                                            soil_type_f%var_3d(:,j,i) )   
741                      ENDDO
742                   ENDDO
743                ENDIF
744             ELSE
745                soil_type_f%from_file = .FALSE.
746             ENDIF
747
748!
749!--          Read pavement type and required attributes
750             IF ( check_existence( var_names, 'pavement_type' ) )  THEN
751                pavement_type_f%from_file = .TRUE. 
752                CALL get_attribute( id_surf, char_fill,                        &
753                                    pavement_type_f%fill, .FALSE.,             &
754                                    'pavement_type' ) 
755!
756!--             PE-wise reading of 2D pavement type.
757                ALLOCATE ( pavement_type_f%var(nys:nyn,nxl:nxr)  )
758                DO  i = nxl, nxr
759                   CALL get_variable( id_surf, 'pavement_type',                &
760                                      i, pavement_type_f%var(:,i) ) 
761                ENDDO
762             ELSE
763                pavement_type_f%from_file = .FALSE.
764             ENDIF
765
766!
767!--          Read water type and required attributes
768             IF ( check_existence( var_names, 'water_type' ) )  THEN
769                water_type_f%from_file = .TRUE. 
770                CALL get_attribute( id_surf, char_fill, water_type_f%fill,     &
771                                    .FALSE., 'water_type' )
772!
773!--             PE-wise reading of 2D water type.
774                ALLOCATE ( water_type_f%var(nys:nyn,nxl:nxr)  )
775                DO  i = nxl, nxr
776                   CALL get_variable( id_surf, 'water_type', i,                &
777                                      water_type_f%var(:,i) ) 
778                ENDDO
779             ELSE
780                water_type_f%from_file = .FALSE.
781             ENDIF
782!
783!--          Read surface fractions and related information
784             IF ( check_existence( var_names, 'surface_fraction' ) )  THEN
785                surface_fraction_f%from_file = .TRUE. 
786                CALL get_attribute( id_surf, char_fill,                        &
787                                    surface_fraction_f%fill,                   &
788                                    .FALSE., 'surface_fraction' )
789!
790!--             Inquire number of surface fractions
791                CALL get_dimension_length( id_surf,                            &
792                                           surface_fraction_f%nf,              &
793                                           'nsurface_fraction' )
794!           
795!--             Allocate dimension array and input array for surface fractions
796                ALLOCATE( surface_fraction_f%nfracs(0:surface_fraction_f%nf-1) )
797                ALLOCATE( surface_fraction_f%frac(0:surface_fraction_f%nf-1,   &
798                                                  nys:nyn,nxl:nxr) )
799!
800!--             Get dimension of surface fractions
801                CALL get_variable( id_surf, 'nsurface_fraction',               &
802                                   surface_fraction_f%nfracs )
803!
804!--             Read surface fractions
805                DO  i = nxl, nxr
806                   DO  j = nys, nyn
807                      CALL get_variable( id_surf, 'surface_fraction', i, j,    &
808                                         surface_fraction_f%frac(:,j,i) ) 
809                   ENDDO
810                ENDDO
811             ELSE
812                surface_fraction_f%from_file = .FALSE. 
813             ENDIF
814!
815!--          Read building parameters and related information
816             IF ( check_existence( var_names, 'building_pars' ) )  THEN
817                building_pars_f%from_file = .TRUE. 
818                CALL get_attribute( id_surf, char_fill,                        &
819                                    building_pars_f%fill,                      &
820                                    .FALSE., 'building_pars' ) 
821!
822!--             Inquire number of building parameters
823                CALL get_dimension_length( id_surf,                            &
824                                           building_pars_f%np,                 &
825                                           'nbuilding_pars' )
826!           
827!--             Allocate dimension array and input array for building parameters
828                ALLOCATE( building_pars_f%pars(0:building_pars_f%np-1) )
829                ALLOCATE( building_pars_f%pars_xy(0:building_pars_f%np-1,      &
830                                                  nys:nyn,nxl:nxr) )
831!
832!--             Get dimension of building parameters
833                CALL get_variable( id_surf, 'nbuilding_pars',                  &
834                                   building_pars_f%pars )
835!
836!--             Read building_pars
837                DO  i = nxl, nxr
838                   DO  j = nys, nyn
839                      CALL get_variable( id_surf, 'building_pars', i, j,       &
840                                         building_pars_f%pars_xy(:,j,i) )   
841                   ENDDO
842                ENDDO
843             ELSE
844                building_pars_f%from_file = .FALSE. 
845             ENDIF
846
847!
848!--          Read albedo type and required attributes
849             IF ( check_existence( var_names, 'albedo_type' ) )  THEN
850                albedo_type_f%from_file = .TRUE. 
851                CALL get_attribute( id_surf, char_fill, albedo_type_f%fill,    &
852                                    .FALSE.,  'albedo_type' ) 
853!
854!--             PE-wise reading of 2D water type.
855                ALLOCATE ( albedo_type_f%var(nys:nyn,nxl:nxr)  )
856                DO  i = nxl, nxr
857                   CALL get_variable( id_surf, 'albedo_type',                  &
858                                      i, albedo_type_f%var(:,i) ) 
859                ENDDO
860             ELSE
861                albedo_type_f%from_file = .FALSE.
862             ENDIF
863!
864!--          Read albedo parameters and related information
865             IF ( check_existence( var_names, 'albedo_pars' ) )  THEN
866                albedo_pars_f%from_file = .TRUE. 
867                CALL get_attribute( id_surf, char_fill, albedo_pars_f%fill,    &
868                                    .FALSE., 'albedo_pars' ) 
869!
870!--             Inquire number of albedo parameters
871                CALL get_dimension_length( id_surf, albedo_pars_f%np,          &
872                                           'nalbedo_pars' )
873!           
874!--             Allocate dimension array and input array for albedo parameters
875                ALLOCATE( albedo_pars_f%pars(0:albedo_pars_f%np-1) )
876                ALLOCATE( albedo_pars_f%pars_xy(0:albedo_pars_f%np-1,          &
877                                                nys:nyn,nxl:nxr) )
878!
879!--             Get dimension of albedo parameters
880                CALL get_variable( id_surf, 'nalbedo_pars', albedo_pars_f%pars )
881
882                DO  i = nxl, nxr
883                   DO  j = nys, nyn
884                      CALL get_variable( id_surf, 'albedo_pars', i, j,         &
885                                         albedo_pars_f%pars_xy(:,j,i) )   
886                   ENDDO
887                ENDDO
888             ELSE
889                albedo_pars_f%from_file = .FALSE. 
890             ENDIF
891
892!
893!--          Read pavement parameters and related information
894             IF ( check_existence( var_names, 'pavement_pars' ) )  THEN
895                pavement_pars_f%from_file = .TRUE. 
896                CALL get_attribute( id_surf, char_fill,                        &
897                                    pavement_pars_f%fill,                      &
898                                    .FALSE., 'pavement_pars' ) 
899!
900!--             Inquire number of pavement parameters
901                CALL get_dimension_length( id_surf, pavement_pars_f%np,        &
902                                           'npavement_pars' )
903!           
904!--             Allocate dimension array and input array for pavement parameters
905                ALLOCATE( pavement_pars_f%pars(0:pavement_pars_f%np-1) )
906                ALLOCATE( pavement_pars_f%pars_xy(0:pavement_pars_f%np-1,      &
907                                                  nys:nyn,nxl:nxr) )
908!
909!--             Get dimension of pavement parameters
910                CALL get_variable( id_surf, 'npavement_pars',                  &
911                                   pavement_pars_f%pars )
912   
913                DO  i = nxl, nxr
914                   DO  j = nys, nyn
915                      CALL get_variable( id_surf, 'pavement_pars', i, j,       &
916                                         pavement_pars_f%pars_xy(:,j,i) )   
917                   ENDDO
918                ENDDO
919             ELSE
920                pavement_pars_f%from_file = .FALSE. 
921             ENDIF
922
923!
924!--          Read pavement subsurface parameters and related information
925             IF ( check_existence( var_names, 'pavement_subsurface_pars' ) )   &
926             THEN
927                pavement_subsurface_pars_f%from_file = .TRUE. 
928                CALL get_attribute( id_surf, char_fill,                        &
929                                    pavement_subsurface_pars_f%fill,           &
930                                    .FALSE., 'pavement_subsurface_pars' ) 
931!
932!--             Inquire number of parameters
933                CALL get_dimension_length( id_surf,                            &
934                                           pavement_subsurface_pars_f%np,      &
935                                           'npavement_subsurface_pars' )
936!
937!--             Inquire number of soil layers
938                CALL get_dimension_length( id_surf,                            &
939                                           pavement_subsurface_pars_f%nz,      &
940                                           'zsoil' )
941!           
942!--             Allocate dimension array and input array for pavement parameters
943                ALLOCATE( pavement_subsurface_pars_f%pars                      &
944                                  (0:pavement_subsurface_pars_f%np-1) )
945                ALLOCATE( pavement_subsurface_pars_f%pars_xyz                  &
946                                  (0:pavement_subsurface_pars_f%np-1,          &
947                                   0:pavement_subsurface_pars_f%nz-1,          &
948                                   nys:nyn,nxl:nxr) )
949!
950!--             Get dimension of pavement parameters
951                CALL get_variable( id_surf, 'npavement_subsurface_pars',       &
952                                   pavement_subsurface_pars_f%pars )
953   
954                DO  i = nxl, nxr
955                   DO  j = nys, nyn
956                      CALL get_variable(                                       &
957                                  id_surf, 'pavement_subsurface_pars',         &
958                                  i, j,                                        &
959                                  pavement_subsurface_pars_f%pars_xyz(:,:,j,i),&
960                                  pavement_subsurface_pars_f%nz,               &
961                                  pavement_subsurface_pars_f%np )   
962                   ENDDO
963                ENDDO
964             ELSE
965                pavement_subsurface_pars_f%from_file = .FALSE. 
966             ENDIF
967
968
969!
970!--          Read vegetation parameters and related information
971             IF ( check_existence( var_names, 'vegetation_pars' ) )  THEN
972                vegetation_pars_f%from_file = .TRUE. 
973                CALL get_attribute( id_surf, char_fill,                        &
974                                    vegetation_pars_f%fill,                    &
975                                    .FALSE.,  'vegetation_pars' )
976!
977!--             Inquire number of vegetation parameters
978                CALL get_dimension_length( id_surf, vegetation_pars_f%np,      &
979                                           'nvegetation_pars' )
980!           
981!--             Allocate dimension array and input array for surface fractions
982                ALLOCATE( vegetation_pars_f%pars(0:vegetation_pars_f%np-1) )
983                ALLOCATE( vegetation_pars_f%pars_xy(0:vegetation_pars_f%np-1,  &
984                                                    nys:nyn,nxl:nxr) )
985!
986!--             Get dimension of the parameters
987                CALL get_variable( id_surf, 'nvegetation_pars',                &
988                                   vegetation_pars_f%pars )
989
990                DO  i = nxl, nxr
991                   DO  j = nys, nyn
992                      CALL get_variable( id_surf, 'vegetation_pars', i, j,     &
993                                         vegetation_pars_f%pars_xy(:,j,i) ) 
994                   ENDDO
995                ENDDO
996             ELSE
997                vegetation_pars_f%from_file = .FALSE. 
998             ENDIF
999
1000!
1001!--          Read root parameters/distribution and related information
1002             IF ( check_existence( var_names, 'soil_pars' ) )  THEN
1003                soil_pars_f%from_file = .TRUE. 
1004                CALL get_attribute( id_surf, char_fill,                        &
1005                                    soil_pars_f%fill,                          & 
1006                                    .FALSE., 'soil_pars' ) 
1007
1008                CALL get_attribute( id_surf, char_lod,                         &
1009                                    soil_pars_f%lod,                           &
1010                                    .FALSE., 'soil_pars' ) 
1011
1012!
1013!--             Inquire number of soil parameters
1014                CALL get_dimension_length( id_surf,                            &
1015                                           soil_pars_f%np,                     &
1016                                           'nsoil_pars' )
1017!
1018!--             Read parameters array
1019                ALLOCATE( soil_pars_f%pars(0:soil_pars_f%np-1) )
1020                CALL get_variable( id_surf, 'nsoil_pars', soil_pars_f%pars )
1021
1022!
1023!--             In case of level of detail 2, also inquire number of vertical
1024!--             soil layers, allocate memory and read the respective dimension
1025                IF ( soil_pars_f%lod == 2 )  THEN
1026                   CALL get_dimension_length( id_surf, soil_pars_f%nz, 'zsoil' )
1027
1028                   ALLOCATE( soil_pars_f%layers(0:soil_pars_f%nz-1) )
1029                   CALL get_variable( id_surf, 'zsoil', soil_pars_f%layers )
1030
1031                ENDIF
1032
1033!
1034!--             Read soil parameters, depending on level of detail
1035                IF ( soil_pars_f%lod == 1 )  THEN     
1036                   ALLOCATE( soil_pars_f%pars_xy(0:soil_pars_f%np-1,           &
1037                                                 nys:nyn,nxl:nxr) )       
1038                   DO  i = nxl, nxr
1039                      DO  j = nys, nyn
1040                         CALL get_variable( id_surf, 'soil_pars', i, j,        &
1041                                            soil_pars_f%pars_xy(:,j,i) ) 
1042                      ENDDO
1043                   ENDDO
1044
1045                ELSEIF ( soil_pars_f%lod == 2 )  THEN
1046                   ALLOCATE( soil_pars_f%pars_xyz(0:soil_pars_f%np-1,          &
1047                                                  0:soil_pars_f%nz-1,          &
1048                                                  nys:nyn,nxl:nxr) )
1049                   DO  i = nxl, nxr
1050                      DO  j = nys, nyn
1051                         CALL get_variable( id_surf, 'soil_pars', i, j,        &
1052                                            soil_pars_f%pars_xyz(:,:,j,i),     &
1053                                            soil_pars_f%nz, soil_pars_f%np ) 
1054                      ENDDO
1055                   ENDDO
1056                ENDIF
1057             ELSE
1058                soil_pars_f%from_file = .FALSE. 
1059             ENDIF
1060
1061!
1062!--          Read water parameters and related information
1063             IF ( check_existence( var_names, 'water_pars' ) )  THEN
1064                water_pars_f%from_file = .TRUE. 
1065                CALL get_attribute( id_surf, char_fill,                        &
1066                                    water_pars_f%fill,                         &
1067                                    .FALSE., 'water_pars' ) 
1068!
1069!--             Inquire number of water parameters
1070                CALL get_dimension_length( id_surf,                            & 
1071                                           water_pars_f%np,                    &
1072                                           'nwater_pars' )
1073!           
1074!--             Allocate dimension array and input array for water parameters
1075                ALLOCATE( water_pars_f%pars(0:water_pars_f%np-1) )
1076                ALLOCATE( water_pars_f%pars_xy(0:water_pars_f%np-1,            &
1077                                               nys:nyn,nxl:nxr) )
1078!
1079!--             Get dimension of water parameters
1080                CALL get_variable( id_surf, 'nwater_pars', water_pars_f%pars )
1081
1082                DO  i = nxl, nxr
1083                   DO  j = nys, nyn
1084                      CALL get_variable( id_surf, 'water_pars', i, j,          &
1085                                         water_pars_f%pars_xy(:,j,i) )   
1086                   ENDDO
1087                ENDDO
1088             ELSE
1089                water_pars_f%from_file = .FALSE. 
1090             ENDIF
1091!
1092!--          Read root area density - parametrized vegetation
1093             IF ( check_existence( var_names, 'root_area_density_lsm' ) )  THEN
1094                root_area_density_lsm_f%from_file = .TRUE.
1095                CALL get_attribute( id_surf, char_fill,                        &
1096                                    root_area_density_lsm_f%fill,              &
1097                                    .FALSE., 'root_area_density_lsm' )
1098!
1099!--             Obtain number of soil layers from file and allocate variable
1100                CALL get_dimension_length( id_surf, root_area_density_lsm_f%nz,&
1101                                           'zsoil' )
1102                ALLOCATE( root_area_density_lsm_f%var                          &
1103                                              (0:root_area_density_lsm_f%nz-1, &
1104                                               nys:nyn,nxl:nxr) )
1105
1106!
1107!--             Read root-area density
1108                DO  i = nxl, nxr
1109                   DO  j = nys, nyn
1110                      CALL get_variable( id_surf, 'root_area_density_lsm',     &
1111                                         i, j,                                 &
1112                                         root_area_density_lsm_f%var(:,j,i) )   
1113                   ENDDO
1114                ENDDO
1115
1116             ELSE
1117                root_area_density_lsm_f%from_file = .FALSE.
1118             ENDIF
1119!
1120!--          Read street type and street crossing
1121             IF ( check_existence( var_names, 'street_type' ) )  THEN
1122                street_type_f%from_file = .TRUE. 
1123                CALL get_attribute( id_surf, char_fill,                        &
1124                                    street_type_f%fill, .FALSE.,               &
1125                                    'street_type' ) 
1126!
1127!--             PE-wise reading of 2D pavement type.
1128                ALLOCATE ( street_type_f%var(nys:nyn,nxl:nxr)  )
1129                DO  i = nxl, nxr
1130                   CALL get_variable( id_surf, 'street_type',                  &
1131                                      i, street_type_f%var(:,i) ) 
1132                ENDDO
1133             ELSE
1134                street_type_f%from_file = .FALSE.
1135             ENDIF
1136
1137             IF ( check_existence( var_names, 'street_crossing' ) )  THEN
1138                street_crossing_f%from_file = .TRUE. 
1139                CALL get_attribute( id_surf, char_fill,                        &
1140                                    street_crossing_f%fill, .FALSE.,           &
1141                                    'street_crossing' ) 
1142!
1143!--             PE-wise reading of 2D pavement type.
1144                ALLOCATE ( street_crossing_f%var(nys:nyn,nxl:nxr)  )
1145                DO  i = nxl, nxr
1146                   CALL get_variable( id_surf, 'street_crossing',              &
1147                                      i, street_crossing_f%var(:,i) ) 
1148                ENDDO
1149             ELSE
1150                street_crossing_f%from_file = .FALSE.
1151             ENDIF
1152!
1153!--          Still missing: root_resolved and building_surface_pars.
1154!--          Will be implemented as soon as they are available.
1155
1156!
1157!--          Finally, close input file
1158             CALL close_input_file( id_surf )
1159#endif
1160          ENDIF
1161#if defined( __parallel )
1162          CALL MPI_BARRIER( comm2d, ierr )
1163#endif
1164       ENDDO
1165!
1166!--    Exchange 1 ghost points for surface variables. Please note, ghost point
1167!--    exchange for 3D parameter lists should be revised by using additional
1168!--    MPI datatypes or rewriting exchange_horiz.
1169!--    Moreover, varialbes will be resized in the following, including ghost
1170!--    points.
1171!--    Start with 2D Integer variables. Please note, for 8-bit integer
1172!--    variables must be swapt to 32-bit integer before calling exchange_horiz.
1173       IF ( albedo_type_f%from_file )  THEN
1174          var_exchange_int                  = INT( albedo_type_f%fill, KIND = 1 )
1175          var_exchange_int(nys:nyn,nxl:nxr) =                                  &
1176                            INT( albedo_type_f%var(nys:nyn,nxl:nxr), KIND = 4 )
1177          CALL exchange_horiz_2d_int( var_exchange_int, nys, nyn, nxl, nxr, nbgp )
1178          DEALLOCATE( albedo_type_f%var )
1179          ALLOCATE( albedo_type_f%var(nysg:nyng,nxlg:nxrg) )
1180          albedo_type_f%var = INT( var_exchange_int, KIND = 1 )
1181       ENDIF
1182       IF ( pavement_type_f%from_file )  THEN
1183          var_exchange_int                  = INT( pavement_type_f%fill, KIND = 1 )
1184          var_exchange_int(nys:nyn,nxl:nxr) =                                  &
1185                          INT( pavement_type_f%var(nys:nyn,nxl:nxr), KIND = 4 )
1186          CALL exchange_horiz_2d_int( var_exchange_int, nys, nyn, nxl, nxr, nbgp )
1187          DEALLOCATE( pavement_type_f%var )
1188          ALLOCATE( pavement_type_f%var(nysg:nyng,nxlg:nxrg) )
1189          pavement_type_f%var = INT( var_exchange_int, KIND = 1 )
1190       ENDIF
1191       IF ( soil_type_f%from_file  .AND.  ALLOCATED( soil_type_f%var_2d ) )  THEN
1192          var_exchange_int                  = INT( soil_type_f%fill, KIND = 1 )
1193          var_exchange_int(nys:nyn,nxl:nxr) =                                  &
1194                            INT( soil_type_f%var_2d(nys:nyn,nxl:nxr), KIND = 4 )
1195          CALL exchange_horiz_2d_int( var_exchange_int, nys, nyn, nxl, nxr, nbgp )
1196          DEALLOCATE( soil_type_f%var_2d )
1197          ALLOCATE( soil_type_f%var_2d(nysg:nyng,nxlg:nxrg) )
1198          soil_type_f%var_2d = INT( var_exchange_int, KIND = 1 )
1199       ENDIF
1200       IF ( vegetation_type_f%from_file )  THEN
1201          var_exchange_int                  = INT( vegetation_type_f%fill, KIND = 1 )
1202          var_exchange_int(nys:nyn,nxl:nxr) =                                  &
1203                        INT( vegetation_type_f%var(nys:nyn,nxl:nxr), KIND = 4 )
1204          CALL exchange_horiz_2d_int( var_exchange_int, nys, nyn, nxl, nxr, nbgp )
1205          DEALLOCATE( vegetation_type_f%var )
1206          ALLOCATE( vegetation_type_f%var(nysg:nyng,nxlg:nxrg) )
1207          vegetation_type_f%var = INT( var_exchange_int, KIND = 1 )
1208       ENDIF
1209       IF ( water_type_f%from_file )  THEN
1210          var_exchange_int                  = INT( water_type_f%fill, KIND = 1 )
1211          var_exchange_int(nys:nyn,nxl:nxr) =                                  &
1212                         INT( water_type_f%var(nys:nyn,nxl:nxr), KIND = 4 )
1213          CALL exchange_horiz_2d_int( var_exchange_int, nys, nyn, nxl, nxr, nbgp )
1214          DEALLOCATE( water_type_f%var )
1215          ALLOCATE( water_type_f%var(nysg:nyng,nxlg:nxrg) )
1216          water_type_f%var = INT( var_exchange_int, KIND = 1 )
1217       ENDIF
1218!
1219!--    Exchange 1 ghost point for 3/4-D variables. For the sake of simplicity,
1220!--    loop further dimensions to use 2D exchange routines.
1221!--    This should be revised later by introducing new MPI datatypes.
1222       IF ( soil_type_f%from_file  .AND.  ALLOCATED( soil_type_f%var_3d ) )    &
1223       THEN
1224          ALLOCATE( var_dum_int_3d(0:nz_soil,nys:nyn,nxl:nxr) ) 
1225          var_dum_int_3d = soil_type_f%var_3d
1226          DEALLOCATE( soil_type_f%var_3d )
1227          ALLOCATE( soil_type_f%var_3d(0:nz_soil,nysg:nyng,nxlg:nxrg) )
1228          soil_type_f%var_3d = soil_type_f%fill
1229
1230          DO  k = 0, nz_soil
1231             var_exchange_int(nys:nyn,nxl:nxr) = var_dum_int_3d(k,nys:nyn,nxl:nxr)
1232             CALL exchange_horiz_2d_int( var_exchange_int, nys, nyn, nxl, nxr, nbgp )
1233             soil_type_f%var_3d(k,:,:) = INT( var_exchange_int(:,:), KIND = 1 )
1234          ENDDO
1235          DEALLOCATE( var_dum_int_3d )
1236       ENDIF
1237
1238       IF ( surface_fraction_f%from_file )  THEN
1239          ALLOCATE( var_dum_real_3d(0:surface_fraction_f%nf-1,nys:nyn,nxl:nxr) )
1240          var_dum_real_3d = surface_fraction_f%frac
1241          DEALLOCATE( surface_fraction_f%frac )
1242          ALLOCATE( surface_fraction_f%frac(0:surface_fraction_f%nf-1,         &
1243                                            nysg:nyng,nxlg:nxrg) )
1244          surface_fraction_f%frac = surface_fraction_f%fill
1245
1246          DO  k = 0, surface_fraction_f%nf-1
1247             var_exchange_real(nys:nyn,nxl:nxr) = var_dum_real_3d(k,nys:nyn,nxl:nxr)
1248             CALL exchange_horiz_2d( var_exchange_real, nbgp )
1249             surface_fraction_f%frac(k,:,:) = var_exchange_real(:,:)
1250          ENDDO
1251          DEALLOCATE( var_dum_real_3d )
1252       ENDIF
1253
1254       IF ( building_pars_f%from_file )  THEN
1255          ALLOCATE( var_dum_real_3d(0:building_pars_f%np-1,nys:nyn,nxl:nxr) )
1256          var_dum_real_3d = building_pars_f%pars_xy
1257          DEALLOCATE( building_pars_f%pars_xy )
1258          ALLOCATE( building_pars_f%pars_xy(0:building_pars_f%np-1,            &
1259                                            nysg:nyng,nxlg:nxrg) )
1260          building_pars_f%pars_xy = building_pars_f%fill
1261          DO  k = 0, building_pars_f%np-1
1262             var_exchange_real(nys:nyn,nxl:nxr) =                              &
1263                                              var_dum_real_3d(k,nys:nyn,nxl:nxr)
1264             CALL exchange_horiz_2d( var_exchange_real, nbgp )
1265             building_pars_f%pars_xy(k,:,:) = var_exchange_real(:,:)
1266          ENDDO
1267          DEALLOCATE( var_dum_real_3d )
1268       ENDIF
1269
1270       IF ( albedo_pars_f%from_file )  THEN
1271          ALLOCATE( var_dum_real_3d(0:albedo_pars_f%np-1,nys:nyn,nxl:nxr) )
1272          var_dum_real_3d = albedo_pars_f%pars_xy
1273          DEALLOCATE( albedo_pars_f%pars_xy )
1274          ALLOCATE( albedo_pars_f%pars_xy(0:albedo_pars_f%np-1,                &
1275                                          nysg:nyng,nxlg:nxrg) )
1276          albedo_pars_f%pars_xy = albedo_pars_f%fill
1277          DO  k = 0, albedo_pars_f%np-1
1278             var_exchange_real(nys:nyn,nxl:nxr) =                              &
1279                                              var_dum_real_3d(k,nys:nyn,nxl:nxr)
1280             CALL exchange_horiz_2d( var_exchange_real, nbgp )
1281             albedo_pars_f%pars_xy(k,:,:) = var_exchange_real(:,:)
1282          ENDDO
1283          DEALLOCATE( var_dum_real_3d )
1284       ENDIF
1285
1286       IF ( pavement_pars_f%from_file )  THEN
1287          ALLOCATE( var_dum_real_3d(0:pavement_pars_f%np-1,nys:nyn,nxl:nxr) )
1288          var_dum_real_3d = pavement_pars_f%pars_xy
1289          DEALLOCATE( pavement_pars_f%pars_xy )
1290          ALLOCATE( pavement_pars_f%pars_xy(0:pavement_pars_f%np-1,            &
1291                                            nysg:nyng,nxlg:nxrg) )
1292          pavement_pars_f%pars_xy = pavement_pars_f%fill
1293          DO  k = 0, pavement_pars_f%np-1
1294             var_exchange_real(nys:nyn,nxl:nxr) =                              &
1295                                              var_dum_real_3d(k,nys:nyn,nxl:nxr)
1296             CALL exchange_horiz_2d( var_exchange_real, nbgp )
1297             pavement_pars_f%pars_xy(k,:,:) = var_exchange_real(:,:)
1298          ENDDO
1299          DEALLOCATE( var_dum_real_3d )
1300       ENDIF
1301
1302       IF ( vegetation_pars_f%from_file )  THEN
1303          ALLOCATE( var_dum_real_3d(0:vegetation_pars_f%np-1,nys:nyn,nxl:nxr) )
1304          var_dum_real_3d = vegetation_pars_f%pars_xy
1305          DEALLOCATE( vegetation_pars_f%pars_xy )
1306          ALLOCATE( vegetation_pars_f%pars_xy(0:vegetation_pars_f%np-1,        &
1307                                            nysg:nyng,nxlg:nxrg) )
1308          vegetation_pars_f%pars_xy = vegetation_pars_f%fill
1309          DO  k = 0, vegetation_pars_f%np-1
1310             var_exchange_real(nys:nyn,nxl:nxr) =                              &
1311                                              var_dum_real_3d(k,nys:nyn,nxl:nxr)
1312             CALL exchange_horiz_2d( var_exchange_real, nbgp )
1313             vegetation_pars_f%pars_xy(k,:,:) = var_exchange_real(:,:)
1314          ENDDO
1315          DEALLOCATE( var_dum_real_3d )
1316       ENDIF
1317
1318       IF ( water_pars_f%from_file )  THEN
1319          ALLOCATE( var_dum_real_3d(0:water_pars_f%np-1,nys:nyn,nxl:nxr) )
1320          var_dum_real_3d = water_pars_f%pars_xy
1321          DEALLOCATE( water_pars_f%pars_xy )
1322          ALLOCATE( water_pars_f%pars_xy(0:water_pars_f%np-1,                  &
1323                                         nysg:nyng,nxlg:nxrg) )
1324          water_pars_f%pars_xy = water_pars_f%fill
1325          DO  k = 0, water_pars_f%np-1
1326             var_exchange_real(nys:nyn,nxl:nxr) =                              &
1327                                              var_dum_real_3d(k,nys:nyn,nxl:nxr)
1328             CALL exchange_horiz_2d( var_exchange_real, nbgp )
1329             water_pars_f%pars_xy(k,:,:) = var_exchange_real(:,:)
1330          ENDDO
1331          DEALLOCATE( var_dum_real_3d )
1332       ENDIF
1333
1334       IF ( root_area_density_lsm_f%from_file )  THEN
1335          ALLOCATE( var_dum_real_3d(0:root_area_density_lsm_f%nz-1,nys:nyn,nxl:nxr) )
1336          var_dum_real_3d = root_area_density_lsm_f%var
1337          DEALLOCATE( root_area_density_lsm_f%var )
1338          ALLOCATE( root_area_density_lsm_f%var(0:root_area_density_lsm_f%nz-1,&
1339                                                nysg:nyng,nxlg:nxrg) )
1340          root_area_density_lsm_f%var = root_area_density_lsm_f%fill
1341
1342          DO  k = 0, root_area_density_lsm_f%nz-1
1343             var_exchange_real(nys:nyn,nxl:nxr) =                              &
1344                                              var_dum_real_3d(k,nys:nyn,nxl:nxr)
1345             CALL exchange_horiz_2d( var_exchange_real, nbgp )
1346             root_area_density_lsm_f%var(k,:,:) = var_exchange_real(:,:)
1347          ENDDO
1348          DEALLOCATE( var_dum_real_3d )
1349       ENDIF
1350
1351       IF ( soil_pars_f%from_file )  THEN
1352          IF ( soil_pars_f%lod == 1 )  THEN
1353
1354             ALLOCATE( var_dum_real_3d(0:soil_pars_f%np-1,nys:nyn,nxl:nxr) )
1355             var_dum_real_3d = soil_pars_f%pars_xy
1356             DEALLOCATE( soil_pars_f%pars_xy )
1357             ALLOCATE( soil_pars_f%pars_xy(0:soil_pars_f%np-1,                 &
1358                                            nysg:nyng,nxlg:nxrg) )
1359             soil_pars_f%pars_xy = soil_pars_f%fill
1360
1361             DO  k = 0, soil_pars_f%np-1
1362                var_exchange_real(nys:nyn,nxl:nxr) =                           &
1363                                              var_dum_real_3d(k,nys:nyn,nxl:nxr)
1364                CALL exchange_horiz_2d( var_exchange_real, nbgp )
1365                soil_pars_f%pars_xy(k,:,:) = var_exchange_real(:,:)
1366             ENDDO
1367             DEALLOCATE( var_dum_real_3d )
1368          ELSEIF ( soil_pars_f%lod == 2 )  THEN
1369             ALLOCATE( var_dum_real_4d(0:soil_pars_f%np-1,                     &
1370                                       0:soil_pars_f%nz-1,                     &
1371                                       nys:nyn,nxl:nxr) )
1372             var_dum_real_4d = soil_pars_f%pars_xyz
1373             DEALLOCATE( soil_pars_f%pars_xyz )
1374             ALLOCATE( soil_pars_f%pars_xyz(0:soil_pars_f%np-1,                &
1375                                            0:soil_pars_f%nz-1,                &
1376                                            nysg:nyng,nxlg:nxrg) )
1377             soil_pars_f%pars_xyz = soil_pars_f%fill
1378
1379             DO  k2 = 0, soil_pars_f%nz-1
1380                DO  k = 0, soil_pars_f%np-1
1381                   var_exchange_real(nys:nyn,nxl:nxr) =                        &
1382                                           var_dum_real_4d(k,k2,nys:nyn,nxl:nxr)
1383                   CALL exchange_horiz_2d( var_exchange_real, nbgp )
1384
1385                   soil_pars_f%pars_xyz(k,k2,:,:) = var_exchange_real(:,:)
1386                ENDDO
1387             ENDDO
1388             DEALLOCATE( var_dum_real_4d )
1389          ENDIF
1390       ENDIF
1391
1392       IF ( pavement_subsurface_pars_f%from_file )  THEN
1393          ALLOCATE( var_dum_real_4d(0:pavement_subsurface_pars_f%np-1,         &
1394                                    0:pavement_subsurface_pars_f%nz-1,         &
1395                                    nys:nyn,nxl:nxr) )
1396          var_dum_real_4d = pavement_subsurface_pars_f%pars_xyz
1397          DEALLOCATE( pavement_subsurface_pars_f%pars_xyz )
1398          ALLOCATE( pavement_subsurface_pars_f%pars_xyz                        &
1399                                          (0:pavement_subsurface_pars_f%np-1,  &
1400                                           0:pavement_subsurface_pars_f%nz-1,  &
1401                                           nysg:nyng,nxlg:nxrg) )
1402          pavement_subsurface_pars_f%pars_xyz = pavement_subsurface_pars_f%fill
1403
1404          DO  k2 = 0, pavement_subsurface_pars_f%nz-1
1405             DO  k = 0, pavement_subsurface_pars_f%np-1
1406                var_exchange_real(nys:nyn,nxl:nxr) =                           &
1407                                          var_dum_real_4d(k,k2,nys:nyn,nxl:nxr)
1408                CALL exchange_horiz_2d( var_exchange_real, nbgp )
1409                pavement_subsurface_pars_f%pars_xyz(k,k2,:,:) =                &
1410                                                        var_exchange_real(:,:)
1411             ENDDO
1412          ENDDO
1413          DEALLOCATE( var_dum_real_4d )
1414       ENDIF
1415
1416!
1417!--    In case of non-cyclic boundary conditions, set Neumann conditions at the
1418!--    lateral boundaries.
1419       IF ( .NOT. bc_ns_cyc )  THEN
1420          IF ( nys == 0  )  THEN
1421             IF ( albedo_type_f%from_file )                                    &
1422                albedo_type_f%var(-1,:) = albedo_type_f%var(0,:)
1423             IF ( pavement_type_f%from_file )                                  &
1424                pavement_type_f%var(-1,:) = pavement_type_f%var(0,:)
1425             IF ( soil_type_f%from_file )  THEN
1426                IF ( ALLOCATED( soil_type_f%var_2d ) )  THEN
1427                   soil_type_f%var_2d(-1,:) = soil_type_f%var_2d(0,:)
1428                ELSE
1429                   soil_type_f%var_3d(:,-1,:) = soil_type_f%var_3d(:,0,:)
1430                ENDIF
1431             ENDIF
1432             IF ( vegetation_type_f%from_file )                                &
1433                vegetation_type_f%var(-1,:) = vegetation_type_f%var(0,:)
1434             IF ( water_type_f%from_file )                                     &
1435                water_type_f%var(-1,:) = water_type_f%var(0,:)
1436             IF ( surface_fraction_f%from_file )                               &
1437                surface_fraction_f%frac(:,-1,:) = surface_fraction_f%frac(:,0,:)
1438             IF ( building_pars_f%from_file )                                  &
1439                building_pars_f%pars_xy(:,-1,:) = building_pars_f%pars_xy(:,0,:)
1440             IF ( albedo_pars_f%from_file )                                    &
1441                albedo_pars_f%pars_xy(:,-1,:) = albedo_pars_f%pars_xy(:,0,:)
1442             IF ( pavement_pars_f%from_file )                                  &
1443                pavement_pars_f%pars_xy(:,-1,:) = pavement_pars_f%pars_xy(:,0,:)
1444             IF ( vegetation_pars_f%from_file )                                &
1445                vegetation_pars_f%pars_xy(:,-1,:) =                            &
1446                                               vegetation_pars_f%pars_xy(:,0,:)
1447             IF ( water_pars_f%from_file )                                     &
1448                water_pars_f%pars_xy(:,-1,:) = water_pars_f%pars_xy(:,0,:)
1449             IF ( root_area_density_lsm_f%from_file )                          &
1450                root_area_density_lsm_f%var(:,-1,:) =                          &
1451                                            root_area_density_lsm_f%var(:,0,:)
1452             IF ( soil_pars_f%from_file )  THEN
1453                IF ( soil_pars_f%lod == 1 )  THEN
1454                   soil_pars_f%pars_xy(:,-1,:) = soil_pars_f%pars_xy(:,0,:)
1455                ELSE
1456                   soil_pars_f%pars_xyz(:,:,-1,:) = soil_pars_f%pars_xyz(:,:,0,:)
1457                ENDIF
1458             ENDIF
1459             IF ( pavement_subsurface_pars_f%from_file )                       &
1460                pavement_subsurface_pars_f%pars_xyz(:,:,-1,:) =                &
1461                                   pavement_subsurface_pars_f%pars_xyz(:,:,0,:)
1462          ENDIF
1463
1464          IF ( nyn == ny )  THEN
1465             IF ( albedo_type_f%from_file )                                    &
1466                albedo_type_f%var(ny+1,:) = albedo_type_f%var(ny,:)
1467             IF ( pavement_type_f%from_file )                                  &
1468                pavement_type_f%var(ny+1,:) = pavement_type_f%var(ny,:)
1469             IF ( soil_type_f%from_file )  THEN
1470                IF ( ALLOCATED( soil_type_f%var_2d ) )  THEN
1471                   soil_type_f%var_2d(ny+1,:) = soil_type_f%var_2d(ny,:)
1472                ELSE
1473                   soil_type_f%var_3d(:,ny+1,:) = soil_type_f%var_3d(:,ny,:)
1474                ENDIF
1475             ENDIF
1476             IF ( vegetation_type_f%from_file )                                &
1477                vegetation_type_f%var(ny+1,:) = vegetation_type_f%var(ny,:)
1478             IF ( water_type_f%from_file )                                     &
1479                water_type_f%var(ny+1,:) = water_type_f%var(ny,:)
1480             IF ( surface_fraction_f%from_file )                               &
1481                surface_fraction_f%frac(:,ny+1,:) =                            &
1482                                             surface_fraction_f%frac(:,ny,:)
1483             IF ( building_pars_f%from_file )                                  &
1484                building_pars_f%pars_xy(:,ny+1,:) =                            &
1485                                             building_pars_f%pars_xy(:,ny,:)
1486             IF ( albedo_pars_f%from_file )                                    &
1487                albedo_pars_f%pars_xy(:,ny+1,:) = albedo_pars_f%pars_xy(:,ny,:)
1488             IF ( pavement_pars_f%from_file )                                  &
1489                pavement_pars_f%pars_xy(:,ny+1,:) =                            &
1490                                             pavement_pars_f%pars_xy(:,ny,:)
1491             IF ( vegetation_pars_f%from_file )                                &
1492                vegetation_pars_f%pars_xy(:,ny+1,:) =                          &
1493                                               vegetation_pars_f%pars_xy(:,ny,:)
1494             IF ( water_pars_f%from_file )                                     &
1495                water_pars_f%pars_xy(:,ny+1,:) = water_pars_f%pars_xy(:,ny,:)
1496             IF ( root_area_density_lsm_f%from_file )                          &
1497                root_area_density_lsm_f%var(:,ny+1,:) =                        &
1498                                            root_area_density_lsm_f%var(:,ny,:)
1499             IF ( soil_pars_f%from_file )  THEN
1500                IF ( soil_pars_f%lod == 1 )  THEN
1501                   soil_pars_f%pars_xy(:,ny+1,:) = soil_pars_f%pars_xy(:,ny,:)
1502                ELSE
1503                   soil_pars_f%pars_xyz(:,:,ny+1,:) =                          &
1504                                              soil_pars_f%pars_xyz(:,:,ny,:)
1505                ENDIF
1506             ENDIF
1507             IF ( pavement_subsurface_pars_f%from_file )                       &
1508                pavement_subsurface_pars_f%pars_xyz(:,:,ny+1,:) =              &
1509                                   pavement_subsurface_pars_f%pars_xyz(:,:,ny,:)
1510          ENDIF
1511       ENDIF
1512
1513       IF ( .NOT. bc_lr_cyc )  THEN
1514          IF ( nxl == 0 )  THEN
1515            IF ( albedo_type_f%from_file )                                     &
1516                albedo_type_f%var(:,-1) = albedo_type_f%var(:,0)
1517             IF ( pavement_type_f%from_file )                                  &
1518                pavement_type_f%var(:,-1) = pavement_type_f%var(:,0)
1519             IF ( soil_type_f%from_file )  THEN
1520                IF ( ALLOCATED( soil_type_f%var_2d ) )  THEN
1521                   soil_type_f%var_2d(:,-1) = soil_type_f%var_2d(:,0)
1522                ELSE
1523                   soil_type_f%var_3d(:,:,-1) = soil_type_f%var_3d(:,:,0)
1524                ENDIF
1525             ENDIF
1526             IF ( vegetation_type_f%from_file )                                &
1527                vegetation_type_f%var(:,-1) = vegetation_type_f%var(:,0)
1528             IF ( water_type_f%from_file )                                     &
1529                water_type_f%var(:,-1) = water_type_f%var(:,0)
1530             IF ( surface_fraction_f%from_file )                               &
1531                surface_fraction_f%frac(:,:,-1) = surface_fraction_f%frac(:,:,0)
1532             IF ( building_pars_f%from_file )                                  &
1533                building_pars_f%pars_xy(:,:,-1) = building_pars_f%pars_xy(:,:,0)
1534             IF ( albedo_pars_f%from_file )                                    &
1535                albedo_pars_f%pars_xy(:,:,-1) = albedo_pars_f%pars_xy(:,:,0)
1536             IF ( pavement_pars_f%from_file )                                  &
1537                pavement_pars_f%pars_xy(:,:,-1) = pavement_pars_f%pars_xy(:,:,0)
1538             IF ( vegetation_pars_f%from_file )                                &
1539                vegetation_pars_f%pars_xy(:,:,-1) =                            &
1540                                               vegetation_pars_f%pars_xy(:,:,0)
1541             IF ( water_pars_f%from_file )                                     &
1542                water_pars_f%pars_xy(:,:,-1) = water_pars_f%pars_xy(:,:,0)
1543             IF ( root_area_density_lsm_f%from_file )                          &
1544                root_area_density_lsm_f%var(:,:,-1) =                          &
1545                                            root_area_density_lsm_f%var(:,:,0)
1546             IF ( soil_pars_f%from_file )  THEN
1547                IF ( soil_pars_f%lod == 1 )  THEN
1548                   soil_pars_f%pars_xy(:,:,-1) = soil_pars_f%pars_xy(:,:,0)
1549                ELSE
1550                   soil_pars_f%pars_xyz(:,:,:,-1) = soil_pars_f%pars_xyz(:,:,:,0)
1551                ENDIF
1552             ENDIF
1553             IF ( pavement_subsurface_pars_f%from_file )                       &
1554                pavement_subsurface_pars_f%pars_xyz(:,:,:,-1) =                &
1555                                    pavement_subsurface_pars_f%pars_xyz(:,:,:,0)
1556          ENDIF
1557
1558          IF ( nxr == nx )  THEN
1559             IF ( albedo_type_f%from_file )                                    &
1560                albedo_type_f%var(:,nx+1) = albedo_type_f%var(:,nx)
1561             IF ( pavement_type_f%from_file )                                  &
1562                pavement_type_f%var(:,nx+1) = pavement_type_f%var(:,nx)
1563             IF ( soil_type_f%from_file )  THEN
1564                IF ( ALLOCATED( soil_type_f%var_2d ) )  THEN
1565                   soil_type_f%var_2d(:,nx+1) = soil_type_f%var_2d(:,nx)
1566                ELSE
1567                   soil_type_f%var_3d(:,:,nx+1) = soil_type_f%var_3d(:,:,nx)
1568                ENDIF
1569             ENDIF
1570             IF ( vegetation_type_f%from_file )                                &
1571                vegetation_type_f%var(:,nx+1) = vegetation_type_f%var(:,nx)
1572             IF ( water_type_f%from_file )                                     &
1573                water_type_f%var(:,nx+1) = water_type_f%var(:,nx)
1574             IF ( surface_fraction_f%from_file )                               &
1575                surface_fraction_f%frac(:,:,nx+1) =                            &
1576                                             surface_fraction_f%frac(:,:,nx)
1577             IF ( building_pars_f%from_file )                                  &
1578                building_pars_f%pars_xy(:,:,nx+1) =                            &
1579                                             building_pars_f%pars_xy(:,:,nx)
1580             IF ( albedo_pars_f%from_file )                                    &
1581                albedo_pars_f%pars_xy(:,:,nx+1) = albedo_pars_f%pars_xy(:,:,nx)
1582             IF ( pavement_pars_f%from_file )                                  &
1583                pavement_pars_f%pars_xy(:,:,nx+1) =                            &
1584                                             pavement_pars_f%pars_xy(:,:,nx)
1585             IF ( vegetation_pars_f%from_file )                                &
1586                vegetation_pars_f%pars_xy(:,:,nx+1) =                          &
1587                                               vegetation_pars_f%pars_xy(:,:,nx)
1588             IF ( water_pars_f%from_file )                                     &
1589                water_pars_f%pars_xy(:,:,nx+1) = water_pars_f%pars_xy(:,:,nx)
1590             IF ( root_area_density_lsm_f%from_file )                          &
1591                root_area_density_lsm_f%var(:,:,nx+1) =                        &
1592                                            root_area_density_lsm_f%var(:,:,nx)
1593             IF ( soil_pars_f%from_file )  THEN
1594                IF ( soil_pars_f%lod == 1 )  THEN
1595                   soil_pars_f%pars_xy(:,:,nx+1) = soil_pars_f%pars_xy(:,:,nx)
1596                ELSE
1597                   soil_pars_f%pars_xyz(:,:,:,nx+1) =                          &
1598                                              soil_pars_f%pars_xyz(:,:,:,nx)
1599                ENDIF
1600             ENDIF
1601             IF ( pavement_subsurface_pars_f%from_file )                       &
1602                pavement_subsurface_pars_f%pars_xyz(:,:,:,nx+1) =              &
1603                                   pavement_subsurface_pars_f%pars_xyz(:,:,:,nx)
1604          ENDIF
1605       ENDIF
1606
1607    END SUBROUTINE netcdf_data_input_surface_data
1608
1609!------------------------------------------------------------------------------!
1610! Description:
1611! ------------
1612!> Reads orography and building information.
1613!------------------------------------------------------------------------------!
1614    SUBROUTINE netcdf_data_input_topo
1615
1616       USE control_parameters,                                                 &
1617           ONLY:  bc_lr_cyc, bc_ns_cyc, message_string, topography
1618
1619       USE indices,                                                            &
1620           ONLY:  nbgp, nx, nxl, nxr, ny, nyn, nys, nzb, nzt
1621                 
1622
1623       IMPLICIT NONE
1624
1625       CHARACTER(LEN=100), DIMENSION(:), ALLOCATABLE ::  var_names  !< variable names in static input file
1626
1627
1628       INTEGER(iwp) ::  i             !< running index along x-direction
1629       INTEGER(iwp) ::  ii            !< running index for IO blocks
1630       INTEGER(iwp) ::  id_topo       !< NetCDF id of topograhy input file
1631       INTEGER(iwp) ::  j             !< running index along y-direction
1632       INTEGER(iwp) ::  k             !< running index along z-direction
1633       INTEGER(iwp) ::  num_vars      !< number of variables in netcdf input file
1634       INTEGER(iwp) ::  skip_n_rows   !< counting variable to skip rows while reading topography file
1635
1636       INTEGER(iwp), DIMENSION(nys-nbgp:nyn+nbgp,nxl-nbgp:nxr+nbgp) ::  var_exchange_int !< dummy variables used to exchange 32-bit Integer arrays
1637
1638       REAL(wp) ::  dum           !< dummy variable to skip columns while reading topography file   
1639
1640       IF ( TRIM( topography ) /= 'read_from_file' )  RETURN
1641
1642       DO  ii = 0, io_blocks-1
1643          IF ( ii == io_group )  THEN
1644!
1645!--          Input via palm-input data standard
1646             IF ( input_pids_static )  THEN
1647#if defined ( __netcdf )
1648!
1649!--             Open file in read-only mode
1650                CALL open_read_file( TRIM( input_file_static ) //              &
1651                                     TRIM( coupling_char ), id_topo ) 
1652
1653!
1654!--             At first, inquire all variable names.
1655!--             This will be used to check whether an  input variable exist
1656!--             or not.
1657                CALL inquire_num_variables( id_topo, num_vars )
1658!
1659!--             Allocate memory to store variable names and inquire them.
1660                ALLOCATE( var_names(1:num_vars) )
1661                CALL inquire_variable_names( id_topo, var_names )
1662!
1663!--             Terrain height. First, get variable-related _FillValue attribute
1664                IF ( check_existence( var_names, 'orography_2D' ) )  THEN
1665                   terrain_height_f%from_file = .TRUE. 
1666                   CALL get_attribute( id_topo, char_fill,                     &
1667                                       terrain_height_f%fill,                  &
1668                                       .FALSE., 'orography_2D' ) 
1669!
1670!--                PE-wise reading of 2D terrain height.
1671                   ALLOCATE ( terrain_height_f%var(nys:nyn,nxl:nxr)  )
1672                   DO  i = nxl, nxr
1673                      CALL get_variable( id_topo, 'orography_2D',              &
1674                                         i, terrain_height_f%var(:,i) ) 
1675                   ENDDO
1676                ELSE
1677                   terrain_height_f%from_file = .FALSE. 
1678                ENDIF
1679
1680!
1681!--             Read building height. First, read its _FillValue attribute,
1682!--             as well as lod attribute
1683                buildings_f%from_file = .FALSE. 
1684                IF ( check_existence( var_names, 'buildings_2D' ) )  THEN
1685                   buildings_f%from_file = .TRUE. 
1686                   CALL get_attribute( id_topo, char_lod, buildings_f%lod,     &
1687                                       .FALSE., 'buildings_2D' )     
1688
1689                   CALL get_attribute( id_topo, char_fill,                     &
1690                                       buildings_f%fill1,                      &
1691                                       .FALSE., 'buildings_2D' )
1692
1693!
1694!--                Read 2D topography
1695                   IF ( buildings_f%lod == 1 )  THEN
1696                      ALLOCATE ( buildings_f%var_2d(nys:nyn,nxl:nxr) )
1697                      DO  i = nxl, nxr
1698                         CALL get_variable( id_topo, 'buildings_2D',           &
1699                                            i, buildings_f%var_2d(:,i) ) 
1700                      ENDDO
1701                   ELSE
1702                      message_string = 'NetCDF attribute lod ' //              &
1703                                       '(level of detail) is not set properly.'
1704                      CALL message( 'netcdf_data_input_mod', 'PA0999',         &
1705                                     1, 2, 0, 6, 0 )
1706                   ENDIF
1707                ENDIF
1708!
1709!--             If available, also read 3D building information. If both are
1710!--             available, use 3D information.
1711                IF ( check_existence( var_names, 'buildings_3D' ) )  THEN
1712                   buildings_f%from_file = .TRUE. 
1713                   CALL get_attribute( id_topo, char_lod, buildings_f%lod,     &
1714                                       .FALSE., 'buildings_3D' )     
1715
1716                   CALL get_attribute( id_topo, char_fill,                     &
1717                                       buildings_f%fill2,                      &
1718                                       .FALSE., 'buildings_3D' )
1719
1720                   CALL get_dimension_length( id_topo, buildings_f%nz, 'z' )
1721 
1722                   IF ( buildings_f%lod == 2 )  THEN
1723                      ALLOCATE( buildings_f%z(nzb:buildings_f%nz-1) )
1724                      CALL get_variable( id_topo, 'z', buildings_f%z )
1725
1726                      ALLOCATE( buildings_f%var_3d(nzb:buildings_f%nz-1,       &
1727                                                   nys:nyn,nxl:nxr) )
1728                      buildings_f%var_3d = 0
1729!
1730!--                   Read data PE-wise. Read yz-slices.
1731                      DO  i = nxl, nxr
1732                         DO  j = nys, nyn
1733                            CALL get_variable( id_topo, 'buildings_3D',        &
1734                                               i, j,                           &
1735                                               buildings_f%var_3d(:,j,i) )
1736                         ENDDO
1737                      ENDDO
1738                   ELSE
1739                      message_string = 'NetCDF attribute lod ' //              &
1740                                       '(level of detail) is not set properly.'
1741                      CALL message( 'netcdf_data_input_mod', 'PA0999',         &
1742                                     1, 2, 0, 6, 0 )
1743                   ENDIF
1744                ENDIF
1745!
1746!--             Read building IDs and its FillValue attribute. Further required
1747!--             for mapping buildings on top of orography.
1748                IF ( check_existence( var_names, 'building_id' ) )  THEN
1749                   building_id_f%from_file = .TRUE. 
1750                   CALL get_attribute( id_topo, char_fill,                     &
1751                                       building_id_f%fill, .FALSE.,            &
1752                                       'building_id' )
1753             
1754
1755                   ALLOCATE ( building_id_f%var(nys:nyn,nxl:nxr) )
1756                   DO  i = nxl, nxr
1757                      CALL get_variable( id_topo, 'building_id',               &
1758                                          i, building_id_f%var(:,i) ) 
1759                   ENDDO
1760                ELSE
1761                   building_id_f%from_file = .FALSE. 
1762                ENDIF
1763!
1764!--             Read building_type and required attributes.
1765                IF ( check_existence( var_names, 'building_type' ) )  THEN
1766                   building_type_f%from_file = .TRUE. 
1767                   CALL get_attribute( id_topo, char_fill,                     &
1768                                       building_type_f%fill, .FALSE.,          &
1769                                       'building_type' ) 
1770               
1771                   ALLOCATE ( building_type_f%var(nys:nyn,nxl:nxr) )
1772                   DO  i = nxl, nxr
1773                      CALL get_variable( id_topo, 'building_type',             &
1774                                         i, building_type_f%var(:,i) )
1775                   ENDDO
1776                ELSE
1777                   building_type_f%from_file = .FALSE.
1778                ENDIF
1779
1780!
1781!--             Close topography input file
1782                CALL close_input_file( id_topo )
1783#else
1784                CONTINUE
1785#endif
1786!
1787!--          ASCII input
1788             ELSE
1789
1790                OPEN( 90, FILE='TOPOGRAPHY_DATA'//TRIM( coupling_char ),       &
1791                      STATUS='OLD', FORM='FORMATTED', ERR=10 )
1792!
1793!--             Read topography PE-wise. Rows are read from nyn to nys, columns
1794!--             are read from nxl to nxr. At first, ny-nyn rows need to be skipped.
1795                skip_n_rows = 0
1796                DO WHILE ( skip_n_rows < ny - nyn )
1797                   READ( 90, * ) 
1798                   skip_n_rows = skip_n_rows + 1
1799                ENDDO
1800!
1801!--             Read data from nyn to nys and nxl to nxr. Therefore, skip
1802!--             column until nxl-1 is reached
1803                ALLOCATE ( buildings_f%var_2d(nys:nyn,nxl:nxr) )
1804                DO  j = nyn, nys, -1
1805                   READ( 90, *, ERR=11, END=11 )                               &
1806                                   ( dum, i = 0, nxl-1 ),                      &
1807                                   ( buildings_f%var_2d(j,i), i = nxl, nxr )
1808                ENDDO
1809
1810                GOTO 12
1811         
1812 10             message_string = 'file TOPOGRAPHY'//TRIM( coupling_char )//    &
1813                                 ' does not exist'
1814                CALL message( 'netcdf_data_input_mod', 'PA0208', 1, 2, 0, 6, 0 )
1815
1816 11             message_string = 'errors in file TOPOGRAPHY_DATA'//            &
1817                                 TRIM( coupling_char )
1818                CALL message( 'netcdf_data_input_mod', 'PA0209', 1, 2, 0, 6, 0 )
1819
1820 12             CLOSE( 90 )
1821                buildings_f%from_file = .TRUE.
1822
1823             ENDIF
1824
1825          ENDIF
1826#if defined( __parallel )
1827          CALL MPI_BARRIER( comm2d, ierr )
1828#endif
1829       ENDDO
1830!
1831!--    Check for minimum requirement to setup building topography. If buildings
1832!--    are provided, also an ID and a type are required.
1833!--    Note, doing this check in check_parameters
1834!--    will be too late (data will be used for grid inititialization before).
1835       IF ( input_pids_static )  THEN
1836          IF ( buildings_f%from_file  .AND.                                    &
1837               .NOT. building_id_f%from_file )  THEN                       
1838             message_string = 'If building heigths are prescribed in ' //      &
1839                              'static input file, also an ID is required.'
1840             CALL message( 'netcdf_data_input_mod', 'PA0999', 1, 2, 0, 6, 0 )
1841          ENDIF
1842       ENDIF
1843!
1844!--    In case no terrain height is provided by static input file, allocate
1845!--    array nevertheless and set terrain height to 0, which simplifies
1846!--    topography initialization.
1847       IF ( .NOT. terrain_height_f%from_file )  THEN
1848          ALLOCATE ( terrain_height_f%var(nys:nyn,nxl:nxr) )
1849          terrain_height_f%var = 0.0_wp
1850       ENDIF
1851!
1852!--    Finally, exchange 1 ghost point for building ID and type.
1853!--    In case of non-cyclic boundary conditions set Neumann conditions at the
1854!--    lateral boundaries.
1855       IF ( building_id_f%from_file )  THEN
1856          var_exchange_int                  = building_id_f%fill
1857          var_exchange_int(nys:nyn,nxl:nxr) = building_id_f%var(nys:nyn,nxl:nxr)
1858          CALL exchange_horiz_2d_int( var_exchange_int, nys, nyn, nxl, nxr, nbgp )
1859          DEALLOCATE( building_id_f%var )
1860          ALLOCATE( building_id_f%var(nys-nbgp:nyn+nbgp,nxl-nbgp:nxr+nbgp) )
1861          building_id_f%var = var_exchange_int
1862
1863          IF ( .NOT. bc_ns_cyc )  THEN
1864             IF ( nys == 0  )  building_id_f%var(-1,:)   = building_id_f%var(0,:)
1865             IF ( nyn == ny )  building_id_f%var(ny+1,:) = building_id_f%var(ny,:)
1866          ENDIF
1867          IF ( .NOT. bc_lr_cyc )  THEN
1868             IF ( nxl == 0  )  building_id_f%var(:,-1)   = building_id_f%var(:,0) 
1869             IF ( nxr == nx )  building_id_f%var(:,nx+1) = building_id_f%var(:,nx)       
1870          ENDIF
1871       ENDIF
1872
1873       IF ( building_type_f%from_file )  THEN
1874          var_exchange_int                  = INT( building_type_f%fill, KIND = 4 )
1875          var_exchange_int(nys:nyn,nxl:nxr) =                                  &
1876                          INT( building_type_f%var(nys:nyn,nxl:nxr), KIND = 4 )
1877          CALL exchange_horiz_2d_int( var_exchange_int, nys, nyn, nxl, nxr, nbgp )
1878          DEALLOCATE( building_type_f%var )
1879          ALLOCATE( building_type_f%var(nys-nbgp:nyn+nbgp,nxl-nbgp:nxr+nbgp) )
1880          building_type_f%var = INT( var_exchange_int, KIND = 1 )
1881
1882          IF ( .NOT. bc_ns_cyc )  THEN
1883             IF ( nys == 0  )  building_type_f%var(-1,:)   = building_type_f%var(0,:)
1884             IF ( nyn == ny )  building_type_f%var(ny+1,:) = building_type_f%var(ny,:)
1885          ENDIF
1886          IF ( .NOT. bc_lr_cyc )  THEN
1887             IF ( nxl == 0  )  building_type_f%var(:,-1)   = building_type_f%var(:,0) 
1888             IF ( nxr == nx )  building_type_f%var(:,nx+1) = building_type_f%var(:,nx)       
1889          ENDIF
1890       ENDIF
1891
1892    END SUBROUTINE netcdf_data_input_topo
1893
1894!------------------------------------------------------------------------------!
1895! Description:
1896! ------------
1897!> Reads initialization data of u, v, w, pt, q, geostrophic wind components,
1898!> as well as soil moisture and soil temperature, derived from larger-scale
1899!> model (COSMO) by Inifor.
1900!------------------------------------------------------------------------------!
1901    SUBROUTINE netcdf_data_input_init_3d
1902
1903       USE arrays_3d,                                                          &
1904           ONLY:  q, pt, u, v, w
1905
1906       USE control_parameters,                                                 &
1907           ONLY:  bc_lr_cyc, bc_ns_cyc, forcing, humidity, land_surface,       &
1908                  message_string, neutral, surface_pressure
1909
1910       USE indices,                                                            &
1911           ONLY:  nx, nxl, nxlu, nxr, ny, nyn, nys, nysv, nzb, nzt
1912
1913       IMPLICIT NONE
1914
1915       CHARACTER(LEN=100), DIMENSION(:), ALLOCATABLE ::  var_names
1916
1917       INTEGER(iwp) ::  i          !< running index along x-direction
1918       INTEGER(iwp) ::  ii         !< running index for IO blocks
1919       INTEGER(iwp) ::  id_dynamic !< NetCDF id of dynamic input file
1920       INTEGER(iwp) ::  j          !< running index along y-direction
1921       INTEGER(iwp) ::  k          !< running index along z-direction
1922       INTEGER(iwp) ::  num_vars   !< number of variables in netcdf input file
1923       INTEGER(iwp) ::  off_i      !< offset in x-direction used for reading the u-component
1924       INTEGER(iwp) ::  off_j      !< offset in y-direction used for reading the v-component
1925
1926!
1927!--    Skip routine if no input file with dynamic input data is available.
1928       IF ( .NOT. input_pids_dynamic )  RETURN
1929!
1930!--    Please note, Inifor is designed to provide initial data for u and v for
1931!--    the prognostic grid points in case of lateral Dirichlet conditions.
1932!--    This means that Inifor provides data from nxlu:nxr (for u) and
1933!--    from nysv:nyn (for v) at the left and south domain boundary, respectively.
1934!--    However, as work-around for the moment, PALM will run with cyclic
1935!--    conditions and will be initialized with data provided by Inifor
1936!--    boundaries in case of Dirichlet.
1937!--    Hence, simply set set nxlu/nysv to 1 (will be reset to its original value
1938!--    at the end of this routine.
1939       IF ( bc_lr_cyc  .AND.  nxl == 0 )  nxlu = 1 
1940       IF ( bc_ns_cyc  .AND.  nys == 0 )  nysv = 1
1941
1942       DO  ii = 0, io_blocks-1
1943          IF ( ii == io_group )  THEN
1944#if defined ( __netcdf )
1945!
1946!--          Open file in read-only mode
1947             CALL open_read_file( TRIM( input_file_dynamic ) //                &
1948                                  TRIM( coupling_char ), id_dynamic ) 
1949
1950!
1951!--          At first, inquire all variable names.
1952             CALL inquire_num_variables( id_dynamic, num_vars )
1953!
1954!--          Allocate memory to store variable names.
1955             ALLOCATE( var_names(1:num_vars) )
1956             CALL inquire_variable_names( id_dynamic, var_names )
1957!
1958!--          Read vertical dimension of scalar und w grid. Will be used for
1959!--          inter- and extrapolation in case of stretched numeric grid.
1960!--          This will be removed when Inifor is able to handle stretched grids.
1961             CALL get_dimension_length( id_dynamic, init_3d%nzu, 'z'     )
1962             CALL get_dimension_length( id_dynamic, init_3d%nzw, 'zw'    )
1963             CALL get_dimension_length( id_dynamic, init_3d%nzs, 'depth' )
1964!
1965!--          Read also the horizontal dimensions. These are used just used fo
1966!--          checking the compatibility with the PALM grid before reading.
1967             CALL get_dimension_length( id_dynamic, init_3d%nx,  'x'  )
1968             CALL get_dimension_length( id_dynamic, init_3d%nxu, 'xu' )
1969             CALL get_dimension_length( id_dynamic, init_3d%ny,  'y'  )
1970             CALL get_dimension_length( id_dynamic, init_3d%nyv, 'yv' )
1971!
1972!--          Check for correct horizontal dimension. Please note, u- and v-grid
1973!--          hase 1 grid point less on Inifor grid.
1974             IF ( init_3d%nx-1 /= nx  .OR.  init_3d%nxu-1 /= nx - 1  .OR.      &
1975                  init_3d%ny-1 /= ny  .OR.  init_3d%nyv-1 /= ny - 1 )  THEN
1976                message_string = 'Number of inifor grid points does not '  //  &
1977                                 'match the number of numeric grid points.'
1978                CALL message( 'netcdf_data_input_mod', 'PA0999', 1, 2, 0, 6, 0 )
1979             ENDIF
1980!
1981!--          Read vertical dimensions. Later, these are required for eventual
1982!--          inter- and extrapolations of the initialization data.
1983             IF ( check_existence( var_names, 'z' ) )  THEN
1984                ALLOCATE( init_3d%zu_atmos(1:init_3d%nzu) )
1985                CALL get_variable( id_dynamic, 'z', init_3d%zu_atmos )
1986             ENDIF
1987             IF ( check_existence( var_names, 'zw' ) )  THEN
1988                ALLOCATE( init_3d%zw_atmos(1:init_3d%nzw) )
1989                CALL get_variable( id_dynamic, 'zw', init_3d%zw_atmos )
1990             ENDIF
1991             IF ( check_existence( var_names, 'depth' ) )  THEN
1992                ALLOCATE( init_3d%z_soil(1:init_3d%nzs) )
1993                CALL get_variable( id_dynamic, 'depth', init_3d%z_soil )
1994             ENDIF
1995!
1996!--          Read geostrophic wind components
1997!              IF ( check_existence( var_names, 'ls_forcing_ug' ) )  THEN
1998!
1999!              ENDIF
2000!              IF ( check_existence( var_names, 'ls_forcing_vg' ) )  THEN
2001!
2002!              ENDIF
2003
2004!
2005!--          Read inital 3D data of u, v, w, pt and q,
2006!--          derived from COSMO model. Read PE-wise yz-slices.
2007!--          Please note, the u-, v- and w-component are defined on different
2008!--          grids with one element less in the x-, y-,
2009!--          and z-direction, respectively. Hence, reading is subdivided
2010!--          into separate loops. Moreover, i and j are used
2011!--          as start index in the NF90 interface.
2012!--          The passed arguments for u, and v are (i,j)-1, respectively,
2013!--          in contrast to the remaining quantities. This is because in case
2014!--          of forcing is applied, the input data for u and v has one
2015!--          element less along the x- and y-direction respectively.
2016!--          Read u-component
2017             IF ( check_existence( var_names, 'init_u' ) )  THEN
2018!
2019!--             Read attributes for the fill value and level-of-detail
2020                CALL get_attribute( id_dynamic, char_fill, init_3d%fill_u,     &
2021                                    .FALSE., 'init_u' )
2022                CALL get_attribute( id_dynamic, char_lod, init_3d%lod_u,       &
2023                                    .FALSE., 'init_u' )
2024!
2025!--             level-of-detail 1 - read initialization profile
2026                IF ( init_3d%lod_u == 1 )  THEN
2027                   ALLOCATE( init_3d%u_init(nzb:nzt+1) )
2028                   init_3d%u_init = 0.0_wp
2029
2030                   CALL get_variable( id_dynamic, 'init_u',                    &
2031                                      init_3d%u_init(nzb+1:nzt+1) )
2032!
2033!--             level-of-detail 2 - read 3D initialization data
2034                ELSEIF ( init_3d%lod_u == 2 )  THEN
2035!
2036!--                Set offset value. In case of Dirichlet conditions at the left
2037!--                domain boundary, the u component starts at nxl+1. This case,
2038!--                the passed start-index for reading the NetCDF data is shifted
2039!--                by -1. 
2040                   off_i = 1 !MERGE( 1, 0, forcing )
2041
2042                   DO  i = nxlu, nxr
2043                      DO  j = nys, nyn   
2044                         CALL get_variable( id_dynamic, 'init_u', i-off_i, j,  &
2045                                            u(nzb+1:nzt+1,j,i) )
2046                      ENDDO
2047                   ENDDO
2048
2049                ENDIF
2050                init_3d%from_file_u = .TRUE.
2051             ENDIF
2052!
2053!--          Read v-component
2054             IF ( check_existence( var_names, 'init_v' ) )  THEN
2055!
2056!--             Read attributes for the fill value and level-of-detail
2057                CALL get_attribute( id_dynamic, char_fill, init_3d%fill_v,     &
2058                                    .FALSE., 'init_v' )
2059                CALL get_attribute( id_dynamic, char_lod, init_3d%lod_v,       &
2060                                    .FALSE., 'init_v' )
2061!
2062!--             level-of-detail 1 - read initialization profile
2063                IF ( init_3d%lod_v == 1 )  THEN
2064                   ALLOCATE( init_3d%v_init(nzb:nzt+1) )
2065                   init_3d%v_init = 0.0_wp
2066
2067                   CALL get_variable( id_dynamic, 'init_v',                    &
2068                                      init_3d%v_init(nzb+1:nzt+1) )
2069
2070!
2071!--             level-of-detail 2 - read 3D initialization data
2072                ELSEIF ( init_3d%lod_v == 2 )  THEN
2073!
2074!--                Set offset value. In case of Dirichlet conditions at the south
2075!--                domain boundary, the v component starts at nys+1. This case,
2076!--                the passed start-index for reading the NetCDF data is shifted
2077!--                by -1. 
2078                   off_j = 1 !MERGE( 1, 0, forcing )
2079
2080                   DO  i = nxl, nxr
2081                      DO  j = nysv, nyn   
2082                         CALL get_variable( id_dynamic, 'init_v', i, j-off_j,  &
2083                                            v(nzb+1:nzt+1,j,i) )
2084                      ENDDO
2085                   ENDDO
2086
2087                ENDIF
2088                init_3d%from_file_v = .TRUE.
2089             ENDIF
2090!
2091!--          Read w-component
2092             IF ( check_existence( var_names, 'init_w' ) )  THEN
2093!
2094!--             Read attributes for the fill value and level-of-detail
2095                CALL get_attribute( id_dynamic, char_fill, init_3d%fill_w,     &
2096                                    .FALSE., 'init_w' )
2097                CALL get_attribute( id_dynamic, char_lod, init_3d%lod_w,       &
2098                                    .FALSE., 'init_w' )
2099!
2100!--             level-of-detail 1 - read initialization profile
2101                IF ( init_3d%lod_w == 1 )  THEN
2102                   ALLOCATE( init_3d%w_init(nzb:nzt+1) )
2103                   init_3d%w_init = 0.0_wp
2104
2105                   CALL get_variable( id_dynamic, 'init_w',                    &
2106                                      init_3d%w_init(nzb+1:nzt) )
2107
2108!
2109!--             level-of-detail 2 - read 3D initialization data
2110                ELSEIF ( init_3d%lod_w == 2 )  THEN
2111                   DO  i = nxl, nxr
2112                      DO  j = nys, nyn
2113                         CALL get_variable( id_dynamic, 'init_w', i, j,        &
2114                                            w(nzb+1:nzt,j,i) )
2115                      ENDDO
2116                   ENDDO
2117
2118                ENDIF
2119                init_3d%from_file_w = .TRUE.
2120             ENDIF
2121!
2122!--          Read potential temperature
2123             IF ( .NOT. neutral )  THEN     
2124                IF ( check_existence( var_names, 'init_pt' ) )  THEN   
2125!
2126!--                Read attributes for the fill value and level-of-detail
2127                   CALL get_attribute( id_dynamic, char_fill, init_3d%fill_pt, &
2128                                       .FALSE., 'init_pt' )
2129                   CALL get_attribute( id_dynamic, char_lod, init_3d%lod_pt,   &
2130                                       .FALSE., 'init_pt' )
2131!
2132!--                level-of-detail 1 - read initialization profile
2133                   IF ( init_3d%lod_pt == 1 )  THEN
2134                      ALLOCATE( init_3d%pt_init(nzb:nzt+1) )
2135
2136                      CALL get_variable( id_dynamic, 'init_pt',                &
2137                                         init_3d%pt_init(nzb+1:nzt+1) )
2138!
2139!--                   Set Neumann surface boundary condition for initial profil
2140                      init_3d%pt_init(nzb) = init_3d%pt_init(nzb+1)
2141!
2142!--                level-of-detail 2 - read 3D initialization data
2143                   ELSEIF ( init_3d%lod_pt == 2 )  THEN
2144                      DO  i = nxl, nxr
2145                         DO  j = nys, nyn
2146                            CALL get_variable( id_dynamic, 'init_pt', i, j,    &
2147                                               pt(nzb+1:nzt+1,j,i) )
2148                         ENDDO
2149                      ENDDO
2150
2151                   ENDIF
2152                   init_3d%from_file_pt = .TRUE.
2153                ENDIF
2154             ENDIF
2155!
2156!--          Read mixing ratio
2157             IF ( humidity )  THEN
2158                IF ( check_existence( var_names, 'init_qv' ) )  THEN
2159!
2160!--                Read attributes for the fill value and level-of-detail
2161                   CALL get_attribute( id_dynamic, char_fill, init_3d%fill_q,  &
2162                                       .FALSE., 'init_qv' )
2163                   CALL get_attribute( id_dynamic, char_lod, init_3d%lod_q,    &
2164                                       .FALSE., 'init_qv' )
2165!
2166!--                level-of-detail 1 - read initialization profile
2167                   IF ( init_3d%lod_q == 1 )  THEN
2168                      ALLOCATE( init_3d%q_init(nzb:nzt+1) )
2169
2170                      CALL get_variable( id_dynamic, 'init_qv',               &
2171                                         init_3d%q_init(nzb+1:nzt+1) )
2172!
2173!--                   Set Neumann surface boundary condition for initial profil
2174                      init_3d%q_init(nzb) = init_3d%q_init(nzb+1)
2175
2176!
2177!--                level-of-detail 2 - read 3D initialization data
2178                   ELSEIF ( init_3d%lod_q == 2 )  THEN
2179                      DO  i = nxl, nxr
2180                         DO  j = nys, nyn
2181                            CALL get_variable( id_dynamic, 'init_qv', i, j,    &
2182                                               q(nzb+1:nzt+1,j,i) )
2183                         ENDDO
2184                      ENDDO
2185
2186                   ENDIF
2187                   init_3d%from_file_q = .TRUE.
2188                ENDIF
2189             ENDIF
2190!
2191!--          Read soil moisture
2192             IF ( land_surface )  THEN
2193                IF ( check_existence( var_names, 'init_soil_m' ) )  THEN
2194!
2195!--                Read attributes for the fill value and level-of-detail
2196                   CALL get_attribute( id_dynamic, char_fill,                  &
2197                                       init_3d%fill_msoil,                     &
2198                                       .FALSE., 'init_soil_m' )
2199                   CALL get_attribute( id_dynamic, char_lod,                   &
2200                                       init_3d%lod_msoil,                      &
2201                                       .FALSE., 'init_soil_m' )
2202!
2203!--                level-of-detail 1 - read initialization profile
2204                   IF ( init_3d%lod_msoil == 1 )  THEN
2205                      ALLOCATE( init_3d%msoil_init(0:init_3d%nzs-1) )
2206
2207                      CALL get_variable( id_dynamic, 'init_soil_m',            &
2208                                         init_3d%msoil_init(0:init_3d%nzs-1) )
2209!
2210!--                level-of-detail 2 - read 3D initialization data
2211                   ELSEIF ( init_3d%lod_msoil == 2 )  THEN  ! need to be corrected
2212                      ALLOCATE ( init_3d%msoil(0:init_3d%nzs-1,nys:nyn,nxl:nxr) )
2213                      DO  i = nxl, nxr
2214                         DO  j = nys, nyn
2215                            CALL get_variable( id_dynamic, 'init_soil_m', i, j,&
2216                                               init_3d%msoil(0:init_3d%nzs-1,j,i) )
2217                         ENDDO
2218                      ENDDO
2219                   ENDIF
2220                   init_3d%from_file_msoil = .TRUE.
2221                ENDIF
2222!
2223!--             Read soil temperature
2224                IF ( check_existence( var_names, 'init_soil_t' ) )  THEN
2225!
2226!--                Read attributes for the fill value and level-of-detail
2227                   CALL get_attribute( id_dynamic, char_fill,                  &
2228                                       init_3d%fill_tsoil,                     &
2229                                       .FALSE., 'init_soil_t' )
2230                   CALL get_attribute( id_dynamic, char_lod,                   &
2231                                       init_3d%lod_tsoil,                      &
2232                                       .FALSE., 'init_soil_t' )
2233!
2234!--                level-of-detail 1 - read initialization profile
2235                   IF ( init_3d%lod_tsoil == 1 )  THEN
2236                      ALLOCATE( init_3d%tsoil_init(0:init_3d%nzs-1) )
2237
2238                      CALL get_variable( id_dynamic, 'init_soil_t',            &
2239                                         init_3d%tsoil_init(0:init_3d%nzs-1) )
2240
2241!
2242!--                level-of-detail 2 - read 3D initialization data
2243                   ELSEIF ( init_3d%lod_tsoil == 2 )  THEN  ! need to be corrected
2244                      ALLOCATE ( init_3d%tsoil(0:init_3d%nzs-1,nys:nyn,nxl:nxr) )
2245                      DO  i = nxl, nxr
2246                         DO  j = nys, nyn
2247                            CALL get_variable( id_dynamic, 'init_soil_t', i, j,&
2248                                               init_3d%tsoil(0:init_3d%nzs-1,j,i) )
2249                         ENDDO
2250                      ENDDO
2251                   ENDIF
2252                   init_3d%from_file_tsoil = .TRUE.
2253                ENDIF
2254             ENDIF
2255!
2256!--          Close input file
2257             CALL close_input_file( id_dynamic )
2258#endif
2259          ENDIF
2260#if defined( __parallel )
2261          CALL MPI_BARRIER( comm2d, ierr )
2262#endif
2263       ENDDO
2264!
2265!--    Finally, check if the input data has any fill values.
2266       IF ( init_3d%from_file_u )  THEN
2267          IF ( ANY( u(nzb+1:nzt+1,nys:nyn,nxlu:nxr) == init_3d%fill_u ) )  THEN
2268             message_string = 'NetCDF input for u_init must not contain ' //   &
2269                              'any _FillValues'
2270             CALL message( 'netcdf_data_input_mod', 'PA0999', 2, 2, 0, 6, 0 )
2271          ENDIF
2272       ENDIF
2273       IF ( init_3d%from_file_v )  THEN
2274          IF ( ANY( v(nzb+1:nzt+1,nysv:nyn,nxl:nxr) == init_3d%fill_v ) )  THEN
2275             message_string = 'NetCDF input for v_init must not contain ' //   &
2276                              'any _FillValues'
2277             CALL message( 'netcdf_data_input_mod', 'PA0999', 2, 2, 0, 6, 0 )
2278          ENDIF
2279       ENDIF
2280       IF ( init_3d%from_file_w )  THEN
2281          IF ( ANY( w(nzb+1:nzt,nys:nyn,nxl:nxr) == init_3d%fill_w ) )  THEN
2282             message_string = 'NetCDF input for w_init must not contain ' //   &
2283                              'any _FillValues'
2284             CALL message( 'netcdf_data_input_mod', 'PA0999', 2, 2, 0, 6, 0 )
2285          ENDIF
2286       ENDIF
2287       IF ( init_3d%from_file_pt )  THEN
2288          IF ( ANY( pt(nzb+1:nzt+1,nys:nyn,nxl:nxr) == init_3d%fill_pt ) )  THEN
2289             message_string = 'NetCDF input for pt_init must not contain ' //  &
2290                              'any _FillValues'
2291             CALL message( 'netcdf_data_input_mod', 'PA0999', 2, 2, 0, 6, 0 )
2292          ENDIF
2293       ENDIF
2294       IF ( init_3d%from_file_q )  THEN
2295          IF ( ANY( q(nzb+1:nzt+1,nys:nyn,nxl:nxr) == init_3d%fill_q ) )  THEN
2296             message_string = 'NetCDF input for q_init must not contain ' //   &
2297                              'any _FillValues'
2298             CALL message( 'netcdf_data_input_mod', 'PA0999', 2, 2, 0, 6, 0 )
2299          ENDIF
2300       ENDIF
2301!
2302!--    Workaround for cyclic conditions. Please see above for further explanation.
2303       IF ( bc_lr_cyc  .AND.  nxl == 0 )  nxlu = nxl 
2304       IF ( bc_ns_cyc  .AND.  nys == 0 )  nysv = nys
2305
2306    END SUBROUTINE netcdf_data_input_init_3d
2307
2308!------------------------------------------------------------------------------!
2309! Description:
2310! ------------
2311!> Reads data at lateral and top boundaries derived from larger-scale model
2312!> (COSMO) by Inifor.
2313!------------------------------------------------------------------------------!
2314    SUBROUTINE netcdf_data_input_lsf
2315
2316       USE control_parameters,                                                 &
2317           ONLY:  force_bound_l, force_bound_n, force_bound_r, force_bound_s,  &
2318                  forcing, humidity, message_string, neutral, simulated_time
2319
2320       USE indices,                                                            &
2321           ONLY:  nx, nxl, nxlu, nxr, ny, nyn, nys, nysv, nzb, nzt
2322
2323       IMPLICIT NONE
2324
2325
2326       INTEGER(iwp) ::  i          !< running index along x-direction
2327       INTEGER(iwp) ::  ii         !< running index for IO blocks
2328       INTEGER(iwp) ::  id_dynamic !< NetCDF id of dynamic input file
2329       INTEGER(iwp) ::  j          !< running index along y-direction
2330       INTEGER(iwp) ::  k          !< running index along z-direction
2331       INTEGER(iwp) ::  num_vars   !< number of variables in netcdf input file
2332       INTEGER(iwp) ::  t          !< running index time dimension
2333
2334       REAL(wp) ::  dum           !< dummy variable to skip columns while reading topography file   
2335
2336       force%from_file = MERGE( .TRUE., .FALSE., input_pids_dynamic ) 
2337!
2338!--    Skip input if no forcing from larger-scale models is applied.
2339       IF ( .NOT. forcing )  RETURN
2340
2341       DO  ii = 0, io_blocks-1
2342          IF ( ii == io_group )  THEN
2343#if defined ( __netcdf )
2344!
2345!--          Open file in read-only mode
2346             CALL open_read_file( TRIM( input_file_dynamic ) //                &
2347                                  TRIM( coupling_char ), id_dynamic ) 
2348!
2349!--          Initialize INIFOR forcing.
2350             IF ( .NOT. force%init )  THEN
2351!
2352!--             At first, inquire all variable names.
2353                CALL inquire_num_variables( id_dynamic, num_vars )
2354!
2355!--             Allocate memory to store variable names.
2356                ALLOCATE( force%var_names(1:num_vars) )
2357                CALL inquire_variable_names( id_dynamic, force%var_names )
2358!
2359!--             Read time dimension, allocate memory and finally read time array
2360                CALL get_dimension_length( id_dynamic, force%nt, 'time' )
2361
2362                IF ( check_existence( force%var_names, 'time' ) )  THEN
2363                   ALLOCATE( force%time(0:force%nt-1) )
2364                   CALL get_variable( id_dynamic, 'time', force%time )
2365                ENDIF
2366!
2367!--             Read vertical dimension of scalar und w grid
2368                CALL get_dimension_length( id_dynamic, force%nzu, 'z' )
2369                CALL get_dimension_length( id_dynamic, force%nzw, 'zw' )
2370
2371                IF ( check_existence( force%var_names, 'z' ) )  THEN
2372                   ALLOCATE( force%zu_atmos(1:force%nzu) )
2373                   CALL get_variable( id_dynamic, 'z', force%zu_atmos )
2374                ENDIF
2375                IF ( check_existence( force%var_names, 'zw' ) )  THEN
2376                   ALLOCATE( force%zw_atmos(1:force%nzw) )
2377                   CALL get_variable( id_dynamic, 'zw', force%zw_atmos )
2378                ENDIF
2379
2380!
2381!--             Read surface pressure
2382                IF ( check_existence( force%var_names,                         &
2383                                  'surface_forcing_surface_pressure' ) )  THEN
2384                   ALLOCATE( force%surface_pressure(0:force%nt-1) )
2385                   CALL get_variable( id_dynamic,                              &
2386                                      'surface_forcing_surface_pressure',      &
2387                                      force%surface_pressure )
2388                ENDIF
2389!
2390!--             Set control flag to indicate that initialization is already done
2391                force%init = .TRUE.
2392
2393             ENDIF
2394
2395!
2396!--          Obtain time index for current input starting at 0.
2397!--          @todo: At the moment time, in INIFOR and simulated time correspond
2398!--                 to each other. If required, adjust to daytime.
2399             force%tind = MINLOC( ABS( force%time - simulated_time ), DIM = 1 )&
2400                          - 1
2401             force%tind_p = force%tind + 1 
2402!
2403!--          Read data at lateral and top boundaries. Please note, at left and
2404!--          right domain boundary, yz-layers are read for u, v, w, pt and q.
2405!--          For the v-component, the data starts at nysv, while for the other
2406!--          quantities the data starts at nys. This is equivalent at the north
2407!--          and south domain boundary for the u-component.
2408!--          The function get_variable_bc assumes the start indices with respect
2409!--          to the netcdf file convention (data starts at index 1). For this
2410!--          reason, nys+1 / nxl+1 are passed instead of nys / nxl. For the
2411!--          the u- and v-component at the north/south, and left/right boundary,
2412!--          nxlu and nysv are passed, respectively, since these always starts
2413!--          at index 1 in case of forcing.
2414
2415             IF ( force_bound_l )  THEN
2416                DO  j = nys, nyn
2417                   DO  t = force%tind, force%tind_p
2418                      CALL get_variable_bc( id_dynamic, 'ls_forcing_left_u',   &
2419                                      t+1,                                     &
2420                                      nzb+1, nzt+1-(nzb+1)+1,                  &
2421                                      j+1, 1,                                  &
2422                                      force%u_left(t-force%tind,nzb+1:nzt+1,j) )
2423                   ENDDO
2424                ENDDO
2425
2426                DO  j = nysv, nyn
2427                   DO  t = force%tind, force%tind_p
2428                      CALL get_variable_bc( id_dynamic, 'ls_forcing_left_v',   &
2429                                      t+1,                                     &
2430                                      nzb+1, nzt+1-(nzb+1)+1,                  &
2431                                      j, 1,                                    &
2432                                      force%v_left(t-force%tind,nzb+1:nzt+1,j) )
2433                   ENDDO
2434                ENDDO
2435                DO  j = nys, nyn
2436                   DO  t = force%tind, force%tind_p
2437                      CALL get_variable_bc( id_dynamic,                        &
2438                                      'ls_forcing_left_w',                     &
2439                                      t+1,                                     &
2440                                      nzb+1, nzt-(nzb+1) + 1,                  &
2441                                      j+1, 1,                                  &
2442                                      force%w_left(t-force%tind,nzb+1:nzt,j) )
2443                   ENDDO
2444                ENDDO
2445                IF ( .NOT. neutral )  THEN
2446                   DO  j = nys, nyn
2447                      DO  t = force%tind, force%tind_p
2448                         CALL get_variable_bc( id_dynamic,                     &
2449                                      'ls_forcing_left_pt',                    &
2450                                      t+1,                                     &
2451                                      nzb+1, nzt+1-(nzb+1)+1,                  &
2452                                      j+1, 1,                                  &
2453                                      force%pt_left(t-force%tind,nzb+1:nzt+1,j) )
2454                      ENDDO
2455                   ENDDO
2456                ENDIF
2457                IF ( humidity )  THEN
2458                   DO  j = nys, nyn
2459                      DO  t = force%tind, force%tind_p
2460                         CALL get_variable_bc( id_dynamic,                     &
2461                                      'ls_forcing_left_qv',                    &
2462                                      t+1,                                     &
2463                                      nzb+1, nzt+1-(nzb+1)+1,                  &
2464                                      j+1, 1,                                  &
2465                                      force%q_left(t-force%tind,nzb+1:nzt+1,j) )
2466                      ENDDO
2467                   ENDDO
2468                ENDIF
2469             ENDIF
2470
2471             IF ( force_bound_r )  THEN
2472                DO  j = nys, nyn
2473                   DO  t = force%tind, force%tind_p
2474                      CALL get_variable_bc( id_dynamic, 'ls_forcing_right_u',  &
2475                                      t+1,                                     &
2476                                      nzb+1, nzt+1-(nzb+1)+1,                  &
2477                                      j+1, 1,                                  &
2478                                      force%u_right(t-force%tind,nzb+1:nzt+1,j) )
2479                   ENDDO
2480                ENDDO
2481                DO  j = nysv, nyn
2482                   DO  t = force%tind, force%tind_p
2483                      CALL get_variable_bc( id_dynamic, 'ls_forcing_right_v',  &
2484                                      t+1,                                     &
2485                                      nzb+1, nzt+1-(nzb+1)+1,                  &
2486                                      j, 1,                                    &
2487                                      force%v_right(t-force%tind,nzb+1:nzt+1,j) )
2488                   ENDDO
2489                ENDDO
2490                DO  j = nys, nyn
2491                   DO  t = force%tind, force%tind_p
2492                      CALL get_variable_bc( id_dynamic, 'ls_forcing_right_w',  &
2493                                      t+1,                                     &
2494                                      nzb+1, nzt-(nzb+1)+1,                    &
2495                                      j+1, 1,                                  &
2496                                      force%w_right(t-force%tind,nzb+1:nzt,j) )
2497                   ENDDO
2498                ENDDO
2499                IF ( .NOT. neutral )  THEN
2500                   DO  j = nys, nyn
2501                      DO  t = force%tind, force%tind_p
2502                         CALL get_variable_bc( id_dynamic,                     &
2503                                      'ls_forcing_right_pt',                   &
2504                                      t+1,                                     &
2505                                      nzb+1, nzt+1-(nzb+1)+1,                  &
2506                                      j+1, 1,                                  &
2507                                      force%pt_right(t-force%tind,nzb+1:nzt+1,j) )
2508                      ENDDO
2509                   ENDDO
2510                ENDIF
2511                IF ( humidity )  THEN
2512                   DO  j = nys, nyn
2513                      DO  t = force%tind, force%tind_p
2514                         CALL get_variable_bc( id_dynamic,                     &
2515                                      'ls_forcing_right_qv',                   &
2516                                      t+1,                                     &
2517                                      nzb+1, nzt+1-(nzb+1)+1,                  &
2518                                      j+1, 1,                                  &
2519                                      force%q_right(t-force%tind,nzb+1:nzt+1,j) )
2520                      ENDDO
2521                   ENDDO
2522                ENDIF
2523             ENDIF
2524
2525             IF ( force_bound_n )  THEN
2526                DO  i = nxlu, nxr
2527                   DO  t = force%tind, force%tind_p
2528                      CALL get_variable_bc( id_dynamic, 'ls_forcing_north_u',  &
2529                                      t+1,                                     &
2530                                      nzb+1, nzt+1-(nzb+1)+1,                  &
2531                                      i, 1,                                    &
2532                                      force%u_north(t-force%tind,nzb+1:nzt+1,i) )
2533                   ENDDO
2534                ENDDO
2535                DO  i = nxl, nxr
2536                   DO  t = force%tind, force%tind_p
2537                      CALL get_variable_bc( id_dynamic, 'ls_forcing_north_v',  &
2538                                      t+1,                                     &
2539                                      nzb+1, nzt+1-(nzb+1)+1,                  &
2540                                      i+1, 1,                                  &
2541                                      force%v_north(t-force%tind,nzb+1:nzt+1,i) )
2542                   ENDDO
2543                ENDDO
2544                DO  i = nxl, nxr
2545                   DO  t = force%tind, force%tind_p
2546                      CALL get_variable_bc( id_dynamic, 'ls_forcing_north_w',  &
2547                                      t+1,                                     &
2548                                      nzb+1, nzt-(nzb+1)+1,                    &
2549                                      i+1, 1,                                  &
2550                                      force%w_north(t-force%tind,nzb+1:nzt,i) )
2551                   ENDDO
2552                ENDDO
2553                IF ( .NOT. neutral )  THEN
2554                   DO  i = nxl, nxr
2555                      DO  t = force%tind, force%tind_p
2556                         CALL get_variable_bc( id_dynamic,                     &
2557                                      'ls_forcing_north_pt',                   &
2558                                      t+1,                                     &
2559                                      nzb+1, nzt+1-(nzb+1)+1,                  &
2560                                      i+1, 1,                                  &
2561                                      force%pt_north(t-force%tind,nzb+1:nzt+1,i) )
2562                      ENDDO
2563                   ENDDO
2564                ENDIF
2565                IF ( humidity )  THEN
2566                   DO  i = nxl, nxr
2567                      DO  t = force%tind, force%tind_p
2568                         CALL get_variable_bc( id_dynamic,                     &
2569                                      'ls_forcing_north_qv',                   &
2570                                      t+1,                                     &
2571                                      nzb+1, nzt+1-(nzb+1)+1,                  &
2572                                      i+1, 1,                                  &
2573                                      force%q_north(t-force%tind,nzb+1:nzt+1,i) )
2574                      ENDDO
2575                   ENDDO
2576                ENDIF
2577             ENDIF
2578
2579             IF ( force_bound_s )  THEN
2580                DO  i = nxlu, nxr
2581                   DO  t = force%tind, force%tind_p
2582                      CALL get_variable_bc( id_dynamic, 'ls_forcing_south_u',  &
2583                                      t+1,                                     &
2584                                      nzb+1, nzt+1-(nzb+1)+1,                  &
2585                                      i, 1,                                    &
2586                                      force%u_south(t-force%tind,nzb+1:nzt+1,i) )
2587                   ENDDO
2588                ENDDO
2589                DO  i = nxl, nxr
2590                   DO  t = force%tind, force%tind_p
2591                      CALL get_variable_bc( id_dynamic, 'ls_forcing_south_v',  &
2592                                      t+1,                                     &
2593                                      nzb+1, nzt+1-(nzb+1)+1,                  &
2594                                      i+1, 1,                                  &
2595                                      force%v_south(t-force%tind,nzb+1:nzt+1,i) )
2596                   ENDDO
2597                ENDDO
2598                DO  i = nxl, nxr
2599                   DO  t = force%tind, force%tind_p
2600                      CALL get_variable_bc( id_dynamic, 'ls_forcing_south_w',  &
2601                                      t+1,                                     &
2602                                      nzb+1, nzt-(nzb+1)+1,                    &
2603                                      i+1, 1,                                  &
2604                                      force%w_south(t-force%tind,nzb+1:nzt,i) )
2605                   ENDDO
2606                ENDDO
2607                IF ( .NOT. neutral )  THEN
2608                   DO  i = nxl, nxr
2609                      DO  t = force%tind, force%tind_p
2610                         CALL get_variable_bc( id_dynamic,                     &
2611                                      'ls_forcing_south_pt',                   &
2612                                      t+1,                                     &
2613                                      nzb+1, nzt+1-(nzb+1)+1,                  &
2614                                      i+1, 1,                                  &
2615                                      force%pt_south(t-force%tind,nzb+1:nzt+1,i) )
2616                      ENDDO
2617                   ENDDO
2618                ENDIF
2619                IF ( humidity )  THEN
2620                   DO  i = nxl, nxr
2621                      DO  t = force%tind, force%tind_p
2622                         CALL get_variable_bc( id_dynamic,                     &
2623                                      'ls_forcing_south_qv',                   &
2624                                      t+1,                                     &
2625                                      nzb+1, nzt+1-(nzb+1)+1,                  &
2626                                      i+1, 1,                                  &
2627                                      force%q_south(t-force%tind,nzb+1:nzt+1,i) )
2628                      ENDDO
2629                   ENDDO
2630                ENDIF
2631             ENDIF
2632!
2633!--          Top boundary
2634             DO  i = nxlu, nxr
2635                DO  t = force%tind, force%tind_p
2636                   CALL get_variable_bc( id_dynamic, 'ls_forcing_top_u',       &
2637                                   t+1,                                        &
2638                                   nys+1, nyn-nys+1,                           &
2639                                   i, 1,                                       &
2640                                   force%u_top(t-force%tind,nys:nyn,i) )
2641                ENDDO
2642             ENDDO
2643             DO  i = nxl, nxr
2644                DO  t = force%tind, force%tind_p
2645                   CALL get_variable_bc( id_dynamic, 'ls_forcing_top_v',       &
2646                                   t+1,                                        &
2647                                   nysv, nyn-nysv+1,                           &
2648                                   i+1, 1,                                     &
2649                                   force%v_top(t-force%tind,nysv:nyn,i) )
2650                ENDDO
2651             ENDDO
2652             DO  i = nxl, nxr
2653                DO  t = force%tind, force%tind_p
2654                   CALL get_variable_bc( id_dynamic, 'ls_forcing_top_w',       &
2655                                   t+1,                                        &
2656                                   nys+1, nyn-nys+1,                           &
2657                                   i+1, 1,                                     &
2658                                   force%w_top(t-force%tind,nys:nyn,i) )
2659                ENDDO
2660             ENDDO
2661             IF ( .NOT. neutral )  THEN
2662                DO  i = nxl, nxr
2663                   DO  t = force%tind, force%tind_p
2664                      CALL get_variable_bc( id_dynamic, 'ls_forcing_top_pt',   &
2665                                      t+1,                                     &
2666                                      nys+1, nyn-nys+1,                        &
2667                                      i+1, 1,                                  &
2668                                      force%pt_top(t-force%tind,nys:nyn,i) )
2669                   ENDDO
2670                ENDDO
2671             ENDIF
2672             IF ( humidity )  THEN
2673                DO  i = nxl, nxr
2674                   DO  t = force%tind, force%tind_p
2675                      CALL get_variable_bc( id_dynamic, 'ls_forcing_top_qv',   &
2676                                      t+1,                                     &
2677                                      nys+1, nyn-nys+1,                        &
2678                                      i+1, 1,                                  &
2679                                      force%q_top(t-force%tind,nys:nyn,i) )
2680                   ENDDO
2681                ENDDO
2682             ENDIF
2683
2684!
2685!--          Close input file
2686             CALL close_input_file( id_dynamic )
2687#endif
2688          ENDIF
2689#if defined( __parallel )
2690          CALL MPI_BARRIER( comm2d, ierr )
2691#endif
2692       ENDDO
2693
2694!
2695!--    Finally, after data input set control flag indicating that vertical
2696!--    inter- and/or extrapolation is required.
2697!--    Please note, inter/extrapolation of INIFOR data is only a workaroud,
2698!--    as long as INIFOR delivers vertically equidistant data.
2699       force%interpolated = .FALSE. 
2700
2701    END SUBROUTINE netcdf_data_input_lsf
2702
2703
2704!------------------------------------------------------------------------------!
2705! Description:
2706! ------------
2707!> Checks input file for consistency and minimum requirements.
2708!------------------------------------------------------------------------------!
2709    SUBROUTINE netcdf_data_input_check_dynamic
2710
2711       USE control_parameters,                                                 &
2712           ONLY:  initializing_actions, forcing, message_string
2713
2714       IMPLICIT NONE
2715
2716!
2717!--    In case of forcing, check whether dynamic input file is present
2718       IF ( .NOT. input_pids_dynamic  .AND.  forcing )  THEN
2719          message_string = 'forcing = .TRUE. requires dynamic input file ' //  &
2720                            TRIM( input_file_dynamic ) // TRIM( coupling_char )
2721          CALL message( 'netcdf_data_input_mod', 'PA0430', 1, 2, 0, 6, 0 )
2722       ENDIF
2723!
2724!--    Dynamic input file must also be present if initialization via inifor is
2725!--    prescribed.
2726       IF ( .NOT. input_pids_dynamic  .AND.                                    &
2727            TRIM( initializing_actions ) == 'inifor' )  THEN
2728          message_string = 'initializing_actions = inifor requires dynamic ' //&
2729                           'input file ' // TRIM( input_file_dynamic ) //      &
2730                           TRIM( coupling_char )
2731          CALL message( 'netcdf_data_input_mod', 'PA0430', 1, 2, 0, 6, 0 )
2732       ENDIF
2733
2734    END SUBROUTINE netcdf_data_input_check_dynamic
2735
2736!------------------------------------------------------------------------------!
2737! Description:
2738! ------------
2739!> Checks input file for consistency and minimum requirements.
2740!------------------------------------------------------------------------------!
2741    SUBROUTINE netcdf_data_input_check_static
2742
2743       USE arrays_3d,                                                          &
2744           ONLY:  zu
2745
2746       USE control_parameters,                                                 &
2747           ONLY:  land_surface, message_string, urban_surface
2748
2749       USE indices,                                                            &
2750           ONLY:  nxl, nxr, nyn, nys
2751
2752       IMPLICIT NONE
2753
2754       INTEGER(iwp) ::  i      !< loop index along x-direction
2755       INTEGER(iwp) ::  j      !< loop index along y-direction
2756       INTEGER(iwp) ::  n_surf !< number of different surface types at given location
2757
2758       LOGICAL      ::  check_passed !< flag indicating if a check passed
2759
2760!
2761!--    Return if no static input file is available
2762       IF ( .NOT. input_pids_static )  RETURN 
2763!
2764!--    Check orography for fill-values. For the moment, give an error message.
2765!--    More advanced methods, e.g. a nearest neighbor algorithm as used in GIS
2766!--    systems might be implemented later.
2767       IF ( ANY( terrain_height_f%var == terrain_height_f%fill ) )  THEN
2768          message_string = 'NetCDF variable orography_2D is not ' //           &
2769                           'allowed to have missing data'
2770          CALL message( 'netcdf_data_input_mod', 'PA0999', 2, 2, 0, 6, 0 )
2771       ENDIF
2772!
2773!--    If 3D buildings are read, check if building information is consistent
2774!--    to numeric grid.
2775       IF ( buildings_f%from_file )  THEN
2776          IF ( buildings_f%lod == 2 )  THEN
2777             IF ( buildings_f%nz > SIZE( zu ) )  THEN
2778                message_string = 'Reading 3D building data - too much ' //     &
2779                                 'data points along the vertical coordinate.' 
2780                CALL message( 'netcdf_data_input_mod', 'PA0999', 2, 2, 0, 6, 0 )
2781             ENDIF
2782
2783             IF ( ANY( buildings_f%z(0:buildings_f%nz-1) /=                    &
2784                       zu(0:buildings_f%nz-1) ) )  THEN
2785                message_string = 'Reading 3D building data - vertical ' //     &
2786                                 'coordinate do not match numeric grid.'
2787                CALL message( 'netcdf_data_input_mod', 'PA0999', 2, 2, 0, 6, 0 )
2788             ENDIF
2789          ENDIF
2790       ENDIF
2791
2792!
2793!--    Skip further checks concerning buildings and natural surface properties
2794!--    if no urban surface and land surface model are applied.
2795       IF (  .NOT. land_surface  .OR.  .NOT. urban_surface )  RETURN 
2796!
2797!--    Check for minimum requirement of surface-classification data in case
2798!--    static input file is used.
2799       IF ( .NOT. vegetation_type_f%from_file  .OR.                            &
2800            .NOT. pavement_type_f%from_file    .OR.                            &
2801            .NOT. building_type_f%from_file    .OR.                            &
2802            .NOT. water_type_f%from_file       .OR.                            &
2803            .NOT. soil_type_f%from_file             )  THEN
2804          message_string = 'Minimum requirement for surface classification ' //&
2805                           'is not fulfilled. At least ' //                    &
2806                           'vegetation_type, pavement_type, ' //               &
2807                           'building_type, soil_type and water_type are '//    &
2808                           'required.'
2809          CALL message( 'netcdf_data_input_mod', 'PA0999', 1, 2, 0, 6, 0 )
2810       ENDIF
2811!
2812!--    Check for general availability of input variables.
2813!--    If vegetation_type is 0 at any location, vegetation_pars as well as
2814!--    root_area_density_lsm are required.
2815       IF ( vegetation_type_f%from_file )  THEN
2816          IF ( ANY( vegetation_type_f%var == 0 ) )  THEN
2817             IF ( .NOT. vegetation_pars_f%from_file )  THEN
2818                message_string = 'If vegegation_type = 0 at any location, ' // &
2819                                 'vegetation_pars is required'
2820                CALL message( 'netcdf_data_input_mod', 'PA0999', 2, 2, 0, 6, 0 )
2821             ENDIF
2822             IF ( .NOT. root_area_density_lsm_f%from_file )  THEN
2823                message_string = 'If vegegation_type = 0 at any location, ' // &
2824                                 'root_area_density_lsm is required'
2825                CALL message( 'netcdf_data_input_mod', 'PA0999', 2, 2, 0, 6, 0 )
2826             ENDIF
2827          ENDIF
2828       ENDIF
2829!
2830!--    If soil_type is zero at any location, soil_pars is required.
2831       IF ( soil_type_f%from_file )  THEN
2832          check_passed = .TRUE. 
2833          IF ( ALLOCATED( soil_type_f%var_2d ) )  THEN
2834             IF ( ANY( soil_type_f%var_2d == 0 ) )  THEN
2835                IF ( .NOT. soil_pars_f%from_file )  check_passed = .FALSE.
2836             ENDIF
2837          ELSE
2838             IF ( ANY( soil_type_f%var_3d == 0 ) )  THEN
2839                IF ( .NOT. soil_pars_f%from_file )  check_passed = .FALSE.
2840             ENDIF
2841          ENDIF
2842          IF ( .NOT. check_passed )  THEN
2843             message_string = 'If soil_type = 0 at any location, ' //          &
2844                              'soil_pars is required'
2845             CALL message( 'netcdf_data_input_mod', 'PA0999', 2, 2, 0, 6, 0 )         
2846          ENDIF
2847       ENDIF
2848!
2849!--    If building_type is zero at any location, building_pars is required.
2850       IF ( building_type_f%from_file )  THEN
2851          IF ( ANY( building_type_f%var == 0 ) )  THEN
2852             IF ( .NOT. building_pars_f%from_file )  THEN
2853                message_string = 'If building_type = 0 at any location, ' //   &
2854                                 'building_pars is required'
2855                CALL message( 'netcdf_data_input_mod', 'PA0999', 2, 2, 0, 6, 0 )         
2856             ENDIF
2857          ENDIF
2858       ENDIF
2859!
2860!--    If albedo_type is zero at any location, albedo_pars is required.
2861       IF ( albedo_type_f%from_file )  THEN
2862          IF ( ANY( albedo_type_f%var == 0 ) )  THEN
2863             IF ( .NOT. albedo_pars_f%from_file )  THEN
2864                message_string = 'If albedo_type = 0 at any location, ' //     &
2865                                 'albedo_pars is required'
2866                CALL message( 'netcdf_data_input_mod', 'PA0999', 2, 2, 0, 6, 0 )     
2867             ENDIF     
2868          ENDIF
2869       ENDIF
2870!
2871!--    If pavement_type is zero at any location, pavement_pars is required.
2872       IF ( pavement_type_f%from_file )  THEN
2873          IF ( ANY( pavement_type_f%var == 0 ) )  THEN
2874             IF ( .NOT. pavement_pars_f%from_file )  THEN
2875                message_string = 'If pavement_type = 0 at any location, ' //   &
2876                                 'pavement_pars is required'
2877                CALL message( 'netcdf_data_input_mod', 'PA0999', 2, 2, 0, 6, 0 )     
2878             ENDIF     
2879          ENDIF
2880       ENDIF
2881!
2882!--    If pavement_type is zero at any location, also pavement_subsurface_pars
2883!--    is required.
2884       IF ( pavement_type_f%from_file )  THEN
2885          IF ( ANY( pavement_type_f%var == 0 ) )  THEN
2886             IF ( .NOT. pavement_subsurface_pars_f%from_file )  THEN
2887                message_string = 'If pavement_type = 0 at any location, ' //   &
2888                                 'pavement_subsurface_pars is required'
2889                CALL message( 'netcdf_data_input_mod', 'PA0999', 2, 2, 0, 6, 0 )     
2890             ENDIF     
2891          ENDIF
2892       ENDIF
2893!
2894!--    If water_type is zero at any location, water_pars is required.
2895       IF ( water_type_f%from_file )  THEN
2896          IF ( ANY( water_type_f%var == 0 ) )  THEN
2897             IF ( .NOT. water_pars_f%from_file )  THEN
2898                message_string = 'If water_type = 0 at any location, ' //      &
2899                                 'water_pars is required'
2900                CALL message( 'netcdf_data_input_mod', 'PA0999', 2, 2, 0, 6, 0 )     
2901             ENDIF     
2902          ENDIF
2903       ENDIF
2904!
2905!--    Check for local consistency of the input data.
2906       DO  i = nxl, nxr
2907          DO  j = nys, nyn
2908!
2909!--          For each (y,x)-location at least one of the parameters
2910!--          vegetation_type, pavement_type, building_type, or water_type
2911!--          must be set to a non­missing value.
2912             IF ( vegetation_type_f%var(j,i) == vegetation_type_f%fill  .AND.  &
2913                  pavement_type_f%var(j,i)   == pavement_type_f%fill    .AND.  &
2914                  building_type_f%var(j,i)   == building_type_f%fill    .AND.  &
2915                  water_type_f%var(j,i)      == water_type_f%fill )  THEN
2916                message_string = 'At least one of the paramters '       //     &
2917                                 'vegetation_type, pavement_type, '     //     &
2918                                 'building_type, or water_type must be set '// &
2919                                 'to a non-missing value'
2920                CALL message( 'netcdf_data_input_mod', 'PA0999', 2, 2, 0, 6, 0 )
2921             ENDIF
2922!
2923!--          Note that a soil_type is required for each location (y,x) where
2924!--          either vegetation_type or pavement_type is a non­missing value.
2925             IF ( ( vegetation_type_f%var(j,i) /= vegetation_type_f%fill  .OR. &
2926                    pavement_type_f%var(j,i)   /= pavement_type_f%fill ) )  THEN
2927                check_passed = .TRUE.
2928                IF ( ALLOCATED( soil_type_f%var_2d ) )  THEN
2929                   IF ( soil_type_f%var_2d(j,i) == soil_type_f%fill )          &
2930                      check_passed = .FALSE.
2931                ELSE
2932                   IF ( ANY( soil_type_f%var_3d(:,j,i) == soil_type_f%fill) )  &
2933                      check_passed = .FALSE.
2934                ENDIF
2935
2936                IF ( .NOT. check_passed )  THEN
2937                   message_string = 'soil_type is required for each '//        &
2938                                 'location (y,x) where vegetation_type or ' // &
2939                                 'pavement_type is a non-missing value.'
2940                   CALL message( 'netcdf_data_input_mod', 'PA0999',            &
2941                                  2, 2, 0, 6, 0 )
2942                ENDIF
2943             ENDIF
2944!
2945!--          Check for consistency of surface fraction. If more than one type
2946!--          is set, surface fraction need to be given and the sum must not
2947!--          be larger than 1.
2948             n_surf = 0
2949             IF ( vegetation_type_f%var(j,i) /= vegetation_type_f%fill )       &
2950                n_surf = n_surf + 1
2951             IF ( water_type_f%var(j,i)      /= water_type_f%fill )            &
2952                n_surf = n_surf + 1
2953             IF ( pavement_type_f%var(j,i)   /= pavement_type_f%fill )         &
2954                n_surf = n_surf + 1
2955             
2956             IF ( n_surf > 1 )  THEN
2957                IF ( ANY ( surface_fraction_f%frac(:,j,i) ==                   &
2958                     surface_fraction_f%fill ) )  THEN
2959                   message_string = 'If more than one surface type is ' //     &
2960                                 'given at a location, surface_fraction ' //   &
2961                                 'must be provided.'
2962                   CALL message( 'netcdf_data_input_mod', 'PA0999',            &
2963                                  2, 2, 0, 6, 0 )
2964                ENDIF
2965                IF ( SUM ( surface_fraction_f%frac(:,j,i) ) > 1.0_wp )  THEN
2966                   message_string = 'surface_fraction must not exceed 1'
2967                   CALL message( 'netcdf_data_input_mod', 'PA0999',            &
2968                                  2, 2, 0, 6, 0 )
2969                ENDIF
2970             ENDIF
2971!
2972!--          Check vegetation_pars. If vegetation_type is 0, all parameters
2973!--          need to be set, otherwise, single parameters set by
2974!--          vegetation_type can be overwritten.
2975             IF ( vegetation_type_f%from_file )  THEN
2976                IF ( vegetation_type_f%var(j,i) == 0 )  THEN
2977                   IF ( ANY( vegetation_pars_f%pars_xy(:,j,i) ==               &
2978                             vegetation_pars_f%fill ) )  THEN
2979                      message_string = 'If vegetation_type(y,x) = 0, all '  // &
2980                                       'parameters of vegetation_pars at '//   & 
2981                                       'this location must be set.' 
2982                      CALL message( 'netcdf_data_input_mod', 'PA0999',         &
2983                                     2, 2, 0, 6, 0 )
2984                   ENDIF
2985                ENDIF
2986             ENDIF
2987!
2988!--          Check root distribution. If vegetation_type is 0, all levels must
2989!--          be set.
2990             IF ( vegetation_type_f%from_file )  THEN
2991                IF ( vegetation_type_f%var(j,i) == 0 )  THEN
2992                   IF ( ANY( root_area_density_lsm_f%var(:,j,i) ==             &
2993                             root_area_density_lsm_f%fill ) )  THEN
2994                      message_string = 'If vegetation_type(y,x) = 0, all ' //  &
2995                                       'levels of root_area_density_lsm ' //   &
2996                                       'must be set at this location.' 
2997                      CALL message( 'netcdf_data_input_mod', 'PA0999',         &
2998                                     2, 2, 0, 6, 0 )
2999                   ENDIF
3000                ENDIF
3001             ENDIF
3002!
3003!--          Check soil parameters. If soil_type is 0, all parameters 
3004!--          must be set.
3005             IF ( soil_type_f%from_file )  THEN
3006                check_passed = .TRUE.
3007                IF ( ALLOCATED( soil_type_f%var_2d ) )  THEN
3008                   IF ( soil_type_f%var_2d(j,i) == 0 )  THEN
3009                      IF ( ANY( soil_pars_f%pars_xy(:,j,i) ==                  &
3010                                soil_pars_f%fill ) )  check_passed = .FALSE.
3011                   ENDIF
3012                ELSE
3013                   IF ( ANY( soil_type_f%var_3d(:,j,i) == 0 ) )  THEN
3014                      IF ( ANY( soil_pars_f%pars_xy(:,j,i) ==                  &
3015                                soil_pars_f%fill ) )  check_passed = .FALSE.
3016                   ENDIF
3017                ENDIF
3018                IF ( .NOT. check_passed )  THEN
3019                   message_string = 'If soil_type(y,x) = 0, all levels of '  //& 
3020                                    'soil_pars at this location must be set.' 
3021                   CALL message( 'netcdf_data_input_mod', 'PA0999',            &
3022                                  2, 2, 0, 6, 0 )
3023                ENDIF
3024             ENDIF
3025
3026!
3027!--          Check building parameters. If building_type is 0, all parameters 
3028!--          must be set.
3029             IF ( building_type_f%from_file )  THEN
3030                IF ( building_type_f%var(j,i) == 0 )  THEN
3031                   IF ( ANY( building_pars_f%pars_xy(:,j,i) ==                 &
3032                             building_pars_f%fill ) )  THEN
3033                      message_string = 'If building_type(y,x) = 0, all ' //    &
3034                                       'parameters of building_pars at this '//&
3035                                       'location must be set.' 
3036                      CALL message( 'netcdf_data_input_mod', 'PA0999',         &
3037                                     2, 2, 0, 6, 0 )
3038                   ENDIF
3039                ENDIF
3040             ENDIF
3041!
3042!--          Check if building_type is set at each building
3043!              IF ( building_type_f%from_file  .AND.  buildings_f%from_file )  THEN
3044!                 IF ( buildings_f%var_2d(j,i) /= buildings_f%fill1  .AND.       &
3045!                      building_type_f%var(j,i) == building_type_f%fill )  THEN
3046!                       WRITE( message_string, * ) 'Each building requires ' //  &
3047!                                        ' a type. i, j = ', i, j
3048!                       CALL message( 'netcdf_data_input_mod', 'PA0999',         &
3049!                                      2, 2, 0, 6, 0 )
3050!                 ENDIF
3051!              ENDIF
3052!
3053!--          Check albedo parameters. If albedo_type is 0, all parameters 
3054!--          must be set.
3055             IF ( albedo_type_f%from_file )  THEN
3056                IF ( albedo_type_f%var(j,i) == 0 )  THEN
3057                   IF ( ANY( albedo_pars_f%pars_xy(:,j,i) ==                   &
3058                             albedo_pars_f%fill ) )  THEN
3059                      message_string = 'If albedo_type(y,x) = 0, all ' //      &
3060                                       'parameters of albedo_pars at this ' // &
3061                                       'location must be set.' 
3062                      CALL message( 'netcdf_data_input_mod', 'PA0999',         &
3063                                     2, 2, 0, 6, 0 )
3064                   ENDIF
3065                ENDIF
3066             ENDIF
3067
3068!
3069!--          Check pavement parameters. If pavement_type is 0, all parameters 
3070!--          of pavement_pars must be set at this location.
3071             IF ( pavement_type_f%from_file )  THEN
3072                IF ( pavement_type_f%var(j,i) == 0 )  THEN
3073                   IF ( ANY( pavement_pars_f%pars_xy(:,j,i) ==                 &
3074                             pavement_pars_f%fill ) )  THEN
3075                      message_string = 'If pavement_type(y,x) = 0, all ' //    &
3076                                       'parameters of pavement_pars at this '//&
3077                                       'location must be set.' 
3078                      CALL message( 'netcdf_data_input_mod', 'PA0999',         &
3079                                     2, 2, 0, 6, 0 )
3080                   ENDIF
3081                ENDIF
3082             ENDIF
3083!
3084!--          Check pavement-subsurface parameters. If pavement_type is 0,
3085!--          all parameters of pavement_subsurface_pars must be set  at this
3086!--          location.
3087             IF ( pavement_type_f%from_file )  THEN
3088                IF ( pavement_type_f%var(j,i) == 0 )  THEN
3089                   IF ( ANY( pavement_subsurface_pars_f%pars_xyz(:,:,j,i) ==   &
3090                             pavement_subsurface_pars_f%fill ) )  THEN
3091                      message_string = 'If pavement_type(y,x) = 0, all ' //    &
3092                                       'parameters of '                  //    &
3093                                       'pavement_subsurface_pars at this '//   &
3094                                       'location must be set.' 
3095                      CALL message( 'netcdf_data_input_mod', 'PA0999',         &
3096                                     2, 2, 0, 6, 0 )
3097                   ENDIF
3098                ENDIF
3099             ENDIF
3100
3101!
3102!--          Check water parameters. If water_type is 0, all parameters 
3103!--          must be set  at this location.
3104             IF ( water_type_f%from_file )  THEN
3105                IF ( water_type_f%var(j,i) == 0 )  THEN
3106                   IF ( ANY( water_pars_f%pars_xy(:,j,i) ==                    &
3107                             water_pars_f%fill ) )  THEN
3108                      message_string = 'If water_type(y,x) = 0, all ' //       &
3109                                       'parameters of water_pars at this ' //  &
3110                                       'location must be set.' 
3111                      CALL message( 'netcdf_data_input_mod', 'PA0999',         &
3112                                     2, 2, 0, 6, 0 )
3113                   ENDIF
3114                ENDIF
3115             ENDIF
3116
3117          ENDDO
3118       ENDDO
3119
3120    END SUBROUTINE netcdf_data_input_check_static
3121
3122!------------------------------------------------------------------------------!
3123! Description:
3124! ------------
3125!> Vertical interpolation and extrapolation of 1D variables.
3126!------------------------------------------------------------------------------!
3127    SUBROUTINE netcdf_data_input_interpolate_1d( var, z_grid, z_file)
3128
3129       IMPLICIT NONE
3130
3131       LOGICAL      ::  top     !< flag indicating extrapolation at model top
3132
3133       INTEGER(iwp) ::  k       !< running index z-direction file
3134       INTEGER(iwp) ::  kk      !< running index z-direction stretched model grid
3135       INTEGER(iwp) ::  kl      !< lower index bound along z-direction
3136       INTEGER(iwp) ::  ku      !< upper index bound along z-direction
3137       INTEGER(iwp) ::  nz_file !< number of vertical levels on file
3138
3139
3140       REAL(wp), DIMENSION(:) ::  z_grid                  !< grid levels on numeric grid
3141       REAL(wp), DIMENSION(:) ::  z_file                  !< grid levels on file grid
3142       REAL(wp), DIMENSION(:), INTENT(INOUT) ::  var      !< treated variable
3143       REAL(wp), DIMENSION(:), ALLOCATABLE   ::  var_tmp  !< temporary variable
3144
3145
3146       kl = LBOUND(var,1)
3147       ku = UBOUND(var,1)
3148       ALLOCATE( var_tmp(kl:ku) )
3149
3150       DO  k = kl, ku
3151
3152          kk = MINLOC( ABS( z_file - z_grid(k) ), DIM = 1 )
3153
3154          IF ( kk < ku )  THEN
3155             IF ( z_file(kk) - z_grid(k) <= 0.0_wp )  THEN
3156                var_tmp(k) = var(kk) +                                         &
3157                                       ( var(kk+1)        - var(kk)    ) /     &
3158                                       ( z_file(kk+1)     - z_file(kk) ) *     &
3159                                       ( z_grid(k)        - z_file(kk) ) 
3160
3161             ELSEIF ( z_file(kk) - z_grid(k) > 0.0_wp )  THEN
3162                var_tmp(k) = var(kk-1) +                                       &
3163                                         ( var(kk)     - var(kk-1)    ) /      &
3164                                         ( z_file(kk)  - z_file(kk-1) ) *      &
3165                                         ( z_grid(k)   - z_file(kk-1) ) 
3166             ENDIF
3167!
3168!--       Extrapolate
3169          ELSE
3170     
3171             var_tmp(k) = var(ku) +   ( var(ku)    - var(ku-1)      ) /        &
3172                                      ( z_file(ku) - z_file(ku-1)   ) *        &
3173                                      ( z_grid(k)  - z_file(ku)     ) 
3174   
3175          ENDIF
3176
3177       ENDDO
3178       var(:) = var_tmp(:)
3179
3180       DEALLOCATE( var_tmp )
3181
3182
3183    END SUBROUTINE netcdf_data_input_interpolate_1d
3184
3185
3186!------------------------------------------------------------------------------!
3187! Description:
3188! ------------
3189!> Vertical interpolation and extrapolation of 1D variables from Inifor grid
3190!> onto Palm grid, where both have same dimension. Please note, the passed
3191!> paramter list in 1D version is different compared to 2D version.
3192!------------------------------------------------------------------------------!
3193    SUBROUTINE netcdf_data_input_interpolate_1d_soil( var, var_file,           &
3194                                                      z_grid, z_file,          &
3195                                                      nzb_var, nzt_var,        & 
3196                                                      nzb_file, nzt_file )
3197
3198       IMPLICIT NONE
3199
3200       INTEGER(iwp) ::  i        !< running index x-direction
3201       INTEGER(iwp) ::  j        !< running index y-direction
3202       INTEGER(iwp) ::  k        !< running index z-direction file
3203       INTEGER(iwp) ::  kk       !< running index z-direction stretched model grid
3204       INTEGER(iwp) ::  ku       !< upper index bound along z-direction for varialbe from file
3205       INTEGER(iwp) ::  nzb_var  !< lower bound of final array
3206       INTEGER(iwp) ::  nzt_var  !< upper bound of final array
3207       INTEGER(iwp) ::  nzb_file !< lower bound of file array
3208       INTEGER(iwp) ::  nzt_file !< upper bound of file array
3209
3210!        LOGICAL, OPTIONAL ::  depth !< flag indicating reverse z-axis, i.e. depth instead of height, e.g. in case of ocean or soil
3211
3212       REAL(wp), DIMENSION(nzb_var:nzt_var)   ::  z_grid   !< grid levels on numeric grid
3213       REAL(wp), DIMENSION(nzb_file:nzt_file) ::  z_file   !< grid levels on file grid
3214       REAL(wp), DIMENSION(nzb_var:nzt_var)   ::  var      !< treated variable
3215       REAL(wp), DIMENSION(nzb_file:nzt_file) ::  var_file !< temporary variable
3216
3217       ku = nzt_file
3218
3219       DO  k = nzb_var, nzt_var
3220!
3221!--       Determine index on Inifor grid which is closest to the actual height
3222          kk = MINLOC( ABS( z_file - z_grid(k) ), DIM = 1 )
3223!
3224!--       If closest index on Inifor grid is smaller than top index,
3225!--       interpolate the data
3226          IF ( kk < nzt_file )  THEN
3227             IF ( z_file(kk) - z_grid(k) <= 0.0_wp )  THEN
3228                var(k) = var_file(kk) + ( var_file(kk+1) - var_file(kk) ) /    &
3229                                        ( z_file(kk+1)   - z_file(kk)   ) *    &
3230                                        ( z_grid(k)      - z_file(kk)   ) 
3231
3232             ELSEIF ( z_file(kk) - z_grid(k) > 0.0_wp )  THEN
3233                var(k) = var_file(kk-1) + ( var_file(kk) - var_file(kk-1) ) /  &
3234                                          ( z_file(kk)   - z_file(kk-1)   ) *  &
3235                                          ( z_grid(k)    - z_file(kk-1)   ) 
3236             ENDIF
3237!
3238!--       Extrapolate if actual height is above the highest Inifor level
3239          ELSE
3240             var(k) = var_file(ku) + ( var_file(ku) - var_file(ku-1) ) /       &
3241                                     ( z_file(ku)   - z_file(ku-1)   ) *       &
3242                                     ( z_grid(k)    - z_file(ku)     ) 
3243
3244          ENDIF
3245
3246       ENDDO
3247
3248    END SUBROUTINE netcdf_data_input_interpolate_1d_soil
3249
3250!------------------------------------------------------------------------------!
3251! Description:
3252! ------------
3253!> Vertical interpolation and extrapolation of 2D variables at lateral boundaries.
3254!------------------------------------------------------------------------------!
3255    SUBROUTINE netcdf_data_input_interpolate_2d( var, z_grid, z_file)
3256
3257       IMPLICIT NONE
3258
3259       LOGICAL      ::  top     !< flag indicating extrapolation at model top
3260
3261       INTEGER(iwp) ::  i       !< running index x- or y -direction
3262       INTEGER(iwp) ::  il      !< lower index bound along x- or y-direction
3263       INTEGER(iwp) ::  iu      !< upper index bound along x- or y-direction
3264       INTEGER(iwp) ::  k       !< running index z-direction file
3265       INTEGER(iwp) ::  kk      !< running index z-direction stretched model grid
3266       INTEGER(iwp) ::  kl      !< lower index bound along z-direction
3267       INTEGER(iwp) ::  ku      !< upper index bound along z-direction
3268       INTEGER(iwp) ::  nz_file !< number of vertical levels on file
3269
3270
3271       REAL(wp), DIMENSION(:) ::  z_grid                  !< grid levels on numeric grid
3272       REAL(wp), DIMENSION(:) ::  z_file                  !< grid levels on file grid
3273       REAL(wp), DIMENSION(:,:), INTENT(INOUT) ::  var    !< treated variable
3274       REAL(wp), DIMENSION(:), ALLOCATABLE   ::  var_tmp  !< temporary variable
3275
3276
3277       il = LBOUND(var,2)
3278       iu = UBOUND(var,2)
3279       kl = LBOUND(var,1)
3280       ku = UBOUND(var,1)
3281       ALLOCATE( var_tmp(kl:ku) )
3282
3283       DO  i = il, iu
3284          DO  k = kl, ku
3285
3286             kk = MINLOC( ABS( z_file - z_grid(k) ), DIM = 1 )
3287
3288             IF ( kk < ku )  THEN
3289                IF ( z_file(kk) - z_grid(k) <= 0.0_wp )  THEN
3290                   var_tmp(k) = var(kk,i) +                                    &
3291                                          ( var(kk+1,i)      - var(kk,i)  ) /  &
3292                                          ( z_file(kk+1)     - z_file(kk) ) *  &
3293                                          ( z_grid(k)        - z_file(kk) ) 
3294
3295                ELSEIF ( z_file(kk) - z_grid(k) > 0.0_wp )  THEN
3296                   var_tmp(k) = var(kk-1,i) +                                  &
3297                                            ( var(kk,i)   - var(kk-1,i)  ) /   &
3298                                            ( z_file(kk)  - z_file(kk-1) ) *   &
3299                                            ( z_grid(k)   - z_file(kk-1) ) 
3300                ENDIF
3301!
3302!--          Extrapolate
3303             ELSE
3304     
3305                var_tmp(k) = var(ku,i) + ( var(ku,i)  - var(ku-1,i)    ) /     &
3306                                         ( z_file(ku) - z_file(ku-1)   ) *     &
3307                                         ( z_grid(k)  - z_file(ku)     ) 
3308   
3309             ENDIF
3310
3311          ENDDO
3312          var(:,i) = var_tmp(:)
3313
3314       ENDDO
3315
3316       DEALLOCATE( var_tmp )
3317
3318
3319    END SUBROUTINE netcdf_data_input_interpolate_2d
3320
3321!------------------------------------------------------------------------------!
3322! Description:
3323! ------------
3324!> Vertical interpolation and extrapolation of 3D variables.
3325!------------------------------------------------------------------------------!
3326    SUBROUTINE netcdf_data_input_interpolate_3d( var, z_grid, z_file )
3327
3328       IMPLICIT NONE
3329
3330       INTEGER(iwp) ::  i       !< running index x-direction
3331       INTEGER(iwp) ::  il      !< lower index bound along x-direction
3332       INTEGER(iwp) ::  iu      !< upper index bound along x-direction
3333       INTEGER(iwp) ::  j       !< running index y-direction
3334       INTEGER(iwp) ::  jl      !< lower index bound along x-direction
3335       INTEGER(iwp) ::  ju      !< upper index bound along x-direction
3336       INTEGER(iwp) ::  k       !< running index z-direction file
3337       INTEGER(iwp) ::  kk      !< running index z-direction stretched model grid
3338       INTEGER(iwp) ::  kl      !< lower index bound along z-direction
3339       INTEGER(iwp) ::  ku      !< upper index bound along z-direction
3340       INTEGER(iwp) ::  nz_file !< number of vertical levels on file
3341
3342       REAL(wp), DIMENSION(:) ::  z_grid                      !< grid levels on numeric grid
3343       REAL(wp), DIMENSION(:) ::  z_file                      !< grid levels on file grid
3344       REAL(wp), DIMENSION(:,:,:), INTENT(INOUT) ::  var      !< treated variable
3345       REAL(wp), DIMENSION(:), ALLOCATABLE       ::  var_tmp  !< temporary variable
3346
3347       il = LBOUND(var,3)
3348       iu = UBOUND(var,3)
3349       jl = LBOUND(var,2)
3350       ju = UBOUND(var,2)
3351       kl = LBOUND(var,1)
3352       ku = UBOUND(var,1)
3353
3354       ALLOCATE( var_tmp(kl:ku) )
3355
3356       DO  i = il, iu
3357          DO  j = jl, ju
3358             DO  k = kl, ku
3359
3360                kk = MINLOC( ABS( z_file - z_grid(k) ), DIM = 1 ) 
3361
3362                IF ( kk < ku )  THEN
3363                   IF ( z_file(kk) - z_grid(k) <= 0.0_wp )  THEN
3364                      var_tmp(k) = var(kk,j,i) +                               &
3365                                             ( var(kk+1,j,i) - var(kk,j,i) ) / &
3366                                             ( z_file(kk+1)  - z_file(kk)  ) * &
3367                                             ( z_grid(k)     - z_file(kk)  ) 
3368
3369                   ELSEIF ( z_file(kk) - z_grid(k) > 0.0_wp )  THEN
3370                      var_tmp(k) = var(kk-1,j,i) +                             &
3371                                             ( var(kk,j,i) - var(kk-1,j,i) ) / &
3372                                             ( z_file(kk)  - z_file(kk-1)  ) * &
3373                                             ( z_grid(k)   - z_file(kk-1)  ) 
3374                   ENDIF
3375!
3376!--             Extrapolate
3377                ELSE
3378                   var_tmp(k) = var(ku,j,i) +                                  &
3379                                       ( var(ku,j,i)  - var(ku-1,j,i)   ) /    &
3380                                       ( z_file(ku)   - z_file(ku-1)    ) *    &
3381                                       ( z_grid(k)    - z_file(ku)      ) 
3382
3383                ENDIF
3384             ENDDO
3385             var(:,j,i) = var_tmp(:)
3386          ENDDO
3387       ENDDO
3388
3389       DEALLOCATE( var_tmp )
3390
3391
3392    END SUBROUTINE netcdf_data_input_interpolate_3d
3393
3394!------------------------------------------------------------------------------!
3395! Description:
3396! ------------
3397!> Checks if a given variables is on file
3398!------------------------------------------------------------------------------!
3399    FUNCTION check_existence( vars_in_file, var_name )
3400
3401       IMPLICIT NONE
3402
3403       CHARACTER(LEN=*) ::  var_name                   !< variable to be checked
3404       CHARACTER(LEN=*), DIMENSION(:) ::  vars_in_file !< list of variables in file
3405
3406       INTEGER(iwp) ::  i                              !< loop variable
3407
3408       LOGICAL ::  check_existence                     !< flag indicating whether a variable exist or not - actual return value
3409
3410       i = 1
3411       check_existence = .FALSE.
3412       DO  WHILE ( i <= SIZE( vars_in_file ) )
3413          check_existence = TRIM( vars_in_file(i) ) == TRIM( var_name )  .OR.  &
3414                            check_existence
3415          i = i + 1
3416       ENDDO
3417
3418       RETURN
3419
3420    END FUNCTION check_existence
3421
3422
3423!------------------------------------------------------------------------------!
3424! Description:
3425! ------------
3426!> Closes an existing netCDF file.
3427!------------------------------------------------------------------------------!
3428    SUBROUTINE close_input_file( id )
3429#if defined( __netcdf )
3430
3431       USE pegrid
3432
3433       IMPLICIT NONE
3434
3435       INTEGER(iwp), INTENT(INOUT)        ::  id        !< file id
3436
3437       nc_stat = NF90_CLOSE( id )
3438       CALL handle_error( 'close', 537 )
3439#endif
3440    END SUBROUTINE close_input_file
3441
3442!------------------------------------------------------------------------------!
3443! Description:
3444! ------------
3445!> Opens an existing netCDF file for reading only and returns its id.
3446!------------------------------------------------------------------------------!
3447    SUBROUTINE open_read_file( filename, id )
3448#if defined( __netcdf )
3449
3450       USE pegrid
3451
3452       IMPLICIT NONE
3453
3454       CHARACTER (LEN=*), INTENT(IN) ::  filename  !< filename
3455       INTEGER(iwp), INTENT(INOUT)   ::  id        !< file id
3456       LOGICAL                       ::  file_open = .FALSE.
3457
3458       nc_stat = NF90_OPEN( filename, NF90_NOWRITE, id )
3459
3460       CALL handle_error( 'open_read_file', 536 )
3461
3462#endif
3463    END SUBROUTINE open_read_file
3464
3465!------------------------------------------------------------------------------!
3466! Description:
3467! ------------
3468!> Reads global or variable-related attributes of type INTEGER (32-bit)
3469!------------------------------------------------------------------------------!
3470     SUBROUTINE get_attribute_int32( id, attribute_name, value, global,        &
3471                                     variable_name )
3472
3473       USE pegrid
3474
3475       IMPLICIT NONE
3476
3477       CHARACTER(LEN=*)            ::  attribute_name   !< attribute name
3478       CHARACTER(LEN=*), OPTIONAL  ::  variable_name    !< variable name
3479
3480       INTEGER(iwp), INTENT(IN)    ::  id               !< file id
3481       INTEGER(iwp)                ::  id_var           !< variable id
3482       INTEGER(iwp), INTENT(INOUT) ::  value            !< read value
3483
3484       LOGICAL, INTENT(IN) ::  global                   !< flag indicating global attribute
3485#if defined( __netcdf )
3486
3487!
3488!--    Read global attribute
3489       IF ( global )  THEN
3490          nc_stat = NF90_GET_ATT( id, NF90_GLOBAL, TRIM( attribute_name ), value )
3491          CALL handle_error( 'get_attribute_int32 global', 522 )
3492!
3493!--    Read attributes referring to a single variable. Therefore, first inquire
3494!--    variable id
3495       ELSE
3496          nc_stat = NF90_INQ_VARID( id, TRIM( variable_name ), id_var )
3497          CALL handle_error( 'get_attribute_int32', 522 )
3498          nc_stat = NF90_GET_ATT( id, id_var, TRIM( attribute_name ), value )
3499          CALL handle_error( 'get_attribute_int32', 522 )       
3500       ENDIF
3501#endif
3502    END SUBROUTINE get_attribute_int32
3503
3504!------------------------------------------------------------------------------!
3505! Description:
3506! ------------
3507!> Reads global or variable-related attributes of type INTEGER (8-bit)
3508!------------------------------------------------------------------------------!
3509     SUBROUTINE get_attribute_int8( id, attribute_name, value, global,         &
3510                                    variable_name )
3511
3512       USE pegrid
3513
3514       IMPLICIT NONE
3515
3516       CHARACTER(LEN=*)            ::  attribute_name   !< attribute name
3517       CHARACTER(LEN=*), OPTIONAL  ::  variable_name    !< variable name
3518
3519       INTEGER(iwp), INTENT(IN)    ::  id               !< file id
3520       INTEGER(iwp)                ::  id_var           !< variable id
3521       INTEGER(KIND=1), INTENT(INOUT) ::  value         !< read value
3522
3523       LOGICAL, INTENT(IN) ::  global                   !< flag indicating global attribute
3524#if defined( __netcdf )
3525
3526!
3527!--    Read global attribute
3528       IF ( global )  THEN
3529          nc_stat = NF90_GET_ATT( id, NF90_GLOBAL, TRIM( attribute_name ), value )
3530          CALL handle_error( 'get_attribute_int8 global', 523 )
3531!
3532!--    Read attributes referring to a single variable. Therefore, first inquire
3533!--    variable id
3534       ELSE
3535          nc_stat = NF90_INQ_VARID( id, TRIM( variable_name ), id_var )
3536          CALL handle_error( 'get_attribute_int8', 523 )
3537          nc_stat = NF90_GET_ATT( id, id_var, TRIM( attribute_name ), value )
3538          CALL handle_error( 'get_attribute_int8', 523 )       
3539       ENDIF
3540#endif
3541    END SUBROUTINE get_attribute_int8
3542
3543!------------------------------------------------------------------------------!
3544! Description:
3545! ------------
3546!> Reads global or variable-related attributes of type REAL
3547!------------------------------------------------------------------------------!
3548     SUBROUTINE get_attribute_real( id, attribute_name, value, global,         &
3549                                    variable_name )
3550
3551       USE pegrid
3552
3553       IMPLICIT NONE
3554
3555       CHARACTER(LEN=*)            ::  attribute_name   !< attribute name
3556       CHARACTER(LEN=*), OPTIONAL  ::  variable_name    !< variable name
3557
3558       INTEGER(iwp), INTENT(IN)    ::  id               !< file id
3559       INTEGER(iwp)                ::  id_var           !< variable id
3560
3561       LOGICAL, INTENT(IN) ::  global                   !< flag indicating global attribute
3562
3563       REAL(wp), INTENT(INOUT)     ::  value            !< read value
3564#if defined( __netcdf )
3565
3566
3567!
3568!-- Read global attribute
3569       IF ( global )  THEN
3570          nc_stat = NF90_GET_ATT( id, NF90_GLOBAL, TRIM( attribute_name ), value )
3571          CALL handle_error( 'get_attribute_real global', 524 )
3572!
3573!-- Read attributes referring to a single variable. Therefore, first inquire
3574!-- variable id
3575       ELSE
3576          nc_stat = NF90_INQ_VARID( id, TRIM( variable_name ), id_var )
3577          CALL handle_error( 'get_attribute_real', 524 )
3578          nc_stat = NF90_GET_ATT( id, id_var, TRIM( attribute_name ), value )
3579          CALL handle_error( 'get_attribute_real', 524 )       
3580       ENDIF
3581#endif
3582    END SUBROUTINE get_attribute_real
3583
3584!------------------------------------------------------------------------------!
3585! Description:
3586! ------------
3587!> Reads global or variable-related attributes of type CHARACTER
3588!> Remark: reading attributes of type NF_STRING return an error code 56 -
3589!> Attempt to convert between text & numbers.
3590!------------------------------------------------------------------------------!
3591     SUBROUTINE get_attribute_string( id, attribute_name, value, global,       &
3592                                      variable_name )
3593
3594       USE pegrid
3595
3596       IMPLICIT NONE
3597
3598       CHARACTER(LEN=*)                ::  attribute_name   !< attribute name
3599       CHARACTER(LEN=*), OPTIONAL      ::  variable_name    !< variable name
3600       CHARACTER(LEN=*), INTENT(INOUT) ::  value            !< read value
3601
3602       INTEGER(iwp), INTENT(IN)    ::  id               !< file id
3603       INTEGER(iwp)                ::  id_var           !< variable id
3604
3605       LOGICAL, INTENT(IN) ::  global                   !< flag indicating global attribute
3606#if defined( __netcdf )
3607
3608!
3609!--    Read global attribute
3610       IF ( global )  THEN
3611          nc_stat = NF90_GET_ATT( id, NF90_GLOBAL, TRIM( attribute_name ), value )
3612          CALL handle_error( 'get_attribute_string global', 525 )
3613!
3614!--    Read attributes referring to a single variable. Therefore, first inquire
3615!--    variable id
3616       ELSE
3617          nc_stat = NF90_INQ_VARID( id, TRIM( variable_name ), id_var )
3618          CALL handle_error( 'get_attribute_string', 525 )
3619
3620          nc_stat = NF90_GET_ATT( id, id_var, TRIM( attribute_name ), value )
3621          CALL handle_error( 'get_attribute_string',525 ) 
3622
3623       ENDIF
3624#endif
3625    END SUBROUTINE get_attribute_string
3626
3627
3628
3629!------------------------------------------------------------------------------!
3630! Description:
3631! ------------
3632!> Get dimension array for a given dimension
3633!------------------------------------------------------------------------------!
3634     SUBROUTINE get_dimension_length( id, dim_len, variable_name )
3635#if defined( __netcdf )
3636
3637       USE pegrid
3638
3639       IMPLICIT NONE
3640
3641       CHARACTER(LEN=*)            ::  variable_name    !< dimension name
3642       CHARACTER(LEN=100)          ::  dum              !< dummy variable to receive return character
3643
3644       INTEGER(iwp)                ::  dim_len          !< dimension size
3645       INTEGER(iwp), INTENT(IN)    ::  id               !< file id
3646       INTEGER(iwp)                ::  id_dim           !< dimension id
3647
3648!
3649!--    First, inquire dimension ID
3650       nc_stat = NF90_INQ_DIMID( id, TRIM( variable_name ), id_dim )
3651       CALL handle_error( 'get_dimension_length', 526 )
3652!
3653!--    Inquire dimension length
3654       nc_stat = NF90_INQUIRE_DIMENSION( id, id_dim, dum, LEN = dim_len )
3655       CALL handle_error( 'get_dimension_length', 526 ) 
3656
3657#endif
3658    END SUBROUTINE get_dimension_length
3659
3660!------------------------------------------------------------------------------!
3661! Description:
3662! ------------
3663!> Reads a 1D integer variable from file.
3664!------------------------------------------------------------------------------!
3665     SUBROUTINE get_variable_1d_int( id, variable_name, var )
3666
3667       USE pegrid
3668
3669       IMPLICIT NONE
3670
3671       CHARACTER(LEN=*)            ::  variable_name    !< variable name
3672
3673       INTEGER(iwp), INTENT(IN)    ::  id               !< file id
3674       INTEGER(iwp)                ::  id_var           !< dimension id
3675
3676       INTEGER(iwp), DIMENSION(:), INTENT(INOUT) ::  var  !< variable to be read
3677#if defined( __netcdf )
3678
3679!
3680!--    First, inquire variable ID
3681       nc_stat = NF90_INQ_VARID( id, TRIM( variable_name ), id_var )
3682       CALL handle_error( 'get_variable_1d_int', 527 )
3683!
3684!--    Inquire dimension length
3685       nc_stat = NF90_GET_VAR( id, id_var, var )
3686       CALL handle_error( 'get_variable_1d_int', 527 ) 
3687
3688#endif
3689    END SUBROUTINE get_variable_1d_int
3690
3691!------------------------------------------------------------------------------!
3692! Description:
3693! ------------
3694!> Reads a 1D float variable from file.
3695!------------------------------------------------------------------------------!
3696     SUBROUTINE get_variable_1d_real( id, variable_name, var )
3697
3698       USE pegrid
3699
3700       IMPLICIT NONE
3701
3702       CHARACTER(LEN=*)            ::  variable_name    !< variable name
3703
3704       INTEGER(iwp), INTENT(IN)    ::  id               !< file id
3705       INTEGER(iwp)                ::  id_var           !< dimension id
3706
3707       REAL(wp), DIMENSION(:), INTENT(INOUT) ::  var  !< variable to be read
3708#if defined( __netcdf )
3709
3710!
3711!--    First, inquire variable ID
3712       nc_stat = NF90_INQ_VARID( id, TRIM( variable_name ), id_var )
3713       CALL handle_error( 'get_variable_1d_real', 527 )
3714!
3715!--    Inquire dimension length
3716       nc_stat = NF90_GET_VAR( id, id_var, var )
3717       CALL handle_error( 'get_variable_1d_real', 527 ) 
3718
3719#endif
3720    END SUBROUTINE get_variable_1d_real
3721
3722!------------------------------------------------------------------------------!
3723! Description:
3724! ------------
3725!> Reads a 2D REAL variable from a file. Reading is done processor-wise,
3726!> i.e. each core reads its own domain in slices along x.
3727!------------------------------------------------------------------------------!
3728    SUBROUTINE get_variable_2d_real( id, variable_name, i, var )
3729
3730       USE indices
3731       USE pegrid
3732
3733       IMPLICIT NONE
3734
3735       CHARACTER(LEN=*)              ::  variable_name   !< variable name
3736
3737       INTEGER(iwp), INTENT(IN)      ::  i               !< index along x direction
3738       INTEGER(iwp), INTENT(IN)      ::  id              !< file id
3739       INTEGER(iwp)                  ::  id_var          !< variable id
3740
3741       REAL(wp), DIMENSION(nys:nyn), INTENT(INOUT) ::  var  !< variable to be read
3742#if defined( __netcdf )
3743!
3744!--    Inquire variable id
3745       nc_stat = NF90_INQ_VARID( id, TRIM( variable_name ), id_var )
3746!
3747!--    Get variable
3748       nc_stat = NF90_GET_VAR( id, id_var, var(nys:nyn),                       &
3749                               start = (/ i+1, nys+1 /),                       &
3750                               count = (/ 1, nyn - nys + 1 /) )
3751
3752       CALL handle_error( 'get_variable_2d_real', 528 )
3753#endif
3754    END SUBROUTINE get_variable_2d_real
3755
3756!------------------------------------------------------------------------------!
3757! Description:
3758! ------------
3759!> Reads a 2D 32-bit INTEGER variable from file. Reading is done processor-wise,
3760!> i.e. each core reads its own domain in slices along x.
3761!------------------------------------------------------------------------------!
3762    SUBROUTINE get_variable_2d_int32( id, variable_name, i, var )
3763
3764       USE indices
3765       USE pegrid
3766
3767       IMPLICIT NONE
3768
3769       CHARACTER(LEN=*)              ::  variable_name   !< variable name
3770
3771       INTEGER(iwp), INTENT(IN)      ::  i               !< index along x direction
3772       INTEGER(iwp), INTENT(IN)      ::  id              !< file id
3773       INTEGER(iwp)                  ::  id_var          !< variable id
3774       INTEGER(iwp), DIMENSION(nys:nyn), INTENT(INOUT) ::  var  !< variable to be read
3775#if defined( __netcdf )
3776!
3777!--    Inquire variable id
3778       nc_stat = NF90_INQ_VARID( id, TRIM( variable_name ), id_var )
3779!
3780!--    Get variable
3781       nc_stat = NF90_GET_VAR( id, id_var, var(nys:nyn),                       &
3782                               start = (/ i+1, nys+1 /),                       &
3783                               count = (/ 1, nyn - nys + 1 /) )
3784
3785       CALL handle_error( 'get_variable_2d_int32', 529 )
3786#endif
3787    END SUBROUTINE get_variable_2d_int32
3788
3789!------------------------------------------------------------------------------!
3790! Description:
3791! ------------
3792!> Reads a 2D 8-bit INTEGER variable from file. Reading is done processor-wise,
3793!> i.e. each core reads its own domain in slices along x.
3794!------------------------------------------------------------------------------!
3795    SUBROUTINE get_variable_2d_int8( id, variable_name, i, var )
3796
3797       USE indices
3798       USE pegrid
3799
3800       IMPLICIT NONE
3801
3802       CHARACTER(LEN=*)              ::  variable_name   !< variable name
3803
3804       INTEGER(iwp), INTENT(IN)      ::  i               !< index along x direction
3805       INTEGER(iwp), INTENT(IN)      ::  id              !< file id
3806       INTEGER(iwp)                  ::  id_var          !< variable id
3807       INTEGER(KIND=1), DIMENSION(nys:nyn), INTENT(INOUT) ::  var  !< variable to be read
3808#if defined( __netcdf )
3809!
3810!--    Inquire variable id
3811       nc_stat = NF90_INQ_VARID( id, TRIM( variable_name ), id_var )
3812!
3813!--    Get variable
3814       nc_stat = NF90_GET_VAR( id, id_var, var(nys:nyn),                       &
3815                               start = (/ i+1, nys+1 /),                       &
3816                               count = (/ 1, nyn - nys + 1 /) )
3817
3818       CALL handle_error( 'get_variable_2d_int8', 530 )
3819#endif
3820    END SUBROUTINE get_variable_2d_int8
3821
3822!------------------------------------------------------------------------------!
3823! Description:
3824! ------------
3825!> Reads a 3D 8-bit INTEGER variable from file.
3826!------------------------------------------------------------------------------!
3827    SUBROUTINE get_variable_3d_int8( id, variable_name, i, j, var )
3828
3829       USE indices
3830       USE pegrid
3831
3832       IMPLICIT NONE
3833
3834       CHARACTER(LEN=*)              ::  variable_name   !< variable name
3835
3836       INTEGER(iwp), INTENT(IN)      ::  i               !< index along x direction
3837       INTEGER(iwp), INTENT(IN)      ::  id              !< file id
3838       INTEGER(iwp)                  ::  id_var          !< variable id
3839       INTEGER(iwp), INTENT(IN)      ::  j               !< index along y direction
3840       INTEGER(iwp)                  ::  n_file          !< number of data-points along 3rd dimension
3841
3842       INTEGER(iwp), DIMENSION(1:3)  ::  id_dim
3843
3844       INTEGER( KIND = 1 ), DIMENSION(nzb:nzt+1), INTENT(INOUT) ::  var  !< variable to be read
3845#if defined( __netcdf )
3846
3847!
3848!--    Inquire variable id
3849       nc_stat = NF90_INQ_VARID( id, TRIM( variable_name ), id_var )
3850!
3851!--    Get length of first dimension, required for the count parameter.
3852!--    Therefore, first inquired dimension ids
3853       nc_stat = NF90_INQUIRE_VARIABLE( id, id_var, DIMIDS = id_dim )
3854       nc_stat = NF90_INQUIRE_DIMENSION( id, id_dim(3), LEN = n_file )
3855!
3856!--    Get variable
3857       nc_stat = NF90_GET_VAR( id, id_var, var,                  &
3858                               start = (/ i+1, j+1, 1 /),                      &
3859                               count = (/ 1, 1, n_file /) )
3860
3861       CALL handle_error( 'get_variable_3d_int8', 531 )
3862#endif
3863    END SUBROUTINE get_variable_3d_int8
3864
3865
3866!------------------------------------------------------------------------------!
3867! Description:
3868! ------------
3869!> Reads a 3D float variable from file. 
3870!------------------------------------------------------------------------------!
3871    SUBROUTINE get_variable_3d_real( id, variable_name, i, j, var )
3872
3873       USE indices
3874       USE pegrid
3875
3876       IMPLICIT NONE
3877
3878       CHARACTER(LEN=*)              ::  variable_name   !< variable name
3879
3880       INTEGER(iwp), INTENT(IN)      ::  i               !< index along x direction
3881       INTEGER(iwp), INTENT(IN)      ::  id              !< file id
3882       INTEGER(iwp)                  ::  id_var          !< variable id
3883       INTEGER(iwp), INTENT(IN)      ::  j               !< index along y direction
3884       INTEGER(iwp)                  ::  n3              !< number of data-points along 3rd dimension
3885
3886       INTEGER(iwp), DIMENSION(3)    ::  id_dim
3887
3888       REAL(wp), DIMENSION(:), INTENT(INOUT) ::  var     !< variable to be read
3889#if defined( __netcdf )
3890
3891!
3892!--    Inquire variable id
3893       nc_stat = NF90_INQ_VARID( id, TRIM( variable_name ), id_var )
3894!
3895!--    Get length of first dimension, required for the count parameter.
3896!--    Therefore, first inquired dimension ids
3897       nc_stat = NF90_INQUIRE_VARIABLE( id, id_var, DIMIDS = id_dim )
3898       nc_stat = NF90_INQUIRE_DIMENSION( id, id_dim(3), LEN = n3 )
3899!
3900!--    Get variable
3901       nc_stat = NF90_GET_VAR( id, id_var, var,                                &
3902                               start = (/ i+1, j+1, 1 /),                      &
3903                               count = (/ 1, 1, n3 /) )
3904
3905       CALL handle_error( 'get_variable_3d_real', 532 )
3906#endif
3907    END SUBROUTINE get_variable_3d_real
3908
3909
3910!------------------------------------------------------------------------------!
3911! Description:
3912! ------------
3913!> Reads a 4D float variable from file. Note, in constrast to 3D versions,
3914!> dimensions are already inquired and passed so that they are known here.
3915!------------------------------------------------------------------------------!
3916    SUBROUTINE get_variable_4d_real( id, variable_name, i, j, var, n3, n4 )
3917
3918       USE indices
3919       USE pegrid
3920
3921       IMPLICIT NONE
3922
3923       CHARACTER(LEN=*)              ::  variable_name   !< variable name
3924
3925       INTEGER(iwp), INTENT(IN)      ::  i               !< index along x direction
3926       INTEGER(iwp), INTENT(IN)      ::  id              !< file id
3927       INTEGER(iwp)                  ::  id_var          !< variable id
3928       INTEGER(iwp), INTENT(IN)      ::  j               !< index along y direction
3929       INTEGER(iwp), INTENT(IN)      ::  n3              !< number of data-points along 3rd dimension
3930       INTEGER(iwp), INTENT(IN)      ::  n4              !< number of data-points along 4th dimension
3931
3932       INTEGER(iwp), DIMENSION(3)    ::  id_dim
3933
3934       REAL(wp), DIMENSION(:,:), INTENT(INOUT) ::  var     !< variable to be read
3935#if defined( __netcdf )
3936
3937!
3938!--    Inquire variable id
3939       nc_stat = NF90_INQ_VARID( id, TRIM( variable_name ), id_var )
3940!
3941!--    Get variable
3942       nc_stat = NF90_GET_VAR( id, id_var, var,                                &
3943                               start = (/ i+1, j+1, 1, 1 /),                   &
3944                               count = (/ 1, 1, n3, n4 /) )
3945
3946       CALL handle_error( 'get_variable_4d_real', 533 )
3947#endif
3948    END SUBROUTINE get_variable_4d_real
3949
3950
3951
3952!------------------------------------------------------------------------------!
3953! Description:
3954! ------------
3955!> Reads a 3D float variable at left, right, north, south and top boundaries.
3956!------------------------------------------------------------------------------!
3957    SUBROUTINE get_variable_bc( id, variable_name, t_start,                    &
3958                                i2_s, count_2, i3_s, count_3,  var )
3959
3960       USE indices
3961       USE pegrid
3962
3963       IMPLICIT NONE
3964
3965       CHARACTER(LEN=*)              ::  variable_name   !< variable name
3966
3967       INTEGER(iwp)                  ::  count_2         !< number of elements in second dimension
3968       INTEGER(iwp)                  ::  count_3         !< number of elements in third dimension (usually 1)
3969       INTEGER(iwp)                  ::  i2_s            !< start index of second dimension
3970       INTEGER(iwp)                  ::  i3_s            !< start index of third dimension
3971       INTEGER(iwp), INTENT(IN)      ::  id              !< file id
3972       INTEGER(iwp)                  ::  id_var          !< variable id
3973       INTEGER(iwp)                  ::  t_start         !< start index at time dimension with respect to netcdf convention
3974
3975       REAL(wp), DIMENSION(:), INTENT(INOUT) ::  var     !< input variable
3976#if defined( __netcdf )
3977
3978!
3979!--    Inquire variable id
3980       nc_stat = NF90_INQ_VARID( id, TRIM( variable_name ), id_var )
3981!
3982!--    Get variable
3983       nc_stat = NF90_GET_VAR( id, id_var, var,                              &
3984                               start = (/ i3_s, i2_s, t_start /),            & 
3985                               count = (/ count_3, count_2, 1 /) )       
3986
3987       CALL handle_error( 'get_variable_bc', 532 )
3988#endif
3989    END SUBROUTINE get_variable_bc
3990
3991
3992
3993!------------------------------------------------------------------------------!
3994! Description:
3995! ------------
3996!> Inquires the number of variables in a file 
3997!------------------------------------------------------------------------------!
3998    SUBROUTINE inquire_num_variables( id, num_vars )
3999
4000       USE indices
4001       USE pegrid
4002
4003       IMPLICIT NONE
4004
4005       INTEGER(iwp), INTENT(IN)      ::  id              !< file id
4006       INTEGER(iwp), INTENT(INOUT)   ::  num_vars        !< number of variables in a file
4007#if defined( __netcdf )
4008
4009       nc_stat = NF90_INQUIRE( id, NVARIABLES = num_vars )
4010       CALL handle_error( 'inquire_num_variables', 534 )
4011
4012#endif
4013    END SUBROUTINE inquire_num_variables
4014
4015
4016!------------------------------------------------------------------------------!
4017! Description:
4018! ------------
4019!> Inquires the variable names belonging to a file.
4020!------------------------------------------------------------------------------!
4021    SUBROUTINE inquire_variable_names( id, var_names )
4022
4023       USE indices
4024       USE pegrid
4025
4026       IMPLICIT NONE
4027
4028       CHARACTER(LEN=*), DIMENSION(:), INTENT(INOUT) ::  var_names   !< return variable - variable names
4029       INTEGER(iwp)                                  ::  i           !< loop variable
4030       INTEGER(iwp), INTENT(IN)                      ::  id          !< file id
4031       INTEGER(iwp)                                  ::  num_vars    !< number of variables (unused return parameter)
4032       INTEGER(iwp), DIMENSION(:), ALLOCATABLE       ::  varids      !< dummy array to strore variable ids temporarily
4033#if defined( __netcdf )
4034
4035       ALLOCATE( varids(1:SIZE(var_names)) )
4036       nc_stat = NF90_INQ_VARIDS( id, NVARS = num_vars, VARIDS = varids )
4037       CALL handle_error( 'inquire_variable_names', 535 )
4038
4039       DO  i = 1, SIZE(var_names)
4040          nc_stat = NF90_INQUIRE_VARIABLE( id, varids(i), NAME = var_names(i) )
4041          CALL handle_error( 'inquire_variable_names', 535 )
4042       ENDDO
4043
4044       DEALLOCATE( varids )
4045#endif
4046    END SUBROUTINE inquire_variable_names
4047
4048!------------------------------------------------------------------------------!
4049! Description:
4050! ------------
4051!> Prints out a text message corresponding to the current status.
4052!------------------------------------------------------------------------------!
4053    SUBROUTINE handle_error( routine_name, errno )
4054
4055       USE control_parameters,                                                 &
4056           ONLY:  message_string
4057
4058       IMPLICIT NONE
4059
4060       CHARACTER(LEN=6) ::  message_identifier
4061       CHARACTER(LEN=*) ::  routine_name
4062
4063       INTEGER(iwp) ::  errno
4064#if defined( __netcdf )
4065
4066       IF ( nc_stat /= NF90_NOERR )  THEN
4067
4068          WRITE( message_identifier, '(''NC'',I4.4)' )  errno
4069          message_string = TRIM( NF90_STRERROR( nc_stat ) )
4070
4071          CALL message( routine_name, message_identifier, 2, 2, 0, 6, 1 )
4072
4073       ENDIF
4074
4075#endif
4076    END SUBROUTINE handle_error
4077
4078
4079 END MODULE netcdf_data_input_mod
Note: See TracBrowser for help on using the repository browser.