source: palm/tags/release-6.0/UTIL/inifor/src/inifor_io.f90 @ 4418

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

inifor: Removed surface forcing and internal arrays from netCDF output

  • Property svn:keywords set to Id
File size: 42.6 KB
Line 
1!> @file src/inifor_io.f90
2!------------------------------------------------------------------------------!
3! This file is part of the PALM model system.
4!
5! PALM is free software: you can redistribute it and/or modify it under the
6! terms of the GNU General Public License as published by the Free Software
7! Foundation, either version 3 of the License, or (at your option) any later
8! version.
9!
10! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
11! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
12! A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
13!
14! You should have received a copy of the GNU General Public License along with
15! PALM. If not, see <http://www.gnu.org/licenses/>.
16!
17! Copyright 2017-2018 Leibniz Universitaet Hannover
18! Copyright 2017-2018 Deutscher Wetterdienst Offenbach
19!------------------------------------------------------------------------------!
20!
21! Current revisions:
22! -----------------
23!
24!
25! Former revisions:
26! -----------------
27! $Id: inifor_io.f90 3456 2018-10-30 14:29:54Z raasch $
28! NetCDf output of internal arrays only with --debug option
29!
30!
31! 3447 2018-10-29 15:52:54Z eckhard
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
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
46!
47! 3183 2018-07-27 14:25:55Z suehring
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
59!
60!
61! 3182 2018-07-27 13:36:03Z suehring
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,                                                                  &
79        ONLY:  DATE, SNAME, PATH, PI, dp, hp, TO_RADIANS, TO_DEGREES, VERSION
80    USE netcdf
81    USE types
82    USE util,                                                                  &
83        ONLY:  reverse, str, real_to_str
84
85    IMPLICIT NONE
86
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
94 CONTAINS
95
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
102       INTEGER               ::  ncid
103       INTEGER, DIMENSION(3) ::  start, count
104
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
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
140       INTEGER               ::  ncid
141       INTEGER, DIMENSION(3) ::  start, count
142
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
169    END SUBROUTINE get_netcdf_variable_real
170
171
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
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))
243        IF ( var % lod .GE. 0 )  THEN
244           CALL check(nf90_put_att(ncid, var % varid, "lod",           var % lod))
245        END IF
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!------------------------------------------------------------------------------!
277    SUBROUTINE parse_command_line_arguments( cfg )
278
279       TYPE(inifor_config), INTENT(INOUT) ::  cfg
280
281       CHARACTER(LEN=PATH)                ::  option, arg
282       INTEGER                            ::  arg_count, i
283
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
293       arg_count = COMMAND_ARGUMENT_COUNT()
294       IF (arg_count .GT. 0)  THEN
295
296          ! Every option should have an argument.
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
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(' ') // &
305             "   - by setting the namelist parameters 'longitude' and 'latitude', or" // NEW_LINE(' ') // &
306             "   - by providing a static driver netCDF file via the -static command-line option."
307
308          i = 1
309          DO WHILE (i .LE. arg_count)
310
311             CALL GET_COMMAND_ARGUMENT( i, option )
312
313             SELECT CASE( TRIM(option) )
314
315             CASE( '--averaging-mode' )
316                CALL get_option_argument( i, arg )
317                cfg % averaging_mode = TRIM(arg)
318
319             CASE( '-date', '-d', '--date' )
320                CALL get_option_argument( i, arg )
321                cfg % start_date = TRIM(arg)
322
323             CASE( '--debug' )
324                cfg % debug = .TRUE.
325
326             ! Elevation of the PALM-4U domain above sea level
327             CASE( '-z0', '-z', '--elevation' )
328                CALL get_option_argument( i, arg )
329                READ(arg, *) cfg % z0
330
331             ! surface pressure, at z0
332             CASE( '-p0', '-r', '--surface-pressure' )
333                cfg % p0_is_set = .TRUE.
334                CALL get_option_argument( i, arg )
335                READ(arg, *) cfg % p0
336
337             ! geostrophic wind in x direction
338             CASE( '-ug', '-u', '--geostrophic-u' )
339                cfg % ug_is_set = .TRUE.
340                CALL get_option_argument( i, arg )
341                READ(arg, *) cfg % ug
342
343             ! geostrophic wind in y direction
344             CASE( '-vg', '-v', '--geostrophic-v' )
345                cfg % vg_is_set = .TRUE.
346                CALL get_option_argument( i, arg )
347                READ(arg, *) cfg % vg
348
349             ! domain centre geographical longitude and latitude
350             CASE( '-clon', '-clat' )
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
357             CASE( '-path', '-p', '--path' )
358                CALL get_option_argument( i, arg )
359                 cfg % input_path = TRIM(arg)
360
361             CASE( '-hhl', '-l', '--hhl-file' )
362                CALL get_option_argument( i, arg )
363                cfg % hhl_file = TRIM(arg)
364
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
374             CASE( '-static', '-t', '--static-driver' )
375                CALL get_option_argument( i, arg )
376                cfg % static_driver_file = TRIM(arg)
377
378             CASE( '-soil', '-s', '--soil-file')
379                CALL get_option_argument( i, arg )
380                cfg % soiltyp_file = TRIM(arg)
381
382             CASE( '--flow-prefix')
383                CALL get_option_argument( i, arg )
384                cfg % flow_prefix = TRIM(arg)
385                cfg % flow_prefix_is_set = .TRUE.
386   
387             CASE( '--radiation-prefix')
388                CALL get_option_argument( i, arg )
389                cfg % radiation_prefix = TRIM(arg)
390                cfg % radiation_prefix_is_set = .TRUE.
391   
392             CASE( '--soil-prefix')
393                CALL get_option_argument( i, arg )
394                cfg % soil_prefix = TRIM(arg)
395                cfg % soil_prefix_is_set = .TRUE.
396   
397             CASE( '--soilmoisture-prefix')
398                CALL get_option_argument( i, arg )
399                cfg % soilmoisture_prefix = TRIM(arg)
400                cfg % soilmoisture_prefix_is_set = .TRUE.
401
402             CASE( '-o', '--output' )
403                CALL get_option_argument( i, arg )
404                cfg % output_file = TRIM(arg)
405
406             CASE( '-n', '--namelist' )
407                CALL get_option_argument( i, arg )
408                cfg % namelist_file = TRIM(arg)
409
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
430             CASE DEFAULT
431                message = "unknown option '" // TRIM(option) // "'."
432                CALL abort('parse_command_line_arguments', message)
433
434             END SELECT
435
436             i = i + 1
437
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
449   
450   SUBROUTINE get_option_argument(i, arg)
451      CHARACTER(LEN=PATH), INTENT(INOUT) ::  arg
452      INTEGER, INTENT(INOUT)             ::  i
453
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.
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')
468
469      ! Only check optional static driver file name, if it has been given.
470      IF (TRIM(cfg % static_driver_file) .NE. '')  THEN
471         all_files_present = all_files_present .AND. file_present(cfg % static_driver_file, 'static driver')
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
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
513
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
520   END SUBROUTINE validate_config
521
522
523   LOGICAL FUNCTION file_present(filename, kind)
524      CHARACTER(LEN=PATH), INTENT(IN) ::  filename
525      CHARACTER(LEN=*), INTENT(IN)    ::  kind
526
527      IF (LEN(TRIM(filename))==0)  THEN
528
529         file_present = .FALSE.
530         message = "No name was given for the " // TRIM(kind) // " file."
531         CALL report('file_present', message)
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
543      END IF
544
545   END FUNCTION file_present
546
547
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!------------------------------------------------------------------------------!
558   SUBROUTINE setup_netcdf_dimensions(output_file, palm_grid,                  &
559                                      start_date_string, origin_lon, origin_lat)
560
561       TYPE(nc_file), INTENT(INOUT)      ::  output_file
562       TYPE(grid_definition), INTENT(IN) ::  palm_grid
563       CHARACTER (LEN=DATE), INTENT(IN)  ::  start_date_string
564       REAL(dp), INTENT(IN)              ::  origin_lon, origin_lat
565
566       CHARACTER (LEN=8)     ::  date_string
567       CHARACTER (LEN=10)    ::  time_string
568       CHARACTER (LEN=5)     ::  zone_string
569       CHARACTER (LEN=SNAME) ::  history_string
570       INTEGER               ::  ncid, nx, ny, nz, nt, dimids(3), dimvarids(3)
571       REAL(dp)              ::  z0
572
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
577       ! Create the NetCDF file. NF90_CLOBBER selects overwrite mode.
578#if defined( __netcdf4 )
579       CALL check(nf90_create(TRIM(output_file % name), OR(NF90_CLOBBER, NF90_HDF5), ncid))
580#else
581       CALL check(nf90_create(TRIM(output_file % name), NF90_CLOBBER, ncid))
582#endif
583
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
594!
595!------------------------------------------------------------------------------
596!- Section 2: Write global NetCDF attributes
597!------------------------------------------------------------------------------
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
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'))
607       CALL check(nf90_put_att(ncid, NF90_GLOBAL, 'history',        TRIM(history_string)))
608       CALL check(nf90_put_att(ncid, NF90_GLOBAL, 'references',     '--'))
609       CALL check(nf90_put_att(ncid, NF90_GLOBAL, 'comment',        '--'))
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)'))))
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)'))))
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)) )
627          CALL check( nf90_def_dim(ncid, "z", nz, dimids(3)) )
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
645       CALL check(nf90_def_dim(ncid, "zsoil", SIZE(palm_grid % depths), dimids(3)) )
646       output_file % dimids_soil = dimids ! save dimids for later
647
648       ! overwrite third dimvarid with the one of depth
649       CALL check(nf90_def_var(ncid, "zsoil", NF90_FLOAT, output_file % dimids_soil(3), dimvarids(3)))
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)) )
661          CALL check(nf90_def_dim(ncid, "zw", nz-1, dimids(3)) )
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"))
688       CALL check(nf90_put_att(ncid, output_file % dimvarid_time, "units",     &
689                               "seconds since " // start_date_string // " UTC"))
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
717    SUBROUTINE setup_netcdf_variables(filename, output_variable_table, debug)
718
719       CHARACTER (LEN=*), INTENT(IN)        ::  filename
720       TYPE(nc_var), INTENT(INOUT), TARGET  ::  output_variable_table(:)
721       LOGICAL, INTENT(IN)                  ::  debug
722
723       TYPE(nc_var), POINTER                ::  var
724       INTEGER                              ::  i, ncid
725       LOGICAL                              ::  to_be_written
726
727       message = "Defining variables in dynamic driver '" // TRIM(filename) // "'."
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
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
741             message = "  variable #" // TRIM(str(i)) // " '" // TRIM(var%name) // "'."
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
753       message = "Dynamic driver '" // TRIM(filename) // "' initialized successfully."
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)
814          CALL get_netcdf_variable(input_file, input_var, buffer(buf_id) % array)
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
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 )
832          ALLOCATE( buffer(nbuffers) )
833 CALL run_control('time', 'alloc')
834         
835          ! Read in all input variables, leave extra buffers-if any-untouched.
836          DO ivar = 1, group % n_inputs
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
842                input_var % name = TRIM( get_pressure_varname(input_file) )
843 CALL run_control('time', 'read')
844             END IF
845
846             CALL get_netcdf_variable(input_file, input_var, buffer(ivar) % array)
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!------------------------------------------------------------------------------!
882    CHARACTER(LEN=2) FUNCTION get_pressure_varname(input_file) RESULT(var)
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
906    END FUNCTION get_pressure_varname
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))
919          CALL check(nf90_close(ncid))
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
932    SUBROUTINE update_output(var, array, iter, output_file, cfg)
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
937       TYPE(inifor_config)       ::  cfg
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
947       ndim = var % ndim
948       IF ( var_is_time_dependent )  ndim = var % ndim - 1
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
1020       CASE ( 'constant scalar profile', 'geostrophic' )
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
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
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
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.