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

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

inifor: bugfix: Fixes issue #815 with geostrophic wind profiles

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