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

Last change on this file since 4313 was 4074, checked in by eckhard, 2 years ago

inifor: Changed initialization mode from 'volume' to 'profile'

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