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

Last change on this file since 4004 was 3997, checked in by eckhard, 5 years ago

inifor: Read origin_z from static driver if given; alert user to warnings

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