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

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

Check dimensions in static input files, further checks for building type.

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