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

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

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

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