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

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

inifor: Read COSMO rotated pole from dataset, check if PALM domain exceeds COSMO domain

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