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

Last change on this file since 3927 was 3866, checked in by eckhard, 6 years ago

inifor: Use PALM's working precision; improved error handling, coding style, and comments

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