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

Last change on this file since 3737 was 3716, checked in by eckhard, 6 years ago

inifor: bugfix: removed dependency on soilmoisture input files; added netcdf preprocessor flag

  • Property svn:keywords set to Id
File size: 49.8 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-2018 Leibniz Universitaet Hannover
18! Copyright 2017-2018 Deutscher Wetterdienst Offenbach
19!------------------------------------------------------------------------------!
20!
21! Current revisions:
22! -----------------
23!
24!
25! Former revisions:
26! -----------------
27! $Id: inifor_io.f90 3716 2019-02-05 17:02:38Z suehring $
28! Removed dependency on soilmoisture input files
29!
30!
31! 3680 2019-01-18 14:54:12Z knoop
32! Moved get_input_file_list() here from grid module, added check for presence of
33!    input files
34!
35!
36!
37!
38! 3618 2018-12-10 13:25:22Z eckhard
39! Prefixed all INIFOR modules with inifor_
40!
41!
42! 3615 2018-12-10 07:21:03Z raasch
43! bugfix: abort replaced by inifor_abort
44!
45! 3557 2018-11-22 16:01:22Z eckhard
46! Updated documentation, removed unused subroutine write_netcdf_variable_2d()
47!
48!
49! 3537 2018-11-20 10:53:14Z eckhard
50! New routine get_netcdf_dim_vector()
51!
52!
53! 3534 2018-11-19 15:35:16Z raasch
54! bugfix: INTENT attribute changed
55!
56! 3456 2018-10-30 14:29:54Z eckhard
57! NetCDf output of internal arrays only with --debug option
58!
59!
60! 3447 2018-10-29 15:52:54Z eckhard
61! Removed INCLUDE statement for get_netcdf_variable()
62! Renamed source files for compatibilty with PALM build system
63!
64!
65! 3395 2018-10-22 17:32:49Z eckhard
66! Added command-line options for configuring the computation of geostrophic
67!     winds (--averagin-mode, --averaging-angle)
68! Added command-line option --input-prefix for setting input file prefixes all
69!     at once
70! Added --debug option for more verbose terminal output
71! Added option-specific *_is_set LOGICALs to indicate invocation from the
72!     command-line
73! Improved error messages in case of empty file-name strings
74! Improved routine naming
75!
76! 3183 2018-07-27 14:25:55Z suehring
77! Introduced new PALM grid stretching
78! Updated variable names and metadata for PIDS v1.9 compatibility
79! Improved handling of the start date string
80! Better compatibility with older Intel compilers:
81! - avoiding implicit array allocation with new get_netcdf_variable()
82!   subroutine instead of function
83! Improved command line interface:
84! - Added configuration validation
85! - New options to configure input file prefixes
86! - GNU-style short and long option names
87! - Added version and copyright output
88!
89!
90! 3182 2018-07-27 13:36:03Z suehring
91! Initial revision
92!
93!
94!
95! Authors:
96! --------
97!> @author Eckhard Kadasch (Deutscher Wetterdienst, Offenbach)
98!
99! Description:
100! ------------
101!> The io module contains the functions needed to read and write netCDF data in
102!> INIFOR.
103!------------------------------------------------------------------------------!
104#if defined ( __netcdf )
105 MODULE inifor_io
106
107    USE inifor_control
108    USE inifor_defs,                                                           &
109        ONLY:  DATE, SNAME, PATH, PI, dp, hp, TO_RADIANS, TO_DEGREES, VERSION
110    USE inifor_types
111    USE inifor_util,                                                           &
112        ONLY:  add_hours_to, reverse, str, real_to_str
113    USE netcdf
114
115    IMPLICIT NONE
116
117!------------------------------------------------------------------------------!
118! Description:
119! ------------
120!> get_netcdf_variable() reads the netCDF data and metadate for the netCDF
121!> variable 'in_var % name' from the file 'in_file'. The netCDF data array is
122!> stored in the 'buffer' array and metadata added to the respective members of
123!> the given 'in_var'.
124!------------------------------------------------------------------------------!
125    INTERFACE get_netcdf_variable
126        MODULE PROCEDURE get_netcdf_variable_int
127        MODULE PROCEDURE get_netcdf_variable_real
128    END INTERFACE get_netcdf_variable
129
130    PRIVATE ::  get_netcdf_variable_int, get_netcdf_variable_real
131
132 CONTAINS
133
134!------------------------------------------------------------------------------!
135! Description:
136! ------------
137!> get_netcdf_variable_int() implements the integer variant for the
138!> get_netcdf_variable interface.
139!------------------------------------------------------------------------------!
140    SUBROUTINE get_netcdf_variable_int(in_file, in_var, buffer)
141
142       CHARACTER(LEN=PATH), INTENT(IN)         ::  in_file
143       TYPE(nc_var), INTENT(INOUT)             ::  in_var
144       INTEGER(hp), ALLOCATABLE, INTENT(INOUT) ::  buffer(:,:,:)
145
146       INTEGER               ::  ncid
147       INTEGER, DIMENSION(3) ::  start, count
148
149       IF ( nf90_open( TRIM(in_file), NF90_NOWRITE, ncid ) .EQ. NF90_NOERR .AND. &
150            nf90_inq_varid( ncid, in_var % name, in_var % varid ) .EQ. NF90_NOERR )  THEN
151
152          CALL get_input_dimensions(in_var, ncid)
153
154          CALL get_netcdf_start_and_count(in_var, start, count)
155 CALL run_control('time', 'read')
156
157          ALLOCATE( buffer( count(1), count(2), count(3) ) )
158 CALL run_control('time', 'alloc')
159
160          CALL check(nf90_get_var( ncid, in_var % varid, buffer,                  &
161                                   start = start,                                 &
162                                   count = count ))
163
164       ELSE
165
166          message = "Failed to read '" // TRIM(in_var % name) // &
167             "' from file '" // TRIM(in_file) // "'."
168          CALL inifor_abort('get_netcdf_variable', message)
169
170       END IF
171
172       CALL check(nf90_close(ncid))
173 CALL run_control('time', 'read')
174
175    END SUBROUTINE get_netcdf_variable_int
176
177
178!------------------------------------------------------------------------------!
179! Description:
180! ------------
181!> get_netcdf_variable_real() implements the real variant for the
182!> get_netcdf_variable interface.
183!------------------------------------------------------------------------------!
184    SUBROUTINE get_netcdf_variable_real(in_file, in_var, buffer)
185
186       CHARACTER(LEN=PATH), INTENT(IN)      ::  in_file
187       TYPE(nc_var), INTENT(INOUT)          ::  in_var
188       REAL(dp), ALLOCATABLE, INTENT(INOUT) ::  buffer(:,:,:)
189
190       INTEGER               ::  ncid
191       INTEGER, DIMENSION(3) ::  start, count
192
193       IF ( nf90_open( TRIM(in_file), NF90_NOWRITE, ncid ) .EQ. NF90_NOERR .AND. &
194            nf90_inq_varid( ncid, in_var % name, in_var % varid ) .EQ. NF90_NOERR )  THEN
195
196          CALL get_input_dimensions(in_var, ncid)
197
198          CALL get_netcdf_start_and_count(in_var, start, count)
199 CALL run_control('time', 'read')
200
201          ALLOCATE( buffer( count(1), count(2), count(3) ) )
202 CALL run_control('time', 'alloc')
203
204          CALL check(nf90_get_var( ncid, in_var % varid, buffer,                  &
205                                   start = start,                                 &
206                                   count = count ))
207
208       ELSE
209
210          message = "Failed to read '" // TRIM(in_var % name) // &
211             "' from file '" // TRIM(in_file) // "'."
212          CALL inifor_abort('get_netcdf_variable', message)
213
214       END IF
215
216       CALL check(nf90_close(ncid))
217 CALL run_control('time', 'read')
218
219    END SUBROUTINE get_netcdf_variable_real
220
221
222!------------------------------------------------------------------------------!
223! Description:
224! ------------
225!> get_netcdf_dim_vector() reads the coordinate array 'coordname' from the
226!> netCDF file 'filename'.
227!------------------------------------------------------------------------------!
228    SUBROUTINE get_netcdf_dim_vector(filename, coordname, coords)
229
230       CHARACTER(LEN=*), INTENT(IN)         ::  filename
231       CHARACTER(LEN=*), INTENT(IN)         ::  coordname
232       REAL(dp), ALLOCATABLE, INTENT(INOUT) ::  coords(:)
233
234       INTEGER ::  ncid, varid, dimlen
235       INTEGER ::  dimids(NF90_MAX_VAR_DIMS)
236
237       IF ( nf90_open( TRIM(filename), NF90_NOWRITE, ncid ) .EQ. NF90_NOERR .AND. &
238            nf90_inq_varid( ncid, coordname, varid ) .EQ. NF90_NOERR )  THEN
239
240          CALL check(nf90_inquire_variable( ncid, varid, dimids = dimids ))
241          CALL check(nf90_inquire_dimension( ncid, dimids(1), len = dimlen ))
242
243          ALLOCATE(coords(dimlen))
244          CALL check(nf90_get_var( ncid, varid, coords))
245
246       ELSE
247
248          message = "Failed to read '" // TRIM(coordname) // &
249             "' from file '" // TRIM(filename) // "'."
250          CALL inifor_abort('get_netcdf_dim_vector', message)
251
252       END IF
253
254    END SUBROUTINE get_netcdf_dim_vector
255
256
257!------------------------------------------------------------------------------!
258! Description:
259! ------------
260!> get_input_dimensions() reads dimensions metadata of the netCDF variable given
261!> by 'in_var % name' into 'in_var' data structure.
262!------------------------------------------------------------------------------!
263    SUBROUTINE get_input_dimensions(in_var, ncid)
264
265       TYPE(nc_var), INTENT(INOUT) ::  in_var
266       INTEGER, INTENT(IN)         ::  ncid
267
268       INTEGER ::  i
269
270       CALL check(nf90_get_att( ncid, in_var % varid, "long_name",             &
271                                in_var % long_name))
272
273       CALL check(nf90_get_att( ncid, in_var % varid, "units",                 &
274                                in_var % units ))
275
276       CALL check(nf90_inquire_variable( ncid, in_var % varid,                 &
277                                         ndims  = in_var % ndim,               &
278                                         dimids = in_var % dimids ))
279
280       DO i = 1, in_var % ndim
281          CALL check(nf90_inquire_dimension( ncid, in_var % dimids(i),         &
282                                             name = in_var % dimname(i),       &
283                                             len  = in_var % dimlen(i) ))
284       END DO
285
286    END SUBROUTINE get_input_dimensions
287
288
289!------------------------------------------------------------------------------!
290! Description:
291! ------------
292!> get_netcdf_start_and_count() gets the start position and element counts for
293!> the given netCDF file. This information is used in get_netcdf_variable_int()
294!> and _real() for reading input variables..
295!------------------------------------------------------------------------------!
296    SUBROUTINE get_netcdf_start_and_count(in_var, start, count)
297
298       TYPE(nc_var), INTENT(INOUT)        ::  in_var
299       INTEGER, DIMENSION(3), INTENT(OUT) ::  start, count
300
301       INTEGER ::  ndim
302
303       IF ( in_var % ndim .LT. 2  .OR.  in_var % ndim .GT. 4 )  THEN
304
305          message = "Failed reading NetCDF variable " //                       &
306             TRIM(in_var % name) // " with " // TRIM(str(in_var % ndim)) //    &
307             " dimensions because only two- and and three-dimensional" //      &
308             " variables are supported."
309          CALL inifor_abort('get_netcdf_start_and_count', message)
310
311       END IF
312
313       start = (/ 1, 1, 1 /)
314       IF ( TRIM(in_var % name) .EQ. 'T_SO' )  THEN
315!
316!--       Skip depth = 0.0 for T_SO and reduce number of depths from 9 to 8
317          in_var % dimlen(3) = in_var % dimlen(3) - 1
318
319!
320!--       Start reading from second level, e.g. depth = 0.005 instead of 0.0
321          start(3) = 2
322       END IF
323
324       IF (in_var % ndim .EQ. 2)  THEN
325          in_var % dimlen(3) = 1
326       ENDIF
327
328       ndim = MIN(in_var % ndim, 3)
329       count = (/ 1, 1, 1 /)
330       count(1:ndim) = in_var % dimlen(1:ndim)
331
332    END SUBROUTINE get_netcdf_start_and_count
333
334
335!------------------------------------------------------------------------------!
336! Description:
337! ------------
338!> Routine for defining netCDF variables in the dynamic driver, INIFOR's netCDF
339!> output.
340!------------------------------------------------------------------------------!
341    SUBROUTINE netcdf_define_variable(var, ncid)
342
343        TYPE(nc_var), INTENT(INOUT) ::  var
344        INTEGER, INTENT(IN)         ::  ncid
345
346        CALL check(nf90_def_var(ncid, var % name, NF90_FLOAT,       var % dimids(1:var % ndim), var % varid))
347        CALL check(nf90_put_att(ncid, var % varid, "long_name",     var % long_name))
348        CALL check(nf90_put_att(ncid, var % varid, "units",         var % units))
349        IF ( var % lod .GE. 0 )  THEN
350           CALL check(nf90_put_att(ncid, var % varid, "lod",           var % lod))
351        END IF
352        CALL check(nf90_put_att(ncid, var % varid, "source",        var % source))
353        CALL check(nf90_put_att(ncid, var % varid, "_FillValue",    NF90_FILL_REAL))
354
355    END SUBROUTINE netcdf_define_variable
356   
357
358!------------------------------------------------------------------------------!
359! Description:
360! ------------
361!> netcdf_get_dimensions() reads in all dimensions and their lengths and stores
362!> them in the given the 'var' data structure. This information is used later
363!> for writing output variables in update_output().
364!------------------------------------------------------------------------------!
365    SUBROUTINE netcdf_get_dimensions(var, ncid)
366
367        TYPE(nc_var), INTENT(INOUT) ::  var
368        INTEGER, INTENT(IN)         ::  ncid
369        INTEGER                     ::  i
370        CHARACTER(SNAME)            ::  null
371
372        DO i = 1, var % ndim
373           CALL check(nf90_inquire_dimension(ncid, var % dimids(i), &
374                                             name = null, &
375                                             len  = var % dimlen(i)  ) )
376        END DO
377
378    END SUBROUTINE netcdf_get_dimensions
379
380
381!------------------------------------------------------------------------------!
382! Description:
383! ------------
384!> This routine parses and interpretes the command-line options and stores the
385!> resulting settings in the 'cfg' data structure.
386!------------------------------------------------------------------------------!
387    SUBROUTINE parse_command_line_arguments( cfg )
388
389       TYPE(inifor_config), INTENT(INOUT) ::  cfg
390
391       CHARACTER(LEN=PATH)                ::  option, arg
392       INTEGER                            ::  arg_count, i
393
394       cfg % p0_is_set = .FALSE.
395       cfg % ug_is_set = .FALSE.
396       cfg % vg_is_set = .FALSE.
397       cfg % flow_prefix_is_set = .FALSE.
398       cfg % input_prefix_is_set = .FALSE.
399       cfg % radiation_prefix_is_set = .FALSE.
400       cfg % soil_prefix_is_set = .FALSE.
401       cfg % soilmoisture_prefix_is_set = .FALSE.
402
403       arg_count = COMMAND_ARGUMENT_COUNT()
404       IF (arg_count .GT. 0)  THEN
405
406          message = "The -clon and -clat command line options are depricated. " // &
407             "Please remove them form your inifor command and specify the " // &
408             "location of the PALM-4U origin either" // NEW_LINE(' ') // &
409             "   - by setting the namelist parameters 'longitude' and 'latitude', or" // NEW_LINE(' ') // &
410             "   - by providing a static driver netCDF file via the -static command-line option."
411
412          i = 1
413          DO WHILE (i .LE. arg_count)
414
415             CALL GET_COMMAND_ARGUMENT( i, option )
416
417             SELECT CASE( TRIM(option) )
418
419             CASE( '--averaging-mode' )
420                CALL get_option_argument( i, arg )
421                cfg % averaging_mode = TRIM(arg)
422
423             CASE( '-date', '-d', '--date' )
424                CALL get_option_argument( i, arg )
425                cfg % start_date = TRIM(arg)
426
427             CASE( '--debug' )
428                cfg % debug = .TRUE.
429
430             CASE( '-z0', '-z', '--elevation' )
431                CALL get_option_argument( i, arg )
432                READ(arg, *) cfg % z0
433
434             CASE( '-p0', '-r', '--surface-pressure' )
435                cfg % p0_is_set = .TRUE.
436                CALL get_option_argument( i, arg )
437                READ(arg, *) cfg % p0
438
439             CASE( '-ug', '-u', '--geostrophic-u' )
440                cfg % ug_is_set = .TRUE.
441                CALL get_option_argument( i, arg )
442                READ(arg, *) cfg % ug
443
444             CASE( '-vg', '-v', '--geostrophic-v' )
445                cfg % vg_is_set = .TRUE.
446                CALL get_option_argument( i, arg )
447                READ(arg, *) cfg % vg
448
449             CASE( '-clon', '-clat' )
450                CALL inifor_abort('parse_command_line_arguments', message)
451
452             CASE( '-path', '-p', '--path' )
453                CALL get_option_argument( i, arg )
454                 cfg % input_path = TRIM(arg)
455
456             CASE( '-hhl', '-l', '--hhl-file' )
457                CALL get_option_argument( i, arg )
458                cfg % hhl_file = TRIM(arg)
459
460             CASE( '--input-prefix')
461                CALL get_option_argument( i, arg )
462                cfg % input_prefix = TRIM(arg)
463                cfg % input_prefix_is_set = .TRUE.
464   
465             CASE( '-a', '--averaging-angle' )
466                CALL get_option_argument( i, arg )
467                READ(arg, *) cfg % averaging_angle
468
469             CASE( '-static', '-t', '--static-driver' )
470                CALL get_option_argument( i, arg )
471                cfg % static_driver_file = TRIM(arg)
472
473             CASE( '-soil', '-s', '--soil-file')
474                CALL get_option_argument( i, arg )
475                cfg % soiltyp_file = TRIM(arg)
476
477             CASE( '--flow-prefix')
478                CALL get_option_argument( i, arg )
479                cfg % flow_prefix = TRIM(arg)
480                cfg % flow_prefix_is_set = .TRUE.
481   
482             CASE( '--radiation-prefix')
483                CALL get_option_argument( i, arg )
484                cfg % radiation_prefix = TRIM(arg)
485                cfg % radiation_prefix_is_set = .TRUE.
486   
487             CASE( '--soil-prefix')
488                CALL get_option_argument( i, arg )
489                cfg % soil_prefix = TRIM(arg)
490                cfg % soil_prefix_is_set = .TRUE.
491   
492             CASE( '--soilmoisture-prefix')
493                CALL get_option_argument( i, arg )
494                cfg % soilmoisture_prefix = TRIM(arg)
495                cfg % soilmoisture_prefix_is_set = .TRUE.
496
497             CASE( '-o', '--output' )
498                CALL get_option_argument( i, arg )
499                cfg % output_file = TRIM(arg)
500
501             CASE( '-n', '--namelist' )
502                CALL get_option_argument( i, arg )
503                cfg % namelist_file = TRIM(arg)
504
505             CASE( '-mode', '-i', '--init-mode' )
506                CALL get_option_argument( i, arg )
507                cfg % ic_mode = TRIM(arg)
508
509             CASE( '-f', '--forcing-mode' )
510                CALL get_option_argument( i, arg )
511                cfg % bc_mode = TRIM(arg)
512
513             CASE( '--version' )
514                CALL print_version()
515                STOP
516
517             CASE( '--help' )
518                CALL print_version()
519                PRINT *, ""
520                PRINT *, "For a list of command-line options have a look at the README file."
521                STOP
522
523             CASE DEFAULT
524                message = "unknown option '" // TRIM(option) // "'."
525                CALL inifor_abort('parse_command_line_arguments', message)
526
527             END SELECT
528
529             i = i + 1
530
531          END DO
532
533       ELSE
534           
535          message = "No arguments present, using default input and output files"
536          CALL report('parse_command_line_arguments', message)
537
538       END IF
539
540   END SUBROUTINE parse_command_line_arguments
541
542   
543
544   SUBROUTINE get_datetime_file_list( start_date_string, start_hour, end_hour, &
545                                      step_hour, input_path, prefix, suffix,   &
546                                      file_list )
547
548      CHARACTER (LEN=DATE), INTENT(IN) ::  start_date_string
549      CHARACTER (LEN=*),    INTENT(IN) ::  prefix, suffix, input_path
550      INTEGER,              INTENT(IN) ::  start_hour, end_hour, step_hour
551      CHARACTER(LEN=*), ALLOCATABLE, INTENT(INOUT) ::  file_list(:)
552
553      INTEGER             ::  number_of_intervals, hour, i
554      CHARACTER(LEN=DATE) ::  date_string
555
556      number_of_intervals = CEILING( REAL(end_hour - start_hour) / step_hour )
557      ALLOCATE( file_list(number_of_intervals + 1) )
558
559      DO i = 0, number_of_intervals
560
561         hour = start_hour + i * step_hour
562         date_string = add_hours_to(start_date_string, hour)
563
564         file_list(i+1) = TRIM(input_path) // TRIM(prefix) //                  &
565                          TRIM(date_string) // TRIM(suffix) // '.nc'
566
567         message = "Set up input file name '" // TRIM(file_list(i+1)) // "'"
568         CALL report('input_file_list', message)
569
570      END DO
571
572   END SUBROUTINE get_datetime_file_list
573
574!------------------------------------------------------------------------------!
575! Description:
576! ------------
577!> Establish a list of files based on the given start and end times and file
578!> prefixes and suffixes.
579!------------------------------------------------------------------------------!
580   SUBROUTINE get_input_file_list( start_date_string, start_hour, end_hour,    &
581                                   step_hour, input_path, prefix, suffix,      &
582                                   file_list, nocheck )
583
584      CHARACTER (LEN=DATE), INTENT(IN) ::  start_date_string
585      CHARACTER (LEN=*),    INTENT(IN) ::  prefix, suffix, input_path
586      INTEGER,              INTENT(IN) ::  start_hour, end_hour, step_hour
587      CHARACTER(LEN=*), ALLOCATABLE, INTENT(INOUT) ::  file_list(:)
588      LOGICAL, OPTIONAL, INTENT(IN)    ::  nocheck
589
590      INTEGER             ::  number_of_intervals, hour, i
591      CHARACTER(LEN=DATE) ::  date_string
592      CHARACTER(LEN=PATH) ::  file_name
593      LOGICAL             ::  check_files
594
595      CALL get_datetime_file_list( start_date_string, start_hour, end_hour,    &
596                                   step_hour, input_path, prefix, suffix,      &
597                                   file_list )
598
599      check_files = .TRUE.
600      IF ( PRESENT ( nocheck ) )  THEN
601         IF ( nocheck )  check_files = .FALSE.
602      END IF
603
604      IF ( check_files )  THEN
605
606         tip = "Please check if you specified the correct file prefix " //     &
607               "using the options --input-prefix, --flow-prefix, etc."
608
609         DO i = 1, SIZE(file_list)
610             CALL verify_file(file_list(i), 'input', tip)
611         END DO
612
613      END IF
614
615   END SUBROUTINE get_input_file_list
616
617
618!------------------------------------------------------------------------------!
619! Description:
620! ------------
621!> Abort INIFOR if the given file is not present.
622!------------------------------------------------------------------------------!
623   SUBROUTINE verify_file(file_name, file_kind, tip)
624
625      CHARACTER(LEN=*), INTENT(IN)           ::  file_name, file_kind
626      CHARACTER(LEN=*), INTENT(IN), OPTIONAL ::  tip
627
628      IF (.NOT. file_present(file_name))  THEN
629
630         IF (LEN(TRIM(file_name)) == 0)  THEN
631
632            message = "No name was given for the " // TRIM(file_kind) // " file."
633
634         ELSE
635
636            message = "The " // TRIM(file_kind) // " file '" //                &
637                      TRIM(file_name) // "' was not found."
638
639            IF (PRESENT(tip))  THEN
640               message = TRIM(message) // " " // TRIM(tip)
641            END IF
642
643         END IF
644
645         CALL inifor_abort('verify_file', message)
646
647      END IF
648
649   END SUBROUTINE verify_file
650
651
652!------------------------------------------------------------------------------!
653! Description:
654! ------------
655!> Get the argument of the i'th command line option, which is at the location
656!> i+1 of the argument list.
657!------------------------------------------------------------------------------!
658   SUBROUTINE get_option_argument(i, arg)
659      CHARACTER(LEN=PATH), INTENT(INOUT) ::  arg
660      INTEGER, INTENT(INOUT)             ::  i
661
662      i = i + 1
663      CALL GET_COMMAND_ARGUMENT(i, arg)
664
665   END SUBROUTINE
666
667
668!------------------------------------------------------------------------------!
669! Description:
670! ------------
671!> Checks the INIFOR configuration 'cfg' for plausibility.
672!------------------------------------------------------------------------------!
673   SUBROUTINE validate_config(cfg)
674      TYPE(inifor_config), INTENT(IN) ::  cfg
675
676      CALL verify_file(cfg % hhl_file, 'HHL')
677      CALL verify_file(cfg % namelist_file, 'NAMELIST')
678      CALL verify_file(cfg % soiltyp_file, 'SOILTYP')
679
680!
681!--   Only check optional static driver file name, if it has been given.
682      IF (TRIM(cfg % static_driver_file) .NE. '')  THEN
683         CALL verify_file(cfg % static_driver_file, 'static driver')
684      END IF
685
686      SELECT CASE( TRIM(cfg % ic_mode) )
687      CASE( 'profile', 'volume')
688      CASE DEFAULT
689         message = "Initialization mode '" // TRIM(cfg % ic_mode) //&
690                   "' is not supported. " //&
691                   "Please select either 'profile' or 'volume', " //&
692                   "or omit the -i/--init-mode/-mode option entirely, which corresponds "//&
693                   "to the latter."
694         CALL inifor_abort( 'validate_config', message )
695      END SELECT
696
697
698      SELECT CASE( TRIM(cfg % bc_mode) )
699      CASE( 'real', 'ideal')
700      CASE DEFAULT
701         message = "Forcing mode '" // TRIM(cfg % bc_mode) //& 
702                   "' is not supported. " //&
703                   "Please select either 'real' or 'ideal', " //&
704                   "or omit the -f/--forcing-mode option entirely, which corresponds "//&
705                   "to the latter."
706         CALL inifor_abort( 'validate_config', message )
707      END SELECT
708
709      SELECT CASE( TRIM(cfg % averaging_mode) )
710      CASE( 'level', 'height')
711      CASE DEFAULT
712         message = "Averaging mode '" // TRIM(cfg % averaging_mode) //&
713                   "' is not supported. " //&
714                   "Please select either 'height' or 'level', " //&
715                   "or omit the --averaging-mode option entirely, which corresponds "//&
716                   "to the latter."
717         CALL inifor_abort( 'validate_config', message )
718      END SELECT
719
720      IF ( cfg % ug_is_set .NEQV. cfg % vg_is_set )  THEN
721         message = "You specified only one component of the geostrophic " // &
722                   "wind. Please specify either both or none."
723         CALL inifor_abort( 'validate_config', message )
724      END IF
725
726   END SUBROUTINE validate_config
727
728
729!------------------------------------------------------------------------------!
730! Description:
731! ------------
732!> Check whether the given file is present on the filesystem.
733!------------------------------------------------------------------------------!
734   LOGICAL FUNCTION file_present(filename)
735      CHARACTER(LEN=PATH), INTENT(IN) ::  filename
736
737      INQUIRE(FILE=filename, EXIST=file_present)
738
739   END FUNCTION file_present
740
741
742!------------------------------------------------------------------------------!
743! Description:
744! ------------
745!> This routine initializes the dynamic driver file, i.e. INIFOR's netCDF output
746!> file.
747!>
748!> Besides writing metadata, such as global attributes, coordinates, variables,
749!> in the netCDF file, various netCDF IDs are saved for later, when INIFOR
750!> writes the actual data.
751!------------------------------------------------------------------------------!
752   SUBROUTINE setup_netcdf_dimensions(output_file, palm_grid,                  &
753                                      start_date_string, origin_lon, origin_lat)
754
755       TYPE(nc_file), INTENT(INOUT)      ::  output_file
756       TYPE(grid_definition), INTENT(IN) ::  palm_grid
757       CHARACTER (LEN=DATE), INTENT(IN)  ::  start_date_string
758       REAL(dp), INTENT(IN)              ::  origin_lon, origin_lat
759
760       CHARACTER (LEN=8)     ::  date_string
761       CHARACTER (LEN=10)    ::  time_string
762       CHARACTER (LEN=5)     ::  zone_string
763       CHARACTER (LEN=SNAME) ::  history_string
764       INTEGER               ::  ncid, nx, ny, nz, nt, dimids(3), dimvarids(3)
765       REAL(dp)              ::  z0
766
767       message = "Initializing PALM-4U dynamic driver file '" //               &
768                 TRIM(output_file % name) // "' and setting up dimensions."
769       CALL report('setup_netcdf_dimensions', message)
770
771!
772!--    Create the netCDF file as in netCDF-4/HDF5 format if __netcdf4 preprocessor flag is given
773#if defined( __netcdf4 )
774       CALL check(nf90_create(TRIM(output_file % name), OR(NF90_CLOBBER, NF90_HDF5), ncid))
775#else
776       CALL check(nf90_create(TRIM(output_file % name), NF90_CLOBBER, ncid))
777#endif
778
779!------------------------------------------------------------------------------
780!- Section 1: Define NetCDF dimensions and coordinates
781!------------------------------------------------------------------------------
782       nt = SIZE(output_file % time)
783       nx = palm_grid % nx
784       ny = palm_grid % ny
785       nz = palm_grid % nz
786       z0 = palm_grid % z0
787
788
789!
790!------------------------------------------------------------------------------
791!- Section 2: Write global NetCDF attributes
792!------------------------------------------------------------------------------
793       CALL date_and_time(DATE=date_string, TIME=time_string, ZONE=zone_string)
794       history_string =                                                        &
795           'Created on '// date_string      //                                 &
796           ' at '       // time_string(1:2) // ':' // time_string(3:4) //      &
797           ' (UTC'      // zone_string // ')'
798
799       CALL check(nf90_put_att(ncid, NF90_GLOBAL, 'title',          'PALM input file for scenario ...'))
800       CALL check(nf90_put_att(ncid, NF90_GLOBAL, 'institution',    'Deutscher Wetterdienst, Offenbach'))
801       CALL check(nf90_put_att(ncid, NF90_GLOBAL, 'author',         'Eckhard Kadasch, eckhard.kadasch@dwd.de'))
802       CALL check(nf90_put_att(ncid, NF90_GLOBAL, 'history',        TRIM(history_string)))
803       CALL check(nf90_put_att(ncid, NF90_GLOBAL, 'references',     '--'))
804       CALL check(nf90_put_att(ncid, NF90_GLOBAL, 'comment',        '--'))
805       CALL check(nf90_put_att(ncid, NF90_GLOBAL, 'origin_lat',     TRIM(real_to_str(origin_lat*TO_DEGREES, '(F18.13)'))))
806       CALL check(nf90_put_att(ncid, NF90_GLOBAL, 'origin_lon',     TRIM(real_to_str(origin_lon*TO_DEGREES, '(F18.13)'))))
807!
808!--    FIXME: This is the elevation relative to COSMO-DE/D2 sea level and does
809!--    FIXME: not necessarily comply with DHHN2016 (c.f. PALM Input Data
810!--    FIXME: Standard v1.9., origin_z)
811       CALL check(nf90_put_att(ncid, NF90_GLOBAL, 'origin_z',       TRIM(real_to_str(z0, '(F18.13)'))))
812       CALL check(nf90_put_att(ncid, NF90_GLOBAL, 'inifor_version', TRIM(VERSION)))
813       CALL check(nf90_put_att(ncid, NF90_GLOBAL, 'palm_version',   '--'))
814
815!
816!
817!------------------------------------------------------------------------------
818!- Section 2a: Define dimensions for cell centers (scalars in soil and atmosph.)
819!------------------------------------------------------------------------------
820!
821!--    reset dimids first
822       dimids = (/0, 0, 0/)
823       CALL check( nf90_def_dim(ncid, "x", nx+1, dimids(1)) )
824       CALL check( nf90_def_dim(ncid, "y", ny+1, dimids(2)) )
825       CALL check( nf90_def_dim(ncid, "z", nz, dimids(3)) )
826!
827!--    save dimids for later
828       output_file % dimids_scl = dimids 
829
830!
831!--    reset dimvarids first
832       dimvarids = (/0, 0, 0/)
833       CALL check(nf90_def_var(ncid, "x", NF90_FLOAT, dimids(1), dimvarids(1)))
834       CALL check(nf90_put_att(ncid, dimvarids(1), "standard_name", "x coordinate of cell centers"))
835       CALL check(nf90_put_att(ncid, dimvarids(1), "units", "m"))
836
837       CALL check(nf90_def_var(ncid, "y", NF90_FLOAT, dimids(2), dimvarids(2)))
838       CALL check(nf90_put_att(ncid, dimvarids(2), "standard_name", "y coordinate of cell centers"))
839       CALL check(nf90_put_att(ncid, dimvarids(2), "units", "m"))
840
841       CALL check(nf90_def_var(ncid, "z", NF90_FLOAT, dimids(3), dimvarids(3)))
842       CALL check(nf90_put_att(ncid, dimvarids(3), "standard_name", "z coordinate of cell centers"))
843       CALL check(nf90_put_att(ncid, dimvarids(3), "units", "m"))
844!
845!--    save dimvarids for later
846       output_file % dimvarids_scl = dimvarids
847
848!
849!--    overwrite third dimid with the one of depth
850       CALL check(nf90_def_dim(ncid, "zsoil", SIZE(palm_grid % depths), dimids(3)) )
851!
852!--    save dimids for later
853       output_file % dimids_soil = dimids
854
855!
856!--    overwrite third dimvarid with the one of depth
857       CALL check(nf90_def_var(ncid, "zsoil", NF90_FLOAT, output_file % dimids_soil(3), dimvarids(3)))
858       CALL check(nf90_put_att(ncid, dimvarids(3), "standard_name", "depth_below_land"))
859       CALL check(nf90_put_att(ncid, dimvarids(3), "positive", "down"))
860       CALL check(nf90_put_att(ncid, dimvarids(3), "units", "m"))
861!
862!--    save dimvarids for later
863       output_file % dimvarids_soil = dimvarids
864!
865!------------------------------------------------------------------------------
866!- Section 2b: Define dimensions for cell faces/velocities
867!------------------------------------------------------------------------------
868!
869!--    reset dimids first
870       dimids = (/0, 0, 0/)
871       CALL check(nf90_def_dim(ncid, "xu", nx, dimids(1)) )
872       CALL check(nf90_def_dim(ncid, "yv", ny, dimids(2)) )
873       CALL check(nf90_def_dim(ncid, "zw", nz-1, dimids(3)) )
874!
875!--    save dimids for later
876       output_file % dimids_vel = dimids
877
878!
879!--    reset dimvarids first
880       dimvarids = (/0, 0, 0/)
881       CALL check(nf90_def_var(ncid, "xu", NF90_FLOAT, dimids(1), dimvarids(1)))
882       CALL check(nf90_put_att(ncid, dimvarids(1), "standard_name", "x coordinate of cell faces"))
883       CALL check(nf90_put_att(ncid, dimvarids(1), "units", "m"))
884
885       CALL check(nf90_def_var(ncid, "yv", NF90_FLOAT, dimids(2), dimvarids(2)))
886       CALL check(nf90_put_att(ncid, dimvarids(2), "standard_name", "y coordinate of cell faces"))
887       CALL check(nf90_put_att(ncid, dimvarids(2), "units", "m"))
888
889       CALL check(nf90_def_var(ncid, "zw", NF90_FLOAT, dimids(3), dimvarids(3)))
890       CALL check(nf90_put_att(ncid, dimvarids(3), "standard_name", "z coordinate of cell faces"))
891       CALL check(nf90_put_att(ncid, dimvarids(3), "units", "m"))
892!
893!--    save dimvarids for later
894       output_file % dimvarids_vel = dimvarids
895
896!
897!------------------------------------------------------------------------------
898!- Section 2c: Define time dimension
899!------------------------------------------------------------------------------
900       CALL check(nf90_def_dim(ncid, "time", nt, output_file % dimid_time) )
901       CALL check(nf90_def_var(ncid, "time", NF90_FLOAT, &
902                                             output_file % dimid_time, &
903                                             output_file % dimvarid_time))
904       CALL check(nf90_put_att(ncid, output_file % dimvarid_time, "standard_name", "time"))
905       CALL check(nf90_put_att(ncid, output_file % dimvarid_time, "long_name", "time"))
906       CALL check(nf90_put_att(ncid, output_file % dimvarid_time, "units",     &
907                               "seconds since " // start_date_string // " UTC"))
908
909       CALL check(nf90_enddef(ncid))
910
911!
912!------------------------------------------------------------------------------
913!- Section 3: Write grid coordinates
914!------------------------------------------------------------------------------
915       CALL check(nf90_put_var(ncid, output_file % dimvarids_scl(1), palm_grid%x))
916       CALL check(nf90_put_var(ncid, output_file % dimvarids_scl(2), palm_grid%y))
917       CALL check(nf90_put_var(ncid, output_file % dimvarids_scl(3), palm_grid%z))
918
919       CALL check(nf90_put_var(ncid, output_file % dimvarids_vel(1), palm_grid%xu))
920       CALL check(nf90_put_var(ncid, output_file % dimvarids_vel(2), palm_grid%yv))
921       CALL check(nf90_put_var(ncid, output_file % dimvarids_vel(3), palm_grid%zw))
922
923!
924!--    TODO Read in soil depths from input file before this.
925       CALL check(nf90_put_var(ncid, output_file % dimvarids_soil(3), palm_grid%depths))
926
927!
928!--    Write time vector
929       CALL check(nf90_put_var(ncid, output_file % dimvarid_time, output_file % time))
930
931!
932!--    Close the file
933       CALL check(nf90_close(ncid))
934
935    END SUBROUTINE setup_netcdf_dimensions
936
937
938!------------------------------------------------------------------------------!
939! Description:
940! ------------
941!> Defines the netCDF variables to be written to the dynamic driver file
942!------------------------------------------------------------------------------!
943    SUBROUTINE setup_netcdf_variables(filename, output_variable_table, debug)
944
945       CHARACTER (LEN=*), INTENT(IN)        ::  filename
946       TYPE(nc_var), INTENT(INOUT), TARGET  ::  output_variable_table(:)
947       LOGICAL, INTENT(IN)                  ::  debug
948
949       TYPE(nc_var), POINTER                ::  var
950       INTEGER                              ::  i, ncid
951       LOGICAL                              ::  to_be_written
952
953       message = "Defining variables in dynamic driver '" // TRIM(filename) // "'."
954       CALL report('setup_netcdf_variables', message)
955
956       CALL check(nf90_open(TRIM(filename), NF90_WRITE, ncid))
957       CALL check(nf90_redef(ncid))
958
959       DO i = 1, SIZE(output_variable_table)
960
961          var => output_variable_table(i)
962
963          to_be_written = ( var % to_be_processed  .AND. .NOT. var % is_internal) .OR.                        &
964                          ( var % is_internal  .AND.  debug )
965
966          IF ( to_be_written )  THEN
967             message = "  variable #" // TRIM(str(i)) // " '" // TRIM(var%name) // "'."
968             CALL report('setup_netcdf_variables', message)
969
970             CALL netcdf_define_variable(var, ncid)
971             CALL netcdf_get_dimensions(var, ncid)
972          END IF
973           
974       END DO
975
976       CALL check(nf90_enddef(ncid))
977       CALL check(nf90_close(ncid))
978
979       message = "Dynamic driver '" // TRIM(filename) // "' initialized successfully."
980       CALL report('setup_netcdf_variables', message)
981
982    END SUBROUTINE setup_netcdf_variables
983
984
985!------------------------------------------------------------------------------!
986! Description:
987! ------------
988!> This routine reads and returns all input variables of the given IO group
989!> It accomodates the data by allocating a container variable that represents a
990!> list of arrays of the same length as the groups variable list. (This list
991!> will typically contain one or two items.) After the container, its members
992!> are allocated one by one with the appropriate, possibly different,
993!> dimensions.
994!>
995!> The 'group' is an INTENT(INOUT) variable so that 'get_netcdf_variable()' can
996!> record netCDF IDs in the 'in_var_list()' member variable.
997!------------------------------------------------------------------------------!
998    SUBROUTINE read_input_variables(group, iter, buffer)
999       TYPE(io_group), INTENT(INOUT), TARGET       ::  group
1000       INTEGER, INTENT(IN)                         ::  iter
1001       TYPE(container), ALLOCATABLE, INTENT(INOUT) ::  buffer(:)
1002       INTEGER                                     ::  hour, buf_id
1003       TYPE(nc_var), POINTER                       ::  input_var
1004       CHARACTER(LEN=PATH), POINTER                ::  input_file
1005       INTEGER                                     ::  ivar, nbuffers
1006
1007       message = "Reading data for I/O group '" // TRIM(group % in_var_list(1) % name) // "'."
1008       CALL report('read_input_variables', message)
1009
1010       input_file => group % in_files(iter)
1011
1012!
1013!------------------------------------------------------------------------------
1014!- Section 1: Load input buffers for accumulated variables
1015!------------------------------------------------------------------------------
1016!
1017!--    radiation budgets, precipitation
1018       IF (group % kind == 'running average' .OR.                              &
1019           group % kind == 'accumulated')  THEN
1020
1021          IF (SIZE(group % in_var_list) .GT. 1 ) THEN
1022             message = "I/O groups may not contain more than one " // & 
1023                       "accumulated variable. Group '" // TRIM(group % kind) //&
1024                       "' contains " //                                        &
1025                       TRIM( str(SIZE(group % in_var_list)) ) // "."
1026             CALL inifor_abort('read_input_variables | accumulation', message)
1027          END IF
1028
1029!
1030!--       use two buffer arrays
1031          nbuffers = 2
1032          IF ( .NOT. ALLOCATED( buffer ) )  ALLOCATE( buffer(nbuffers) )
1033
1034!
1035!--       hour of the day
1036          hour = iter - 1
1037!
1038!--       chose correct buffer array
1039          buf_id = select_buffer(hour)
1040
1041 CALL run_control('time', 'read')
1042          IF ( ALLOCATED(buffer(buf_id) % array) )  DEALLOCATE(buffer(buf_id) % array)
1043 CALL run_control('time', 'alloc')
1044
1045          input_var => group % in_var_list(1)
1046          CALL get_netcdf_variable(input_file, input_var, buffer(buf_id) % array)
1047          CALL report('read_input_variables', "Read accumulated " // TRIM(group % in_var_list(1) % name)) 
1048
1049          IF ( input_var % is_upside_down )  CALL reverse(buffer(buf_id) % array)
1050 CALL run_control('time', 'comp')
1051         
1052!------------------------------------------------------------------------------
1053!- Section 2: Load input buffers for normal I/O groups
1054!------------------------------------------------------------------------------
1055       ELSE
1056
1057!
1058!--       Allocate one input buffer per input_variable. If more quantities
1059!--       have to be computed than input variables exist in this group,
1060!--       allocate more buffers. For instance, in the thermodynamics group,
1061!--       there are three input variabels (PP, T, Qv) and four quantities
1062!--       necessart (P, Theta, Rho, qv) for the corresponding output fields
1063!--       (p0, Theta, qv, ug, and vg)
1064          nbuffers = MAX( group % n_inputs, group % n_output_quantities )
1065          ALLOCATE( buffer(nbuffers) )
1066 CALL run_control('time', 'alloc')
1067         
1068!
1069!--       Read in all input variables, leave extra buffers-if any-untouched.
1070          DO ivar = 1, group % n_inputs
1071
1072             input_var => group % in_var_list(ivar)
1073
1074!
1075!            Check wheather P or PP is present in input file
1076             IF (input_var % name == 'P')  THEN
1077                input_var % name = TRIM( get_pressure_varname(input_file) )
1078 CALL run_control('time', 'read')
1079             END IF
1080
1081             CALL get_netcdf_variable(input_file, input_var, buffer(ivar) % array)
1082
1083             IF ( input_var % is_upside_down )  CALL reverse(buffer(ivar) % array)
1084 CALL run_control('time', 'comp')
1085
1086          END DO
1087       END IF
1088
1089    END SUBROUTINE read_input_variables
1090
1091
1092!------------------------------------------------------------------------------!
1093! Description:
1094! ------------
1095!> Select the appropriate buffer ID for accumulated COSMO input variables
1096!> depending on the current hour.
1097!------------------------------------------------------------------------------!
1098    INTEGER FUNCTION select_buffer(hour)
1099       INTEGER, INTENT(IN) ::  hour
1100       INTEGER             ::  step
1101
1102       select_buffer = 0
1103       step = MODULO(hour, 3) + 1
1104
1105       SELECT CASE(step)
1106       CASE(1, 3)
1107           select_buffer = 1
1108       CASE(2)
1109           select_buffer = 2
1110       CASE DEFAULT
1111           message = "Invalid step '" // TRIM(str(step))
1112           CALL inifor_abort('select_buffer', message)
1113       END SELECT
1114    END FUNCTION select_buffer
1115
1116
1117!------------------------------------------------------------------------------!
1118! Description:
1119! ------------
1120!> Checks if the input_file contains the absolute pressure, 'P', or the pressure
1121!> perturbation, 'PP', and returns the appropriate string.
1122!------------------------------------------------------------------------------!
1123    CHARACTER(LEN=2) FUNCTION get_pressure_varname(input_file) RESULT(var)
1124       CHARACTER(LEN=*) ::  input_file
1125       INTEGER          ::  ncid, varid
1126
1127       CALL check(nf90_open( TRIM(input_file), NF90_NOWRITE, ncid ))
1128       IF ( nf90_inq_varid( ncid, 'P', varid ) .EQ. NF90_NOERR )  THEN
1129
1130          var = 'P'
1131
1132       ELSE IF ( nf90_inq_varid( ncid, 'PP', varid ) .EQ. NF90_NOERR )  THEN
1133
1134          var = 'PP'
1135          CALL report('get_pressure_var', 'Using PP instead of P')
1136
1137       ELSE
1138
1139          message = "Failed to read '" // TRIM(var) // &
1140                    "' from file '" // TRIM(input_file) // "'."
1141          CALL inifor_abort('get_pressure_var', message)
1142
1143       END IF
1144
1145       CALL check(nf90_close(ncid))
1146
1147    END FUNCTION get_pressure_varname
1148
1149
1150!------------------------------------------------------------------------------!
1151! Description:
1152! ------------
1153!> Read the given global attribute form the given netCDF file.
1154!------------------------------------------------------------------------------!
1155    FUNCTION get_netcdf_attribute(filename, attribute) RESULT(attribute_value)
1156
1157       CHARACTER(LEN=*), INTENT(IN) ::  filename, attribute
1158       REAL(dp)                     ::  attribute_value
1159
1160       INTEGER                      ::  ncid
1161
1162       IF ( nf90_open( TRIM(filename), NF90_NOWRITE, ncid ) == NF90_NOERR )  THEN
1163
1164          CALL check(nf90_get_att(ncid, NF90_GLOBAL, TRIM(attribute), attribute_value))
1165          CALL check(nf90_close(ncid))
1166
1167       ELSE
1168
1169          message = "Failed to read '" // TRIM(attribute) // &
1170                    "' from file '" // TRIM(filename) // "'."
1171          CALL inifor_abort('get_netcdf_attribute', message)
1172
1173       END IF
1174
1175    END FUNCTION get_netcdf_attribute
1176
1177
1178
1179!------------------------------------------------------------------------------!
1180! Description:
1181! ------------
1182!> Updates the dynamic driver with the interpolated field of the current
1183!> variable at the current time step.
1184!------------------------------------------------------------------------------!
1185    SUBROUTINE update_output(var, array, iter, output_file, cfg)
1186       TYPE(nc_var), INTENT(IN)  ::  var
1187       REAL(dp), INTENT(IN)      ::  array(:,:,:)
1188       INTEGER, INTENT(IN)       ::  iter
1189       TYPE(nc_file), INTENT(IN) ::  output_file
1190       TYPE(inifor_config)       ::  cfg
1191
1192       INTEGER ::  ncid, ndim, start(4), count(4)
1193       LOGICAL ::  var_is_time_dependent
1194
1195       var_is_time_dependent = (                                               &
1196          var % dimids( var % ndim ) == output_file % dimid_time               &
1197       )
1198
1199!
1200!--    Skip time dimension for output
1201       ndim = var % ndim
1202       IF ( var_is_time_dependent )  ndim = var % ndim - 1
1203
1204       start(:)      = (/1,1,1,1/)
1205       start(ndim+1) = iter
1206       count(1:ndim) = var%dimlen(1:ndim)
1207
1208       CALL check(nf90_open(output_file % name, NF90_WRITE, ncid))
1209
1210!
1211!--    Reduce dimension of output array according to variable kind
1212       SELECT CASE (TRIM(var % kind))
1213       
1214       CASE ( 'init scalar profile', 'init u profile', 'init v profile',       &
1215              'init w profile' )
1216
1217          CALL check(nf90_put_var( ncid, var%varid, array(1,1,:) ) )
1218
1219       CASE ( 'init soil', 'init scalar', 'init u', 'init v', 'init w' )
1220
1221          CALL check(nf90_put_var( ncid, var%varid, array(:,:,:) ) )
1222
1223       CASE( 'left scalar', 'right scalar', 'left w', 'right w' ) 
1224
1225          CALL check(nf90_put_var( ncid, var%varid, array(1,:,:),              &
1226                                   start=start(1:ndim+1),                      &
1227                                   count=count(1:ndim) ) )
1228         
1229
1230          IF (.NOT. SIZE(array, 2) .EQ. var % dimlen(1))  THEN
1231             PRINT *, "inifor: update_output: Dimension ", 1, " of variable ", &
1232                 TRIM(var % name), " (", var % dimlen(1),                      &
1233                 ") does not match the dimension of the output array (",       &
1234                 SIZE(array, 2), ")."
1235             STOP
1236          END IF
1237         
1238
1239       CASE( 'north scalar', 'south scalar', 'north w', 'south w' )
1240
1241          CALL check(nf90_put_var( ncid, var%varid, array(:,1,:),              &
1242                                   start=start(1:ndim+1),                      &
1243                                   count=count(1:ndim) ) )
1244         
1245
1246       CASE( 'surface forcing', 'top scalar', 'top w' )
1247
1248          CALL check(nf90_put_var( ncid, var%varid, array(:,:,1),              &
1249                                   start=start(1:ndim+1),                      &
1250                                   count=count(1:ndim) ) )
1251         
1252       CASE ( 'left u', 'right u', 'left v', 'right v' )
1253
1254          CALL check(nf90_put_var( ncid, var%varid, array(1,:,:),              &
1255                                   start=start(1:ndim+1),                      &
1256                                   count=count(1:ndim) ) )
1257
1258       CASE ( 'north u', 'south u', 'north v', 'south v' )
1259
1260          CALL check(nf90_put_var( ncid, var%varid, array(:,1,:),              &
1261                                   start=start(1:ndim+1),                      &
1262                                   count=count(1:ndim) ) )
1263
1264       CASE ( 'top u', 'top v' )
1265
1266          CALL check(nf90_put_var( ncid, var%varid, array(:,:,1),              &
1267                                   start=start(1:ndim+1),                      &
1268                                   count=count(1:ndim) ) )
1269
1270       CASE ( 'time series' )
1271
1272          CALL check(nf90_put_var( ncid, var%varid, array(1,1,1),              &
1273                                   start=start(1:ndim+1) ) )
1274
1275       CASE ( 'constant scalar profile', 'geostrophic' )
1276
1277          CALL check(nf90_put_var( ncid, var%varid, array(1,1,:),              &
1278                                   start=start(1:ndim+1),                      &
1279                                   count=count(1:ndim) ) )
1280
1281       CASE ( 'internal profile' )
1282
1283          IF ( cfg % debug )  THEN
1284             CALL check(nf90_put_var( ncid, var%varid, array(1,1,:),           &
1285                                      start=start(1:ndim+1),                   &
1286                                      count=count(1:ndim) ) )
1287          END IF
1288
1289       CASE ( 'large-scale scalar forcing', 'large-scale w forcing' )
1290
1291           message = "Doing nothing in terms of writing large-scale forings."
1292           CALL report('update_output', message)
1293
1294       CASE DEFAULT
1295
1296           message = "Variable kind '" // TRIM(var % kind) //                  &
1297                    "' not recognized."
1298           CALL inifor_abort('update_output', message)
1299
1300       END SELECT
1301
1302       CALL check(nf90_close(ncid))
1303
1304    END SUBROUTINE update_output
1305
1306
1307!------------------------------------------------------------------------------!
1308! Description:
1309! ------------
1310!> Checks the status of a netCDF API call and aborts if an error occured
1311!------------------------------------------------------------------------------!
1312    SUBROUTINE check(status)
1313
1314       INTEGER, INTENT(IN) ::  status
1315
1316       IF (status /= nf90_noerr)  THEN
1317          message = "NetCDF API call failed with error: " //                     &
1318                    TRIM( nf90_strerror(status) )
1319          CALL inifor_abort('io.check', message)
1320       END IF
1321
1322    END SUBROUTINE check
1323
1324 END MODULE inifor_io
1325#endif
Note: See TracBrowser for help on using the repository browser.