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

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

inifor: Removed surface forcing and internal arrays from netCDF output

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