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

Last change on this file since 3785 was 3785, checked in by eckhard, 3 years ago

inifor: Removed unused variables, improved coding style

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