source: palm/trunk/UTIL/inifor/src/inifor_io.f90 @ 4719

Last change on this file since 4719 was 4675, checked in by eckhard, 4 years ago

Support for homogeneous (domain-averaged) boundary conditions and soil profile initialization

  • Property svn:keywords set to Id
File size: 65.9 KB
Line 
1!> @file src/inifor_io.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 2017-2020 Leibniz Universitaet Hannover
18! Copyright 2017-2020 Deutscher Wetterdienst Offenbach
19!------------------------------------------------------------------------------!
20!
21! Current revisions:
22! -----------------
23!
24!
25! Former revisions:
26! -----------------
27! $Id: inifor_io.f90 4675 2020-09-11 10:00:26Z pavelkrc $
28! New command-line options for soil profile initialization and homogeneous
29!    boundary conditions
30!
31!
32! 4659 2020-08-31 11:21:17Z eckhard
33! New command-line option '--precipitation' enables output of precipitation
34!    surface forcing
35! Only define netCDF variables in enabled IO groups
36! Improved code formatting
37!
38!
39! 4569 2020-06-19 12:46:17Z eckhard
40! Removed unused variable
41!
42!
43! 4568 2020-06-19 11:56:30Z eckhard
44! Handle COSMO soil data with and without additional surface temperature
45!
46!
47! 4553 2020-06-03 16:34:15Z eckhard
48! Option --help now points the user to INIFOR's wiki page
49! Remove deprecated command-line options -clon and -clat
50!
51!
52! 4538 2020-05-18 13:45:35Z eckhard
53! Make setting the vertical PALM origin mandatory
54!
55!
56! 4523 2020-05-07 15:58:16Z eckhard
57! respect integer working precision (iwp) specified in inifor_defs.f90
58!
59!
60! 4481 2020-03-31 18:55:54Z maronga
61! Pass hhl_file directly instead of entire INIFOR configuration
62!
63!
64! 3997 2019-05-23 12:35:57Z eckhard
65! Added boolean indicator for --elevation option invocation
66! Stop INIFOR if no command-line options given
67!
68!
69! 3866 2019-04-05 14:25:01Z eckhard
70! Use PALM's working precision
71! Improved coding style
72!
73!
74! 3801 2019-03-15 17:14:25Z eckhard
75! Added routine get_cosmo_grid() to read in COSMO rotated pole from COSMO domain
76! Moved get_soil_layer_thickness() here from inifor_grid
77!
78! 3785 2019-03-06 10:41:14Z eckhard
79! Temporariliy disabled height-based geostrophic wind averaging
80! Improved variable naming
81!
82!
83! 3764 2019-02-26 13:42:09Z eckhard
84! Removed dependency on radiation input files
85!
86!
87! 3716 2019-02-05 17:02:38Z eckhard
88! Removed dependency on soilmoisture input files
89!
90!
91! 3680 2019-01-18 14:54:12Z knoop
92! Moved get_input_file_list() here from grid module, added check for presence of
93!    input files
94!
95!
96! 3618 2018-12-10 13:25:22Z eckhard
97! Prefixed all INIFOR modules with inifor_
98!
99!
100! 3615 2018-12-10 07:21:03Z raasch
101! bugfix: abort replaced by inifor_abort
102!
103! 3557 2018-11-22 16:01:22Z eckhard
104! Updated documentation, removed unused subroutine write_netcdf_variable_2d()
105!
106!
107! 3537 2018-11-20 10:53:14Z eckhard
108! New routine get_netcdf_dim_vector()
109!
110!
111! 3534 2018-11-19 15:35:16Z raasch
112! bugfix: INTENT attribute changed
113!
114! 3456 2018-10-30 14:29:54Z eckhard
115! NetCDf output of internal arrays only with --debug option
116!
117!
118! 3447 2018-10-29 15:52:54Z eckhard
119! Removed INCLUDE statement for get_netcdf_variable()
120! Renamed source files for compatibilty with PALM build system
121!
122!
123! 3395 2018-10-22 17:32:49Z eckhard
124! Added command-line options for configuring the computation of geostrophic
125!     winds (--averagin-mode, --averaging-angle)
126! Added command-line option --input-prefix for setting input file prefixes all
127!     at once
128! Added --debug option for more verbose terminal output
129! Added option-specific *_is_set LOGICALs to indicate invocation from the
130!     command-line
131! Improved error messages in case of empty file-name strings
132! Improved routine naming
133!
134! 3183 2018-07-27 14:25:55Z suehring
135! Introduced new PALM grid stretching
136! Updated variable names and metadata for PIDS v1.9 compatibility
137! Improved handling of the start date string
138! Better compatibility with older Intel compilers:
139! - avoiding implicit array allocation with new get_netcdf_variable()
140!   subroutine instead of function
141! Improved command line interface:
142! - Added configuration validation
143! - New options to configure input file prefixes
144! - GNU-style short and long option names
145! - Added version and copyright output
146!
147!
148! 3182 2018-07-27 13:36:03Z suehring
149! Initial revision
150!
151!
152!
153! Authors:
154! --------
155!> @author Eckhard Kadasch (Deutscher Wetterdienst, Offenbach)
156!
157! Description:
158! ------------
159!> The io module contains the functions needed to read and write netCDF data in
160!> INIFOR.
161!------------------------------------------------------------------------------!
162#if defined ( __netcdf )
163 MODULE inifor_io
164
165    USE inifor_control
166    USE inifor_defs,                                                           &
167        ONLY:  CFG_INIT_PROFILE, CFG_INIT_VOLUME,                              &
168               CFG_INIT_SOIL_PROFILE, CFG_INIT_SOIL_VOLUME,                    &
169               CFG_FORCING_HETERO, CFG_FORCING_HOMO, CFG_FORCING_NUDGING,      &
170               DATE, SNAME, PATH, PI, TO_RADIANS, TO_DEGREES, VERSION,         &
171               NC_DEPTH_DIM_IDX, NC_DEPTH_NAME, NC_HHL_NAME, NC_RLAT_NAME,     &
172               NC_RLON_NAME, NC_ROTATED_POLE_NAME, NC_POLE_LATITUDE_NAME,      &
173               NC_POLE_LONGITUDE_NAME, RHO_L, iwp, wp,                         &
174               PIDS_ORIGIN_LON, PIDS_ORIGIN_LAT, PIDS_ORIGIN_Z
175    USE inifor_types
176    USE inifor_util,                                                           &
177        ONLY:  add_hours_to, nearly_equal, reverse, str, real_to_str
178    USE netcdf
179
180    IMPLICIT NONE
181
182!------------------------------------------------------------------------------!
183! Description:
184! ------------
185!> get_netcdf_variable() reads the netCDF data and metadate for the netCDF
186!> variable 'in_var%name' from the file 'in_file'. The netCDF data array is
187!> stored in the 'buffer' array and metadata added to the respective members of
188!> the given 'in_var'.
189!------------------------------------------------------------------------------!
190    INTERFACE get_netcdf_variable
191       MODULE PROCEDURE get_netcdf_variable_int
192       MODULE PROCEDURE get_netcdf_variable_real
193    END INTERFACE get_netcdf_variable
194
195    PRIVATE ::  get_netcdf_variable_int, get_netcdf_variable_real
196
197 CONTAINS
198
199!------------------------------------------------------------------------------!
200! Description:
201! ------------
202!> get_netcdf_variable_int() implements the integer variant for the
203!> get_netcdf_variable interface.
204!------------------------------------------------------------------------------!
205 SUBROUTINE get_netcdf_variable_int(in_file, in_var, buffer)
206
207    CHARACTER(LEN=PATH), INTENT(IN)          ::  in_file
208    TYPE(nc_var), INTENT(INOUT)              ::  in_var
209    INTEGER(iwp), ALLOCATABLE, INTENT(INOUT) ::  buffer(:,:,:)
210
211    INTEGER               ::  ncid
212    INTEGER, DIMENSION(3) ::  start, count
213
214    IF ( nf90_open( TRIM(in_file), NF90_NOWRITE, ncid ) .EQ. NF90_NOERR .AND. &
215         nf90_inq_varid( ncid, in_var%name, in_var%varid ) .EQ. NF90_NOERR )  THEN
216
217       CALL get_input_dimensions(in_var, ncid)
218
219       CALL get_netcdf_start_and_count(in_var, start, count)
220       CALL log_runtime('time', 'read')
221
222       ALLOCATE( buffer( count(1), count(2), count(3) ) )
223       CALL log_runtime('time', 'alloc')
224
225       CALL check(nf90_get_var( ncid, in_var%varid, buffer,                  &
226                                start = start,                                 &
227                                count = count ))
228
229    ELSE
230
231       message = "Failed to read '" // TRIM(in_var%name) // &
232          "' from file '" // TRIM(in_file) // "'."
233       CALL inifor_abort('get_netcdf_variable', message)
234
235    ENDIF
236
237    CALL check(nf90_close(ncid))
238    CALL log_runtime('time', 'read')
239
240 END SUBROUTINE get_netcdf_variable_int
241
242
243!------------------------------------------------------------------------------!
244! Description:
245! ------------
246!> get_netcdf_variable_real() implements the real variant for the
247!> get_netcdf_variable interface.
248!------------------------------------------------------------------------------!
249 SUBROUTINE get_netcdf_variable_real(in_file, in_var, buffer)
250
251    CHARACTER(LEN=PATH), INTENT(IN)      ::  in_file
252    TYPE(nc_var), INTENT(INOUT)          ::  in_var
253    REAL(wp), ALLOCATABLE, INTENT(INOUT) ::  buffer(:,:,:)
254
255    INTEGER               ::  ncid
256    INTEGER, DIMENSION(3) ::  start, count
257
258    IF ( nf90_open( TRIM(in_file), NF90_NOWRITE, ncid ) .EQ. NF90_NOERR .AND. &
259         nf90_inq_varid( ncid, in_var%name, in_var%varid ) .EQ. NF90_NOERR )  THEN
260
261       CALL get_input_dimensions(in_var, ncid)
262
263       CALL get_netcdf_start_and_count(in_var, start, count)
264       CALL log_runtime('time', 'read')
265
266       ALLOCATE( buffer( count(1), count(2), count(3) ) )
267       CALL log_runtime('time', 'alloc')
268
269       CALL check(nf90_get_var( ncid, in_var%varid, buffer,                  &
270                                start = start,                                 &
271                                count = count ))
272
273    ELSE
274
275       message = "Failed to read '" // TRIM(in_var%name) // &
276          "' from file '" // TRIM(in_file) // "'."
277       CALL inifor_abort('get_netcdf_variable', message)
278
279    ENDIF
280
281    CALL check(nf90_close(ncid))
282    CALL log_runtime('time', 'read')
283
284 END SUBROUTINE get_netcdf_variable_real
285
286
287!------------------------------------------------------------------------------!
288! Description:
289! ------------
290!> get_netcdf_dim_vector() reads the coordinate array 'coordname' from the
291!> netCDF file 'filename'.
292!------------------------------------------------------------------------------!
293 SUBROUTINE get_netcdf_dim_vector(filename, coordname, coords)
294
295    CHARACTER(LEN=*), INTENT(IN)         ::  filename
296    CHARACTER(LEN=*), INTENT(IN)         ::  coordname
297    REAL(wp), ALLOCATABLE, INTENT(INOUT) ::  coords(:)
298
299    INTEGER ::  ncid, varid, dimlen
300    INTEGER ::  dimids(NF90_MAX_VAR_DIMS)
301
302    IF ( nf90_open( TRIM(filename), NF90_NOWRITE, ncid ) .EQ. NF90_NOERR .AND. &
303         nf90_inq_varid( ncid, coordname, varid ) .EQ. NF90_NOERR )  THEN
304
305       CALL check(nf90_inquire_variable( ncid, varid, dimids = dimids ))
306       CALL check(nf90_inquire_dimension( ncid, dimids(1), len = dimlen ))
307
308       ALLOCATE(coords(dimlen))
309       CALL check(nf90_get_var( ncid, varid, coords))
310
311    ELSE
312
313       message = "Failed to read '" // TRIM(coordname) // &
314          "' from file '" // TRIM(filename) // "'."
315       CALL inifor_abort('get_netcdf_dim_vector', message)
316
317    ENDIF
318
319 END SUBROUTINE get_netcdf_dim_vector
320
321
322!------------------------------------------------------------------------------!
323! Description:
324! ------------
325!> get_input_dimensions() reads dimensions metadata of the netCDF variable given
326!> by 'in_var%name' into 'in_var' data structure.
327!------------------------------------------------------------------------------!
328 SUBROUTINE get_input_dimensions(in_var, ncid)
329
330    TYPE(nc_var), INTENT(INOUT) ::  in_var
331    INTEGER, INTENT(IN)         ::  ncid
332
333    INTEGER ::  i
334
335    CALL check(nf90_get_att( ncid, in_var%varid, "long_name",             &
336                             in_var%long_name))
337
338    CALL check(nf90_get_att( ncid, in_var%varid, "units",                 &
339                             in_var%units ))
340
341    CALL check(nf90_inquire_variable( ncid, in_var%varid,                 &
342                                      ndims  = in_var%ndim,               &
343                                      dimids = in_var%dimids ))
344
345    DO  i = 1, in_var%ndim
346       CALL check(nf90_inquire_dimension( ncid, in_var%dimids(i),         &
347                                          name = in_var%dimname(i),       &
348                                          len  = in_var%dimlen(i) ))
349    ENDDO
350
351 END SUBROUTINE get_input_dimensions
352
353
354!------------------------------------------------------------------------------!
355! Description:
356! ------------
357!> get_netcdf_start_and_count() gets the start position and element counts for
358!> the given netCDF file. This information is used in get_netcdf_variable_int()
359!> and _real() for reading input variables..
360!------------------------------------------------------------------------------!
361 SUBROUTINE get_netcdf_start_and_count(in_var, start, count)
362
363    TYPE(nc_var), INTENT(INOUT)        ::  in_var
364    INTEGER, DIMENSION(3), INTENT(OUT) ::  start, count
365
366    INTEGER ::  ndim
367
368    IF ( in_var%ndim .LT. 2  .OR.  in_var%ndim .GT. 4 )  THEN
369
370       message = "Failed reading NetCDF variable " //                       &
371          TRIM(in_var%name) // " with " //                                  &
372          TRIM(str(INT(in_var%ndim, kind=iwp))) //                          &
373          " dimensions because only two- and and three-dimensional" //      &
374          " variables are supported."
375       CALL inifor_abort('get_netcdf_start_and_count', message)
376
377    ENDIF
378
379    start = (/ 1, 1, 1 /)
380    IF ( TRIM(in_var%name) .EQ. 'T_SO' .AND.                                &
381         in_var%has_redundant_first_level )  THEN
382       
383!
384!--    Skip depth = 0.0 for T_SO and reduce number of depths from 9 to 8
385       in_var%dimlen(3) = in_var%dimlen(3) - 1
386
387!
388!--    Start reading from second level, e.g. depth = 0.005 instead of 0.0
389       start(3) = 2
390    ENDIF
391
392    IF (in_var%ndim .EQ. 2)  THEN
393       in_var%dimlen(3) = 1
394    ENDIF
395
396    ndim = MIN(in_var%ndim, 3)
397    count = (/ 1, 1, 1 /)
398    count(1:ndim) = in_var%dimlen(1:ndim)
399
400 END SUBROUTINE get_netcdf_start_and_count
401
402
403!------------------------------------------------------------------------------!
404! Description:
405! ------------
406!> Routine for defining netCDF variables in the dynamic driver, INIFOR's netCDF
407!> output.
408!------------------------------------------------------------------------------!
409 SUBROUTINE netcdf_define_variable(var, ncid)
410
411     TYPE(nc_var), INTENT(INOUT) ::  var
412     INTEGER, INTENT(IN)         ::  ncid
413
414     CALL check(nf90_def_var(ncid, var%name, NF90_FLOAT,       var%dimids(1:var%ndim), var%varid))
415     CALL check(nf90_put_att(ncid, var%varid, "long_name",     var%long_name))
416     CALL check(nf90_put_att(ncid, var%varid, "units",         var%units))
417     IF ( var%lod .GE. 0 )  THEN
418        CALL check(nf90_put_att(ncid, var%varid, "lod",           var%lod))
419     ENDIF
420     CALL check(nf90_put_att(ncid, var%varid, "source",        var%source))
421     CALL check(nf90_put_att(ncid, var%varid, "_FillValue",    NF90_FILL_REAL))
422
423 END SUBROUTINE netcdf_define_variable
424   
425
426!------------------------------------------------------------------------------!
427! Description:
428! ------------
429!> netcdf_get_dimensions() reads in all dimensions and their lengths and stores
430!> them in the given the 'var' data structure. This information is used later
431!> for writing output variables in update_output().
432!------------------------------------------------------------------------------!
433 SUBROUTINE netcdf_get_dimensions(var, ncid)
434
435     TYPE(nc_var), INTENT(INOUT) ::  var
436     INTEGER, INTENT(IN)         ::  ncid
437     INTEGER                     ::  i
438     CHARACTER(SNAME)            ::  null
439
440     DO  i = 1, var%ndim
441        CALL check(nf90_inquire_dimension(ncid, var%dimids(i), &
442                                          name = null, &
443                                          len  = var%dimlen(i)  ) )
444     ENDDO
445
446 END SUBROUTINE netcdf_get_dimensions
447
448
449!------------------------------------------------------------------------------!
450! Description:
451! ------------
452!> This routine parses and interpretes the command-line options and stores the
453!> resulting settings in the 'cfg' data structure.
454!------------------------------------------------------------------------------!
455 SUBROUTINE parse_command_line_arguments( cfg )
456
457    TYPE(inifor_config), INTENT(INOUT) ::  cfg
458
459    CHARACTER(LEN=PATH)                ::  option, arg
460    INTEGER                            ::  arg_count, i
461
462    cfg%flow_prefix_is_set = .FALSE.
463    cfg%input_prefix_is_set = .FALSE.
464    cfg%p0_is_set = .FALSE.
465    cfg%radiation_prefix_is_set = .FALSE.
466    cfg%soil_prefix_is_set = .FALSE.
467    cfg%precipitation_prefix_is_set = .FALSE.
468    cfg%process_precipitation = .FALSE.
469    cfg%static_driver_is_set = .FALSE.
470    cfg%ug_defined_by_user = .FALSE.
471    cfg%vg_defined_by_user = .FALSE.
472    cfg%z0_is_set = .FALSE.
473
474    arg_count = COMMAND_ARGUMENT_COUNT()
475    IF (arg_count .GT. 0)  THEN
476
477       i = 1
478       DO  WHILE (i .LE. arg_count)
479
480          CALL GET_COMMAND_ARGUMENT( i, option )
481
482          SELECT CASE( TRIM(option) )
483
484             CASE( '--averaging-mode' )
485                CALL get_option_argument( i, arg )
486                cfg%averaging_mode = TRIM(arg)
487
488             CASE( '-date', '-d', '--date' )
489                CALL get_option_argument( i, arg )
490                cfg%start_date = TRIM(arg)
491
492             CASE( '--debug' )
493                cfg%debug = .TRUE.
494
495             CASE( '-z0', '-z', '--elevation' )
496                cfg%z0_is_set = .TRUE.
497                CALL get_option_argument( i, arg )
498                READ(arg, *) cfg%z0
499
500             CASE( '-p0', '-r', '--surface-pressure' )
501                cfg%p0_is_set = .TRUE.
502                CALL get_option_argument( i, arg )
503                READ(arg, *) cfg%p0
504
505             CASE( '-ug', '-u', '--geostrophic-u' )
506                cfg%ug_defined_by_user = .TRUE.
507                CALL get_option_argument( i, arg )
508                READ(arg, *) cfg%ug
509
510             CASE( '-vg', '-v', '--geostrophic-v' )
511                cfg%vg_defined_by_user = .TRUE.
512                CALL get_option_argument( i, arg )
513                READ(arg, *) cfg%vg
514
515             CASE( '-path', '-p', '--path' )
516                CALL get_option_argument( i, arg )
517                 cfg%input_path = TRIM(arg)
518
519             CASE( '-hhl', '-l', '--hhl-file' )
520                CALL get_option_argument( i, arg )
521                cfg%hhl_file = TRIM(arg)
522
523             CASE( '--input-prefix')
524                CALL get_option_argument( i, arg )
525                cfg%input_prefix = TRIM(arg)
526                cfg%input_prefix_is_set = .TRUE.
527   
528             CASE( '-a', '--averaging-angle' )
529                CALL get_option_argument( i, arg )
530                READ(arg, *) cfg%averaging_angle
531
532             CASE( '-static', '-t', '--static-driver' )
533                cfg%static_driver_is_set = .TRUE.
534                CALL get_option_argument( i, arg )
535                cfg%static_driver_file = TRIM(arg)
536
537             CASE( '-soil', '-s', '--soil-file')
538                CALL get_option_argument( i, arg )
539                cfg%soiltyp_file = TRIM(arg)
540
541             CASE( '--flow-prefix')
542                CALL get_option_argument( i, arg )
543                cfg%flow_prefix = TRIM(arg)
544                cfg%flow_prefix_is_set = .TRUE.
545   
546             CASE( '--radiation-prefix')
547                CALL get_option_argument( i, arg )
548                cfg%radiation_prefix = TRIM(arg)
549                cfg%radiation_prefix_is_set = .TRUE.
550   
551             CASE( '--soil-prefix')
552                CALL get_option_argument( i, arg )
553                cfg%soil_prefix = TRIM(arg)
554                cfg%soil_prefix_is_set = .TRUE.
555   
556             CASE( '--precipitation-prefix')
557                CALL get_option_argument( i, arg )
558                cfg%precipitation_prefix = TRIM(arg)
559                cfg%precipitation_prefix_is_set = .TRUE.
560
561             CASE( '--precipitation')
562                cfg%process_precipitation = .TRUE.
563
564             CASE( '-o', '--output' )
565                CALL get_option_argument( i, arg )
566                cfg%output_file = TRIM(arg)
567
568             CASE( '-n', '--namelist' )
569                CALL get_option_argument( i, arg )
570                cfg%namelist_file = TRIM(arg)
571
572             CASE( '-i', '--init-mode' )
573                CALL get_option_argument( i, arg )
574                cfg%ic_mode = TRIM(arg)
575
576             CASE( '--soil-init-mode' )
577                CALL get_option_argument( i, arg )
578                cfg%isc_mode = TRIM(arg)
579
580             CASE( '-f', '--forcing-mode' )
581                CALL get_option_argument( i, arg )
582                cfg%bc_mode = TRIM(arg)
583
584             CASE( '--version' )
585                CALL print_version
586                STOP
587
588             CASE( '--help' )
589                CALL print_version
590                PRINT *, ""
591                PRINT *, &
592                   "For documentation and a list of available command-line options " // NEW_LINE( " " ) // &
593                   " please visit https://palm.muk.uni-hannover.de/trac/wiki/doc/app/iofiles/inifor."
594                STOP
595
596             CASE DEFAULT
597                message = "unknown option '" // TRIM(option) // "'."
598                CALL inifor_abort('parse_command_line_arguments', message)
599
600          END SELECT
601
602          i = i + 1
603
604       ENDDO
605
606    ELSE
607         
608       CALL print_version
609       CALL report( 'parse_command_line_arguments', 'No arguments present, exiting.' )
610       STOP
611
612    ENDIF
613
614 END SUBROUTINE parse_command_line_arguments
615
616
617
618 SUBROUTINE get_datetime_file_list( start_date_string, start_hour, end_hour, &
619                                    step_hour, input_path, prefix, suffix,   &
620                                    file_list )
621
622    CHARACTER (LEN=DATE), INTENT(IN) ::  start_date_string
623    CHARACTER (LEN=*),    INTENT(IN) ::  prefix, suffix, input_path
624    INTEGER(iwp),         INTENT(IN) ::  start_hour, end_hour, step_hour
625    CHARACTER(LEN=*), ALLOCATABLE, INTENT(INOUT) ::  file_list(:)
626
627    INTEGER(iwp)        ::  number_of_intervals, hour, i
628    CHARACTER(LEN=DATE) ::  date_string
629
630    number_of_intervals = CEILING( REAL(end_hour - start_hour) / step_hour )
631    ALLOCATE( file_list(number_of_intervals + 1) )
632
633    DO  i = 0, number_of_intervals
634
635       hour = start_hour + i * step_hour
636       date_string = add_hours_to(start_date_string, hour)
637
638       file_list(i+1) = TRIM(input_path) // TRIM(prefix) //                  &
639                        TRIM(date_string) // TRIM(suffix) // '.nc'
640
641    ENDDO
642
643 END SUBROUTINE get_datetime_file_list
644
645!------------------------------------------------------------------------------!
646! Description:
647! ------------
648!> Establish a list of files based on the given start and end times and file
649!> prefixes and suffixes.
650!------------------------------------------------------------------------------!
651 SUBROUTINE get_input_file_list( start_date_string, start_hour, end_hour,    &
652                                 step_hour, input_path, prefix, suffix,      &
653                                 file_list )
654
655    CHARACTER (LEN=DATE), INTENT(IN) ::  start_date_string
656    CHARACTER (LEN=*),    INTENT(IN) ::  prefix, suffix, input_path
657    INTEGER(iwp),         INTENT(IN) ::  start_hour, end_hour, step_hour
658    CHARACTER(LEN=*), ALLOCATABLE, INTENT(INOUT) ::  file_list(:)
659
660    CALL get_datetime_file_list( start_date_string, start_hour, end_hour,    &
661                                 step_hour, input_path, prefix, suffix,      &
662                                 file_list )
663
664 END SUBROUTINE get_input_file_list
665
666
667!------------------------------------------------------------------------------!
668! Description:
669! ------------
670!> Abort INIFOR if the given file is not present.
671!------------------------------------------------------------------------------!
672LOGICAL FUNCTION file_is_present( filename, file_kind, message )
673
674    CHARACTER(LEN=*), INTENT(IN)            ::  filename, file_kind
675    CHARACTER(LEN=*), INTENT(OUT)           ::  message
676
677    file_is_present = file_present(filename)
678
679    IF (.NOT. file_is_present)  THEN
680
681       IF (LEN( TRIM( filename ) ) == 0)  THEN
682          message = "No name was given for the " // TRIM( file_kind ) // " file."
683       ELSE
684          message = "The " // TRIM( file_kind ) // " file '" //                  &
685                    TRIM( filename ) // "' was not found."
686       ENDIF
687
688    ENDIF
689
690END FUNCTION file_is_present
691
692
693!------------------------------------------------------------------------------!
694! Description:
695! ------------
696!> Abort INIFOR if the given file is not present.
697!------------------------------------------------------------------------------!
698 SUBROUTINE verify_file( file_name, file_kind )
699
700    CHARACTER(LEN=*), INTENT(IN)           ::  file_name, file_kind
701
702    IF (.NOT. file_is_present( file_name, file_kind, message ))  THEN
703
704       CALL inifor_abort( 'verify_file', message )
705
706    ENDIF
707
708    message = "Set up input file name '" // TRIM( file_name ) // "'"
709    CALL report( 'verify_file', message )
710
711 END SUBROUTINE verify_file
712!------------------------------------------------------------------------------!
713! Description:
714! ------------
715!> Get the argument of the i'th command line option, which is at the location
716!> i+1 of the argument list.
717!------------------------------------------------------------------------------!
718 SUBROUTINE get_option_argument( i, arg )
719    CHARACTER(LEN=PATH), INTENT(INOUT) ::  arg
720    INTEGER, INTENT(INOUT)             ::  i
721
722    i = i + 1
723    CALL GET_COMMAND_ARGUMENT( i, arg )
724
725 END SUBROUTINE
726
727
728!------------------------------------------------------------------------------!
729! Description:
730! ------------
731!> Checks the INIFOR configuration 'cfg' for plausibility.
732!------------------------------------------------------------------------------!
733 SUBROUTINE validate_config( cfg )
734    TYPE(inifor_config), INTENT(IN) ::  cfg
735
736    CALL verify_file( cfg%hhl_file, 'HHL' )
737    CALL verify_file( cfg%namelist_file, 'NAMELIST' )
738    CALL verify_file( cfg%soiltyp_file, 'SOILTYP' )
739
740!
741!-- Only check optional static driver file name, if it has been given.
742    IF (TRIM( cfg%static_driver_file ) .NE. '')  THEN
743       CALL verify_file( cfg%static_driver_file, 'static driver' )
744    ENDIF
745
746    SELECT CASE( TRIM( cfg%ic_mode ) )
747       CASE( CFG_INIT_PROFILE, CFG_INIT_VOLUME )
748       CASE DEFAULT
749          message = "Initialization mode '" // TRIM( cfg%ic_mode ) //&
750                    "' is not supported. " //&
751                    "Please select either '" // CFG_INIT_PROFILE //"' or '" // &
752                    CFG_INIT_VOLUME //"', " //&
753                    "or omit the -i/--init-mode option entirely, which corresponds "//&
754                    "to the latter."
755          CALL inifor_abort( 'validate_config', message )
756    END SELECT
757
758    SELECT CASE( TRIM( cfg%isc_mode ) )
759       CASE( CFG_INIT_SOIL_PROFILE, CFG_INIT_SOIL_VOLUME )
760       CASE DEFAULT
761          message = "Soil initialization mode '" // TRIM( cfg%isc_mode ) //&
762                    "' is not supported. " //&
763                    "Please select either '" // CFG_INIT_SOIL_PROFILE //"' or '" // &
764                    CFG_INIT_SOIL_VOLUME //"', " //&
765                    "or omit the --soil-init-mode option entirely, which corresponds "//&
766                    "to the latter."
767          CALL inifor_abort( 'validate_config', message )
768    END SELECT
769
770    SELECT CASE( TRIM(cfg%bc_mode) )
771       CASE( CFG_FORCING_HOMO, CFG_FORCING_HETERO, CFG_FORCING_NUDGING )
772       CASE DEFAULT
773          message = "Forcing mode '" // TRIM( cfg%bc_mode ) //& 
774                    "' is not supported. " //&
775                    "Please select either '" // CFG_FORCING_NUDGING //          &
776                    "', '" // CFG_FORCING_HOMO // "', or '" //                 &
777                    CFG_FORCING_HETERO // "' " //&
778                    "or omit the -f/--forcing-mode option entirely, which corresponds "//&
779                    "to the latter."
780          CALL inifor_abort( 'validate_config', message )
781    END SELECT
782
783    SELECT CASE( TRIM( cfg%averaging_mode ) )
784       CASE( 'level' )
785       CASE( 'height' )
786          message = "Averaging mode '" // TRIM( cfg%averaging_mode ) //&
787                    "' is currently not supported. " //&
788                    "Please use level-based averaging by selecting 'level', " //&
789                    "or by omitting the --averaging-mode option entirely."
790          CALL inifor_abort( 'validate_config', message )
791       CASE DEFAULT
792          message = "Averaging mode '" // TRIM( cfg%averaging_mode ) //&
793                    "' is not supported. " //&
794          !          "Please select either 'height' or 'level', " //&
795          !          "or omit the --averaging-mode option entirely, which corresponds "//&
796          !          "to the latter."
797                    "Please use level-based averaging by selecting 'level', " //&
798                    "or by omitting the --averaging-mode option entirely."
799          CALL inifor_abort( 'validate_config', message )
800    END SELECT
801
802    IF ( cfg%ug_defined_by_user .NEQV. cfg%vg_defined_by_user )  THEN
803       message = "You specified only one component of the geostrophic " // &
804                 "wind. Please specify either both or none."
805       CALL inifor_abort( 'validate_config', message )
806    ENDIF
807
808    IF ( .NOT. cfg%static_driver_is_set .AND. .NOT. cfg%z0_is_set )  THEN
809       message =                                                               &
810          "The vertical origin of the PALM grid has not been defined. " // NEW_LINE( " " ) // &
811          "Please specify the right value for your setup by either " // NEW_LINE( " " ) // &
812          "  - using the command-line option --elevation <height above sea level>, or by" // NEW_LINE( " " ) // &
813          "  - specifying a static driver file using --static <filename> in order to use " // NEW_LINE( " " ) // &
814          "    use the value of origin_z (and origin_lon and origin_lat) specifed therein."
815       CALL inifor_abort( 'validate_config', message )
816    ENDIF
817
818 END SUBROUTINE validate_config
819
820
821 SUBROUTINE get_cosmo_grid( hhl_file, soil_file, rlon, rlat, hhl, hfl, depths, &
822                            d_depth, d_depth_rho_inv, phi_n, lambda_n,       &
823                            phi_equat,                                       &
824                            lonmin_cosmo, lonmax_cosmo,                      &
825                            latmin_cosmo, latmax_cosmo,                      &
826                            nlon, nlat, nlev, ndepths )
827
828    CHARACTER(LEN=PATH), INTENT(IN)                      ::  hhl_file  !< path to file containing the HHL variable (height of half layers)
829    CHARACTER(LEN=PATH), INTENT(IN)                      ::  soil_file !< path to one of the soil input files for reading soil layer depths
830    REAL(wp), DIMENSION(:), ALLOCATABLE, INTENT(OUT)     ::  rlon      !< longitudes of COSMO-DE's rotated-pole grid
831    REAL(wp), DIMENSION(:), ALLOCATABLE, INTENT(OUT)     ::  rlat      !< latitudes of COSMO-DE's rotated-pole grid
832    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, INTENT(OUT) ::  hhl       !< heights of half layers (cell faces) above sea level in COSMO-DE, read in from external file
833    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, INTENT(OUT) ::  hfl       !< heights of full layers (cell centres) above sea level in COSMO-DE, computed as arithmetic average of hhl
834    REAL(wp), DIMENSION(:), ALLOCATABLE, INTENT(OUT)     ::  depths    !< COSMO-DE's TERRA-ML soil layer depths
835    REAL(wp), DIMENSION(:), ALLOCATABLE, INTENT(OUT)     ::  d_depth
836    REAL(wp), DIMENSION(:), ALLOCATABLE, INTENT(OUT)     ::  d_depth_rho_inv
837    REAL(wp), INTENT(OUT)                                ::  phi_n
838    REAL(wp), INTENT(OUT)                                ::  phi_equat
839    REAL(wp), INTENT(OUT)                                ::  lambda_n
840    REAL(wp), INTENT(OUT)                                ::  lonmin_cosmo !< Minimunm longitude of COSMO-DE's rotated-pole grid [COSMO rotated-pole rad]
841    REAL(wp), INTENT(OUT)                                ::  lonmax_cosmo !< Maximum longitude of COSMO-DE's rotated-pole grid [COSMO rotated-pole rad]
842    REAL(wp), INTENT(OUT)                                ::  latmin_cosmo !< Minimunm latitude of COSMO-DE's rotated-pole grid [COSMO rotated-pole rad]
843    REAL(wp), INTENT(OUT)                                ::  latmax_cosmo !< Maximum latitude of COSMO-DE's rotated-pole grid [COSMO rotated-pole rad]
844    INTEGER(iwp), INTENT(OUT)                            ::  nlon, nlat, nlev, ndepths
845
846    TYPE(nc_var) ::  cosmo_var !< COSMO dummy variable, used for reading HHL, rlon, rlat
847    INTEGER(iwp) ::  k
848
849!
850!-- Read in COSMO's heights of half layers (vertical cell faces)
851    cosmo_var%name = NC_HHL_NAME
852    CALL get_netcdf_variable( hhl_file, cosmo_var, hhl )
853    CALL get_netcdf_dim_vector( hhl_file, NC_RLON_NAME, rlon )
854    CALL get_netcdf_dim_vector( hhl_file, NC_RLAT_NAME, rlat )
855    CALL get_netcdf_dim_vector( soil_file, NC_DEPTH_NAME, depths)
856    CALL log_runtime( 'time', 'read' )
857
858    CALL reverse( hhl )
859    nlon = SIZE( hhl, 1 )
860    nlat = SIZE( hhl, 2 )
861    nlev = SIZE( hhl, 3 )
862    ndepths = SIZE( depths )
863
864    CALL log_runtime( 'time', 'comp' )
865
866    ALLOCATE( hfl( nlon, nlat, nlev-1 ) )
867    ALLOCATE( d_depth( ndepths ), d_depth_rho_inv( ndepths ) )
868    CALL log_runtime('time', 'alloc')
869
870    CALL get_soil_layer_thickness( depths, d_depth )
871    d_depth_rho_inv = 1.0_wp / ( d_depth * RHO_L )
872
873!
874!-- Compute COSMO's heights of full layers (cell centres)
875    DO  k = 1, nlev-1
876       hfl(:,:,k) = 0.5_wp * ( hhl(:,:,k) +                                  &
877                               hhl(:,:,k+1) )
878    ENDDO
879!
880!-- COSMO rotated pole coordinates
881    phi_n = TO_RADIANS                                                       &
882          * get_netcdf_variable_attribute( hhl_file,                         &
883                                           NC_ROTATED_POLE_NAME,             &
884                                           NC_POLE_LATITUDE_NAME )
885
886    lambda_n = TO_RADIANS                                                    &
887             * get_netcdf_variable_attribute( hhl_file,                      &
888                                              NC_ROTATED_POLE_NAME,          &
889                                              NC_POLE_LONGITUDE_NAME )
890
891    phi_equat = 90.0_wp * TO_RADIANS - phi_n
892
893    lonmin_cosmo = MINVAL( rlon ) * TO_RADIANS
894    lonmax_cosmo = MAXVAL( rlon ) * TO_RADIANS
895    latmin_cosmo = MINVAL( rlat ) * TO_RADIANS
896    latmax_cosmo = MAXVAL( rlat ) * TO_RADIANS
897    CALL log_runtime('time', 'comp')
898
899 END SUBROUTINE get_cosmo_grid
900
901
902!------------------------------------------------------------------------------!
903! Description:
904! ------------
905!> Fills the thickness array of the COSMO soil layers. Since COSMO's (i.e.
906!> TERRA_ML's [1]) soil layer boundaries follow the rule
907!>
908!>    depth(0) = 0.0, and
909!>    depth(k) = 0.01 * 3**(k-1), k in [1,2,3,...,7]
910!>
911!> and full levels are defined as the midpoints between two layer boundaries,
912!> all except the first layer thicknesses equal the depth of the midpoint.
913!>
914!> [1] A Description of the Nonhydrostatic Regional COSMO Model Part II :
915!>     Physical Parameterization*, Sect. 11 TERRA_ML.
916!>     http://www.cosmo-model.org/content/model/documentation/core/cosmoPhysParamtr.pdf)
917!>
918!> Input parameters:
919!> -----------------
920!>
921!> depths: array of full soil layer depths (cell midpoints)
922!>
923!>
924!> Output parameters:
925!> ------------------
926!>
927!> d_depth: array of soil layer thicknesses
928!>
929!------------------------------------------------------------------------------!
930 SUBROUTINE get_soil_layer_thickness( depths, d_depth )
931
932    REAL(wp), INTENT(IN)  ::  depths(:)
933    REAL(wp), INTENT(OUT) ::  d_depth(:)
934
935    d_depth(:) = depths(:)
936    d_depth(1) = 2.0_wp * depths(1)
937
938 END SUBROUTINE get_soil_layer_thickness
939!------------------------------------------------------------------------------!
940! Description:
941! ------------
942!> Check whether the given file is present on the filesystem.
943!------------------------------------------------------------------------------!
944 LOGICAL FUNCTION file_present( filename )
945    CHARACTER(LEN=PATH), INTENT(IN) ::  filename
946
947    INQUIRE( FILE=filename, EXIST=file_present )
948
949 END FUNCTION file_present
950
951
952!------------------------------------------------------------------------------!
953! Description:
954! ------------
955!> This routine initializes the dynamic driver file, i.e. INIFOR's netCDF output
956!> file.
957!>
958!> Besides writing metadata, such as global attributes, coordinates, variables,
959!> in the netCDF file, various netCDF IDs are saved for later, when INIFOR
960!> writes the actual data.
961!------------------------------------------------------------------------------!
962 SUBROUTINE setup_netcdf_dimensions( output_file, palm_grid,                  &
963                                     start_date_string, origin_lon, origin_lat )
964
965    TYPE(nc_file), INTENT(INOUT)      ::  output_file
966    TYPE(grid_definition), INTENT(IN) ::  palm_grid
967    CHARACTER (LEN=DATE), INTENT(IN)  ::  start_date_string
968    REAL(wp), INTENT(IN)              ::  origin_lon, origin_lat
969
970    CHARACTER (LEN=8)     ::  date_string
971    CHARACTER (LEN=10)    ::  time_string
972    CHARACTER (LEN=5)     ::  zone_string
973    CHARACTER (LEN=SNAME) ::  history_string
974    INTEGER               ::  nx, ny, nz, nt
975    INTEGER               ::  ncid, dimids(3), dimvarids(3)
976    REAL(wp)              ::  z0
977
978    message = "Initializing PALM-4U dynamic driver file '" //               &
979              TRIM(output_file%name) // "' and setting up dimensions."
980    CALL report('setup_netcdf_dimensions', message)
981
982!
983!-- Create the netCDF file as in netCDF-4/HDF5 format if __netcdf4 preprocessor flag is given
984#if defined( __netcdf4 )
985    CALL check(nf90_create(TRIM(output_file%name), OR(NF90_CLOBBER, NF90_HDF5), ncid))
986#else
987    CALL check(nf90_create(TRIM(output_file%name), NF90_CLOBBER, ncid))
988#endif
989
990!------------------------------------------------------------------------------
991!- Section 1: Define NetCDF dimensions and coordinates
992!------------------------------------------------------------------------------
993    nt = SIZE(output_file%time)
994    nx = palm_grid%nx
995    ny = palm_grid%ny
996    nz = palm_grid%nz
997    z0 = palm_grid%z0
998
999
1000!
1001!------------------------------------------------------------------------------
1002!- Section 2: Write global NetCDF attributes
1003!------------------------------------------------------------------------------
1004    CALL date_and_time(DATE=date_string, TIME=time_string, ZONE=zone_string)
1005    history_string =                                                        &
1006        'Created on '// date_string      //                                 &
1007        ' at '       // time_string(1:2) // ':' // time_string(3:4) //      &
1008        ' (UTC'      // zone_string // ')'
1009
1010    CALL check(nf90_put_att(ncid, NF90_GLOBAL, 'title',          'PALM input file for scenario ...'))
1011    CALL check(nf90_put_att(ncid, NF90_GLOBAL, 'institution',    'Deutscher Wetterdienst, Offenbach'))
1012    CALL check(nf90_put_att(ncid, NF90_GLOBAL, 'author',         'Eckhard Kadasch, eckhard.kadasch@dwd.de'))
1013    CALL check(nf90_put_att(ncid, NF90_GLOBAL, 'history',        TRIM(history_string)))
1014    CALL check(nf90_put_att(ncid, NF90_GLOBAL, 'references',     '--'))
1015    CALL check(nf90_put_att(ncid, NF90_GLOBAL, 'comment',        '--'))
1016    CALL check(nf90_put_att(ncid, NF90_GLOBAL, 'origin_lat',     TRIM(real_to_str(origin_lat*TO_DEGREES, '(F18.13)'))))
1017    CALL check(nf90_put_att(ncid, NF90_GLOBAL, 'origin_lon',     TRIM(real_to_str(origin_lon*TO_DEGREES, '(F18.13)'))))
1018!
1019!-- FIXME: This is the elevation relative to COSMO-DE/D2 sea level and does
1020!-- FIXME: not necessarily comply with DHHN2016 (c.f. PALM Input Data
1021!-- FIXME: Standard v1.9., origin_z)
1022    CALL check(nf90_put_att(ncid, NF90_GLOBAL, 'origin_z',       TRIM(real_to_str(z0, '(F18.13)'))))
1023    CALL check(nf90_put_att(ncid, NF90_GLOBAL, 'inifor_version', TRIM(VERSION)))
1024    CALL check(nf90_put_att(ncid, NF90_GLOBAL, 'palm_version',   '--'))
1025
1026!
1027!
1028!------------------------------------------------------------------------------
1029!- Section 2a: Define dimensions for cell centers (scalars in soil and atmosph.)
1030!------------------------------------------------------------------------------
1031!
1032!-- reset dimids first
1033    dimids = (/0, 0, 0/)
1034    CALL check( nf90_def_dim(ncid, "x", nx+1, dimids(1)) )
1035    CALL check( nf90_def_dim(ncid, "y", ny+1, dimids(2)) )
1036    CALL check( nf90_def_dim(ncid, "z", nz, dimids(3)) )
1037!
1038!-- save dimids for later
1039    output_file%dimids_scl = dimids 
1040
1041!
1042!-- reset dimvarids first
1043    dimvarids = (/0, 0, 0/)
1044    CALL check(nf90_def_var(ncid, "x", NF90_FLOAT, dimids(1), dimvarids(1)))
1045    CALL check(nf90_put_att(ncid, dimvarids(1), "standard_name", "x coordinate of cell centers"))
1046    CALL check(nf90_put_att(ncid, dimvarids(1), "units", "m"))
1047
1048    CALL check(nf90_def_var(ncid, "y", NF90_FLOAT, dimids(2), dimvarids(2)))
1049    CALL check(nf90_put_att(ncid, dimvarids(2), "standard_name", "y coordinate of cell centers"))
1050    CALL check(nf90_put_att(ncid, dimvarids(2), "units", "m"))
1051
1052    CALL check(nf90_def_var(ncid, "z", NF90_FLOAT, dimids(3), dimvarids(3)))
1053    CALL check(nf90_put_att(ncid, dimvarids(3), "standard_name", "z coordinate of cell centers"))
1054    CALL check(nf90_put_att(ncid, dimvarids(3), "units", "m"))
1055!
1056!-- save dimvarids for later
1057    output_file%dimvarids_scl = dimvarids
1058
1059!
1060!-- overwrite third dimid with the one of depth
1061    CALL check(nf90_def_dim(ncid, "zsoil", SIZE(palm_grid%depths), dimids(3)) )
1062!
1063!-- save dimids for later
1064    output_file%dimids_soil = dimids
1065
1066!
1067!-- overwrite third dimvarid with the one of depth
1068    CALL check(nf90_def_var(ncid, "zsoil", NF90_FLOAT, output_file%dimids_soil(3), dimvarids(3)))
1069    CALL check(nf90_put_att(ncid, dimvarids(3), "standard_name", "depth_below_land"))
1070    CALL check(nf90_put_att(ncid, dimvarids(3), "positive", "down"))
1071    CALL check(nf90_put_att(ncid, dimvarids(3), "units", "m"))
1072!
1073!-- save dimvarids for later
1074    output_file%dimvarids_soil = dimvarids
1075!
1076!------------------------------------------------------------------------------
1077!- Section 2b: Define dimensions for cell faces/velocities
1078!------------------------------------------------------------------------------
1079!
1080!-- reset dimids first
1081    dimids = (/0, 0, 0/)
1082    CALL check(nf90_def_dim(ncid, "xu", nx, dimids(1)) )
1083    CALL check(nf90_def_dim(ncid, "yv", ny, dimids(2)) )
1084    CALL check(nf90_def_dim(ncid, "zw", nz-1, dimids(3)) )
1085!
1086!-- save dimids for later
1087    output_file%dimids_vel = dimids
1088
1089!
1090!-- reset dimvarids first
1091    dimvarids = (/0, 0, 0/)
1092    CALL check(nf90_def_var(ncid, "xu", NF90_FLOAT, dimids(1), dimvarids(1)))
1093    CALL check(nf90_put_att(ncid, dimvarids(1), "standard_name", "x coordinate of cell faces"))
1094    CALL check(nf90_put_att(ncid, dimvarids(1), "units", "m"))
1095
1096    CALL check(nf90_def_var(ncid, "yv", NF90_FLOAT, dimids(2), dimvarids(2)))
1097    CALL check(nf90_put_att(ncid, dimvarids(2), "standard_name", "y coordinate of cell faces"))
1098    CALL check(nf90_put_att(ncid, dimvarids(2), "units", "m"))
1099
1100    CALL check(nf90_def_var(ncid, "zw", NF90_FLOAT, dimids(3), dimvarids(3)))
1101    CALL check(nf90_put_att(ncid, dimvarids(3), "standard_name", "z coordinate of cell faces"))
1102    CALL check(nf90_put_att(ncid, dimvarids(3), "units", "m"))
1103!
1104!-- save dimvarids for later
1105    output_file%dimvarids_vel = dimvarids
1106
1107!
1108!------------------------------------------------------------------------------
1109!- Section 2c: Define time dimension
1110!------------------------------------------------------------------------------
1111    CALL check(nf90_def_dim(ncid, "time", nt, output_file%dimid_time) )
1112    CALL check(nf90_def_var(ncid, "time", NF90_FLOAT, &
1113                                          output_file%dimid_time, &
1114                                          output_file%dimvarid_time))
1115    CALL check(nf90_put_att(ncid, output_file%dimvarid_time, "standard_name", "time"))
1116    CALL check(nf90_put_att(ncid, output_file%dimvarid_time, "long_name", "time"))
1117    CALL check(nf90_put_att(ncid, output_file%dimvarid_time, "units",     &
1118                            "seconds since " // start_date_string // " UTC"))
1119
1120    CALL check(nf90_enddef(ncid))
1121
1122!
1123!------------------------------------------------------------------------------
1124!- Section 3: Write grid coordinates
1125!------------------------------------------------------------------------------
1126    CALL check(nf90_put_var(ncid, output_file%dimvarids_scl(1), palm_grid%x))
1127    CALL check(nf90_put_var(ncid, output_file%dimvarids_scl(2), palm_grid%y))
1128    CALL check(nf90_put_var(ncid, output_file%dimvarids_scl(3), palm_grid%z))
1129
1130    CALL check(nf90_put_var(ncid, output_file%dimvarids_vel(1), palm_grid%xu))
1131    CALL check(nf90_put_var(ncid, output_file%dimvarids_vel(2), palm_grid%yv))
1132    CALL check(nf90_put_var(ncid, output_file%dimvarids_vel(3), palm_grid%zw))
1133
1134!
1135!-- TODO Read in soil depths from input file before this.
1136    CALL check(nf90_put_var(ncid, output_file%dimvarids_soil(3), palm_grid%depths))
1137
1138!
1139!-- Write time vector
1140    CALL check(nf90_put_var(ncid, output_file%dimvarid_time, output_file%time))
1141
1142!
1143!-- Close the file
1144    CALL check(nf90_close(ncid))
1145
1146 END SUBROUTINE setup_netcdf_dimensions
1147
1148
1149!------------------------------------------------------------------------------!
1150! Description:
1151! ------------
1152!> Defines the netCDF variables to be written to the dynamic driver file
1153!------------------------------------------------------------------------------!
1154 SUBROUTINE setup_netcdf_variables(filename, io_group_list)
1155
1156    CHARACTER (LEN=*), INTENT(IN)       ::  filename
1157    TYPE(io_group), INTENT(IN), TARGET  ::  io_group_list(:)
1158
1159    TYPE(io_group), POINTER             ::  group
1160    TYPE(nc_var), POINTER               ::  var
1161    INTEGER(iwp)                        ::  group_idx, var_idx, n_var
1162    INTEGER                             ::  ncid
1163    LOGICAL                             ::  to_be_written
1164
1165    message = "Defining variables in dynamic driver '" // TRIM(filename) // "'."
1166    CALL report('setup_netcdf_variables', message)
1167
1168    CALL check(nf90_open(TRIM(filename), NF90_WRITE, ncid))
1169    CALL check(nf90_redef(ncid))
1170
1171    n_var = 0
1172    DO  group_idx = 1, SIZE( io_group_list )
1173
1174       group => io_group_list(group_idx)
1175       DO var_idx = 1, SIZE( group%out_vars )
1176
1177          to_be_written = .FALSE.
1178
1179          IF (ALLOCATED( group%out_vars ))  THEN
1180             var => group%out_vars(var_idx)
1181             n_var = n_var + 1
1182
1183             to_be_written = (                                                 &
1184                group%to_be_processed  .AND.  var%to_be_processed              &
1185                .AND. .NOT. var%is_internal                                    &
1186             )
1187          ENDIF
1188
1189          IF ( to_be_written )  THEN
1190             message = "  variable #" // TRIM(str(n_var)) // " '" // TRIM(var%name) // "'."
1191             CALL report('setup_netcdf_variables', message)
1192
1193             CALL netcdf_define_variable(var, ncid)
1194             CALL netcdf_get_dimensions(var, ncid)
1195          ENDIF
1196
1197       ENDDO
1198       
1199    ENDDO
1200
1201    CALL check(nf90_enddef(ncid))
1202    CALL check(nf90_close(ncid))
1203
1204    message = "Dynamic driver '" // TRIM(filename) // "' initialized successfully."
1205    CALL report('setup_netcdf_variables', message)
1206
1207 END SUBROUTINE setup_netcdf_variables
1208
1209
1210!------------------------------------------------------------------------------!
1211! Description:
1212! ------------
1213!> This routine reads and returns all input variables of the given IO group
1214!> It accomodates the data by allocating a container variable that represents a
1215!> list of arrays of the same length as the groups variable list. (This list
1216!> will typically contain one or two items.) After the container, its members
1217!> are allocated one by one with the appropriate, possibly different,
1218!> dimensions.
1219!>
1220!> The 'group' is an INTENT(INOUT) variable so that 'get_netcdf_variable()' can
1221!> record netCDF IDs in the 'in_var_list()' member variable.
1222!------------------------------------------------------------------------------!
1223 SUBROUTINE read_input_variables(group, iter, buffer)
1224    TYPE(io_group), INTENT(INOUT), TARGET       ::  group
1225    INTEGER(iwp), INTENT(IN)                    ::  iter
1226    TYPE(container), ALLOCATABLE, INTENT(INOUT) ::  buffer(:)
1227    INTEGER(iwp)                                ::  hour, buf_id
1228    TYPE(nc_var), POINTER                       ::  input_var
1229    CHARACTER(LEN=PATH), POINTER                ::  input_file
1230    INTEGER(iwp)                                ::  ivar, nbuffers
1231
1232    message = "Reading data for I/O group '" // TRIM(group%in_var_list(1)%name) // "'."
1233    CALL report('read_input_variables', message)
1234
1235    input_file => group%in_files(iter)
1236
1237!
1238!------------------------------------------------------------------------------
1239!- Section 1: Load input buffers for accumulated variables
1240!------------------------------------------------------------------------------
1241!
1242!-- radiation budgets, precipitation
1243    IF ( group%kind == 'running average' .OR.                                  &
1244         group%kind == 'accumulated' )  THEN
1245
1246       IF ( SIZE( group%in_var_list ) .GT. 1 ) THEN
1247          message = "I/O groups may not contain more than one " // & 
1248                    "accumulated variable. Group '" // TRIM(group%kind) //&
1249                    "' contains " //                                           &
1250                    TRIM( str( SIZE( group%in_var_list, kind=iwp ) ) ) // "."
1251          CALL inifor_abort( 'read_input_variables | accumulation', message )
1252       ENDIF
1253
1254!
1255!--    use two buffer arrays
1256       nbuffers = 2
1257       IF ( .NOT. ALLOCATED( buffer ) )  ALLOCATE( buffer(nbuffers) )
1258
1259!
1260!--    hour of the day
1261       hour = iter - 1
1262!
1263!--    chose correct buffer array
1264       buf_id = select_buffer(hour)
1265
1266       CALL log_runtime('time', 'read')
1267       IF ( ALLOCATED(buffer(buf_id)%array) )  DEALLOCATE(buffer(buf_id)%array)
1268       CALL log_runtime('time', 'alloc')
1269
1270       input_var => group%in_var_list(1)
1271       CALL get_netcdf_variable(input_file, input_var, buffer(buf_id)%array)
1272       CALL report('read_input_variables', "Read accumulated " // TRIM(group%in_var_list(1)%name)) 
1273
1274       IF ( input_var%is_upside_down )  CALL reverse(buffer(buf_id)%array)
1275       CALL log_runtime('time', 'comp')
1276         
1277!------------------------------------------------------------------------------
1278!- Section 2: Load input buffers for normal I/O groups
1279!------------------------------------------------------------------------------
1280    ELSE
1281
1282!
1283!--    Allocate one input buffer per output quantity. If more quantities
1284!--    have to be computed than input variables exist in this group,
1285!--    allocate more buffers. For instance, in the thermodynamics group,
1286!--    there are three input variabels (PP, T, Qv) and four quantities
1287!--    necessart (P, Theta, Rho, qv) for the corresponding output fields
1288!--    (p0, Theta, qv, ug, and vg)
1289       ALLOCATE( buffer(group%n_output_quantities) )
1290       CALL log_runtime('time', 'alloc')
1291       
1292!
1293!--    Read in all input variables, leave extra buffers-if any-untouched.
1294       DO  ivar = 1, group%n_inputs
1295
1296          input_var => group%in_var_list(ivar)
1297
1298          IF ( input_var%to_be_processed )  THEN
1299!            Check wheather P or PP is present in input file
1300             IF (input_var%name == 'P')  THEN
1301                input_var%name = TRIM( get_pressure_varname(input_file) )
1302             CALL log_runtime('time', 'read')
1303             ENDIF
1304
1305             CALL get_netcdf_variable(input_file, input_var, buffer(ivar)%array)
1306
1307             IF ( input_var%is_upside_down )  CALL reverse(buffer(ivar)%array)
1308             CALL log_runtime('time', 'comp')
1309          ENDIF
1310
1311       ENDDO
1312    ENDIF
1313
1314 END SUBROUTINE read_input_variables
1315
1316
1317!------------------------------------------------------------------------------!
1318! Description:
1319! ------------
1320!> Select the appropriate buffer ID for accumulated COSMO input variables
1321!> depending on the current hour.
1322!------------------------------------------------------------------------------!
1323 INTEGER(iwp) FUNCTION select_buffer(hour)
1324    INTEGER(iwp), INTENT(IN) ::  hour
1325    INTEGER(iwp)             ::  step
1326
1327    select_buffer = 0_iwp
1328    step = MODULO(hour, 3_iwp) + 1_iwp
1329
1330    SELECT CASE(step)
1331       CASE(1, 3)
1332           select_buffer = 1_iwp
1333       CASE(2)
1334           select_buffer = 2_iwp
1335       CASE DEFAULT
1336           message = "Invalid step '" // TRIM(str(step))
1337           CALL inifor_abort('select_buffer', message)
1338    END SELECT
1339 END FUNCTION select_buffer
1340
1341
1342!------------------------------------------------------------------------------!
1343! Description:
1344! ------------
1345!> Checks if the input_file contains the absolute pressure, 'P', or the pressure
1346!> perturbation, 'PP', and returns the appropriate string.
1347!------------------------------------------------------------------------------!
1348 CHARACTER(LEN=2) FUNCTION get_pressure_varname(input_file) RESULT(var)
1349    CHARACTER(LEN=*) ::  input_file
1350    INTEGER          ::  ncid, varid
1351
1352    CALL check(nf90_open( TRIM(input_file), NF90_NOWRITE, ncid ))
1353    IF ( nf90_inq_varid( ncid, 'P', varid ) .EQ. NF90_NOERR )  THEN
1354
1355       var = 'P'
1356
1357    ELSE IF ( nf90_inq_varid( ncid, 'PP', varid ) .EQ. NF90_NOERR )  THEN
1358
1359       var = 'PP'
1360       CALL report('get_pressure_var', 'Using PP instead of P')
1361
1362    ELSE
1363
1364       message = "Failed to read '" // TRIM(var) // &
1365                 "' from file '" // TRIM(input_file) // "'."
1366       CALL inifor_abort('get_pressure_var', message)
1367
1368    ENDIF
1369
1370    CALL check(nf90_close(ncid))
1371
1372 END FUNCTION get_pressure_varname
1373
1374
1375!------------------------------------------------------------------------------!
1376! Description:
1377! ------------
1378!> Read the given global attribute form the given netCDF file.
1379!------------------------------------------------------------------------------!
1380 FUNCTION get_netcdf_attribute(filename, attribute) RESULT(attribute_value)
1381
1382    CHARACTER(LEN=*), INTENT(IN) ::  filename, attribute
1383    REAL(wp)                     ::  attribute_value
1384
1385    INTEGER                      ::  ncid
1386
1387    IF ( nf90_open( TRIM(filename), NF90_NOWRITE, ncid ) == NF90_NOERR )  THEN
1388
1389       CALL check(nf90_get_att(ncid, NF90_GLOBAL, TRIM(attribute), attribute_value))
1390       CALL check(nf90_close(ncid))
1391
1392    ELSE
1393
1394       message = "Failed to read '" // TRIM(attribute) // &
1395                 "' from file '" // TRIM(filename) // "'."
1396       CALL inifor_abort('get_netcdf_attribute', message)
1397
1398    ENDIF
1399
1400 END FUNCTION get_netcdf_attribute
1401
1402
1403!------------------------------------------------------------------------------!
1404! Description:
1405! ------------
1406!> Read the attribute of the given variable form the given netCDF file.
1407!------------------------------------------------------------------------------!
1408 FUNCTION get_netcdf_variable_attribute(filename, varname, attribute)       &
1409    RESULT(attribute_value)
1410
1411    CHARACTER(LEN=*), INTENT(IN) ::  filename, varname, attribute
1412    REAL(wp)                     ::  attribute_value
1413
1414    INTEGER                      ::  ncid, varid
1415
1416    IF ( nf90_open( TRIM(filename), NF90_NOWRITE, ncid ) == NF90_NOERR )  THEN
1417
1418       CALL check( nf90_inq_varid( ncid, TRIM( varname ), varid ) )
1419       CALL check( nf90_get_att( ncid, varid, TRIM( attribute ),            &
1420                   attribute_value ) )
1421       CALL check( nf90_close( ncid ) )
1422
1423    ELSE
1424
1425       message = "Failed to read '" // TRIM( varname ) // ":" //            &
1426                 TRIM( attribute ) // "' from file '" // TRIM(filename) // "'."
1427       CALL inifor_abort('get_netcdf_variable_attribute', message)
1428
1429    ENDIF
1430
1431 END FUNCTION get_netcdf_variable_attribute
1432
1433!------------------------------------------------------------------------------!
1434! Description:
1435! ------------
1436!> Updates the dynamic driver with the interpolated field of the current
1437!> variable at the current time step.
1438!------------------------------------------------------------------------------!
1439 SUBROUTINE update_output(var, array, iter, output_file, cfg)
1440    TYPE(nc_var), INTENT(IN)  ::  var
1441    REAL(wp), INTENT(IN)      ::  array(:,:,:)
1442    INTEGER(iwp), INTENT(IN)  ::  iter
1443    TYPE(nc_file), INTENT(IN) ::  output_file
1444    TYPE(inifor_config)       ::  cfg
1445
1446    INTEGER      ::  ncid, ndim, start(4), count(4)
1447    LOGICAL      ::  var_is_time_dependent
1448
1449    var_is_time_dependent = (                                                  &
1450       var%dimids( var%ndim ) == output_file%dimid_time                        &
1451    )
1452
1453!
1454!-- Skip time dimension for output
1455    ndim = var%ndim
1456    IF ( var_is_time_dependent )  ndim = var%ndim - 1
1457
1458    start(:)      = (/1,1,1,1/)
1459    start(ndim+1) = iter
1460    count(1:ndim) = var%dimlen(1:ndim)
1461
1462    CALL check(nf90_open(output_file%name, NF90_WRITE, ncid))
1463
1464!
1465!-- Reduce dimension of output array according to variable kind
1466    SELECT CASE (TRIM(var%kind))
1467       
1468       CASE( 'init scalar profile', 'init u profile', 'init v profile',        &
1469             'init w profile', 'init soil profile' )
1470
1471          CALL check(nf90_put_var( ncid, var%varid, array(1,1,:) ) )
1472
1473       CASE( 'init soil', 'init scalar', 'init u', 'init v', 'init w' )
1474
1475          CALL check(nf90_put_var( ncid, var%varid, array(:,:,:) ) )
1476
1477       CASE( 'left scalar', 'right scalar', 'left w', 'right w' ) 
1478
1479          CALL check(nf90_put_var( ncid, var%varid, array(1,:,:),              &
1480                                   start=start(1:ndim+1),                      &
1481                                   count=count(1:ndim) ) )
1482         
1483
1484          IF (.NOT. SIZE(array, 2) .EQ. var%dimlen(1))  THEN
1485             PRINT *, "inifor: update_output: Dimension ", 1, " of variable ", &
1486                 TRIM(var%name), " (", var%dimlen(1),                          &
1487                 ") does not match the dimension of the output array (",       &
1488                 SIZE(array, 2), ")."
1489             STOP
1490          ENDIF
1491         
1492
1493       CASE( 'north scalar', 'south scalar', 'north w', 'south w' )
1494
1495          CALL check(nf90_put_var( ncid, var%varid, array(:,1,:),              &
1496                                   start=start(1:ndim+1),                      &
1497                                   count=count(1:ndim) ) )
1498         
1499
1500       CASE( 'surface forcing', 'top scalar', 'top w' )
1501
1502          CALL check(nf90_put_var( ncid, var%varid, array(:,:,1),              &
1503                                   start=start(1:ndim+1),                      &
1504                                   count=count(1:ndim) ) )
1505         
1506       CASE( 'left u', 'right u', 'left v', 'right v' )
1507
1508          CALL check(nf90_put_var( ncid, var%varid, array(1,:,:),              &
1509                                   start=start(1:ndim+1),                      &
1510                                   count=count(1:ndim) ) )
1511
1512       CASE( 'north u', 'south u', 'north v', 'south v' )
1513
1514          CALL check(nf90_put_var( ncid, var%varid, array(:,1,:),              &
1515                                   start=start(1:ndim+1),                      &
1516                                   count=count(1:ndim) ) )
1517
1518       CASE( 'top u', 'top v' )
1519
1520          CALL check(nf90_put_var( ncid, var%varid, array(:,:,1),              &
1521                                   start=start(1:ndim+1),                      &
1522                                   count=count(1:ndim) ) )
1523
1524       CASE( 'time series' )
1525
1526          CALL check(nf90_put_var( ncid, var%varid, array(1,1,1),              &
1527                                   start=start(1:ndim+1) ) )
1528
1529       CASE( 'constant scalar profile', 'geostrophic',                         &
1530             'left scalar profile', 'right scalar profile',                    &
1531             'north scalar profile', 'south scalar profile',                   &
1532             'left u profile', 'right u profile',                              &
1533             'north u profile', 'south u profile',                             &
1534             'left v profile', 'right v profile',                              &
1535             'north v profile', 'south v profile',                             &
1536             'top scalar profile', 'top u profile', 'top v profile',           &
1537             'top w profile',                                                  &
1538             'left w profile', 'right w profile',                              &
1539             'north w profile', 'south w profile' )
1540
1541          CALL check(nf90_put_var( ncid, var%varid, array(1,1,:),              &
1542                                   start=start(1:ndim+1),                      &
1543                                   count=count(1:ndim) ) )
1544
1545       CASE ( 'internal profile' )
1546
1547          IF ( cfg%debug )  THEN
1548             CALL check(nf90_put_var( ncid, var%varid, array(1,1,:),           &
1549                                      start=start(1:ndim+1),                   &
1550                                      count=count(1:ndim) ) )
1551          ENDIF
1552
1553       CASE ( 'large-scale scalar forcing', 'large-scale w forcing' )
1554
1555           message = "Doing nothing in terms of writing large-scale forings."
1556           CALL report('update_output', message)
1557
1558       CASE DEFAULT
1559
1560           message = "Variable kind '" // TRIM(var%kind) //                  &
1561                    "' not recognized."
1562           CALL inifor_abort('update_output', message)
1563
1564    END SELECT
1565
1566    CALL check(nf90_close(ncid))
1567
1568 END SUBROUTINE update_output
1569
1570
1571!------------------------------------------------------------------------------!
1572! Description:
1573! ------------
1574!> Checks the status of a netCDF API call and aborts if an error occured
1575!------------------------------------------------------------------------------!
1576 SUBROUTINE check(status)
1577
1578    INTEGER, INTENT(IN) ::  status
1579
1580    IF (status /= nf90_noerr)  THEN
1581       message = "NetCDF API call failed with error: " //                     &
1582                 TRIM( nf90_strerror(status) )
1583       CALL inifor_abort('io.check', message)
1584    ENDIF
1585
1586 END SUBROUTINE check
1587
1588
1589!------------------------------------------------------------------------------!
1590! Description:
1591! ------------
1592!> Setup the origin of the PALM coordinate system based on what is given in the
1593!> INIFOR namelist file and via an optional static driver.
1594!------------------------------------------------------------------------------!
1595 SUBROUTINE set_palm_origin( cfg, namelist_longitude, namelist_latitude,       &
1596                             origin_lon, origin_lat, z0 )
1597
1598    TYPE(inifor_config), INTENT(IN) ::  cfg
1599    REAL(wp), INTENT(IN)            ::  namelist_longitude, namelist_latitude
1600    REAL(wp), INTENT(OUT)           ::  origin_lon, origin_lat, z0
1601
1602    message = 'Reading PALM-4U origin from'
1603    IF ( TRIM( cfg%static_driver_file ) .NE. '' )  THEN
1604
1605       origin_lon = get_netcdf_attribute( cfg%static_driver_file, TRIM( PIDS_ORIGIN_LON ) )
1606       origin_lat = get_netcdf_attribute( cfg%static_driver_file, TRIM( PIDS_ORIGIN_LAT ) )
1607       z0         = get_netcdf_attribute( cfg%static_driver_file, TRIM( PIDS_ORIGIN_Z ) )
1608
1609       message = TRIM(message) // " static driver file '"                      &
1610                               // TRIM( cfg%static_driver_file ) // "'"
1611
1612
1613    ELSE
1614
1615       origin_lon = namelist_longitude
1616       origin_lat = namelist_latitude
1617
1618       message = TRIM( message ) // " namlist file '"                          &
1619                                 // TRIM( cfg%namelist_file ) // "'"
1620
1621    ENDIF
1622    origin_lon = origin_lon * TO_RADIANS
1623    origin_lat = origin_lat * TO_RADIANS
1624
1625    CALL report('set_palm_origin', message)
1626
1627    IF ( cfg%z0_is_set )  THEN
1628       z0 = cfg%z0
1629       IF ( TRIM( cfg%static_driver_file ) .NE. '' )  THEN
1630          message = 'You specified both --static-driver/-t and --elevation/-z0. ' // &
1631                    'Using the command line value (' // TRIM( real_to_str_f( cfg%z0 ) ) // &
1632                    ') instead of static driver value (' // TRIM( real_to_str_f( z0 ) ) // ').'
1633          CALL warn( 'set_palm_origin', message )
1634       ENDIF
1635    ENDIF
1636
1637 END SUBROUTINE set_palm_origin
1638
1639
1640!------------------------------------------------------------------------------!
1641! Description:
1642! ------------
1643! This function is meant to check weather a COSMO soil variable has an
1644! additional and redunant surface value at depth = 0.0. For instance operational
1645! DWD COSMO output contains the surface temperature in T_SO as a copy of the
1646! values in the first soil layer.
1647!------------------------------------------------------------------------------!
1648 LOGICAL FUNCTION has_surface_value( soil_var, filename ) 
1649
1650    TYPE(nc_var), INTENT(IN)     ::  soil_var
1651    CHARACTER(LEN=*), INTENT(IN) ::  filename
1652   
1653    REAL(wp), ALLOCATABLE        ::  depths(:)
1654
1655    CALL get_dimension_vector_of_variable(                                     &
1656       soil_var%name,                                                          &
1657       dim_idx = NC_DEPTH_DIM_IDX,                                             &
1658       filename = filename,                                                    &
1659       dim_vector = depths                                                     &
1660    )
1661
1662    has_surface_value = nearly_equal( depths(1), 0.0_wp, 10 * EPSILON(1.0_wp) )
1663
1664 END FUNCTION has_surface_value 
1665
1666
1667!------------------------------------------------------------------------------!
1668! Description:
1669! ------------
1670! This routine reads the dim_idx-th dimension vector of the variable varname
1671! from netCDF file filename. It is used for finding the depth coordinate vector
1672! of COSMO soil variables without knowing its name.
1673!------------------------------------------------------------------------------!
1674 SUBROUTINE get_dimension_vector_of_variable( varname, dim_idx, filename, dim_vector )
1675    CHARACTER(LEN=*), INTENT(IN)              ::  varname, filename
1676    INTEGER, INTENT(IN)                       ::  dim_idx
1677
1678    REAL(wp), INTENT(OUT), ALLOCATABLE ::  dim_vector(:)
1679
1680    INTEGER                            ::  dimids(NF90_MAX_VAR_DIMS)
1681    INTEGER                            ::  varid
1682    CHARACTER(LEN=NF90_MAX_NAME)       ::  dimname
1683
1684    INTEGER ::  ncid
1685
1686    IF ( nf90_open( TRIM( filename ), NF90_NOWRITE, ncid ) .EQ. NF90_NOERR )  THEN
1687
1688!
1689!--    get id of variable varname
1690       CALL check( nf90_inq_varid( ncid, TRIM( varname ), varid ) )
1691       
1692!
1693!--    get dimension ids of variable with varid
1694       CALL check( nf90_inquire_variable( ncid, varid, dimids = dimids ) )
1695
1696!
1697!--    get name of dim_idx-th dimension variable
1698       CALL check( nf90_inquire_dimension( ncid, dimids(dim_idx), name = dimname ) )
1699       CALL check( nf90_close( ncid ) )
1700
1701    ELSE
1702
1703       message = "Failed to open file '" // TRIM(filename) // "'."
1704       CALL inifor_abort('get_netcdf_variable', message)
1705
1706    ENDIF
1707
1708    ! get dimension vector with dimname
1709    CALL get_netcdf_dim_vector( filename, dimname, dim_vector )
1710
1711 END SUBROUTINE get_dimension_vector_of_variable
1712
1713
1714 LOGICAL FUNCTION netcdf_variable_present_in_file( varname, filename )
1715    CHARACTER(LEN=SNAME), INTENT(IN) ::  varname
1716    CHARACTER(LEN=PATH), INTENT(IN)  ::  filename
1717
1718    INTEGER ::  ncid, varid
1719   
1720    netcdf_variable_present_in_file = (                                        &
1721       nf90_open( TRIM( filename ), NF90_NOWRITE, ncid ) .EQ. NF90_NOERR .AND. &
1722       nf90_inq_varid( ncid, varname, varid ) .EQ. NF90_NOERR .AND.            &
1723       nf90_close( ncid ) .EQ. NF90_NOERR                                      &
1724    )
1725 
1726 END FUNCTION netcdf_variable_present_in_file
1727
1728
1729 END MODULE inifor_io
1730#endif
Note: See TracBrowser for help on using the repository browser.