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

Last change on this file since 3534 was 3534, checked in by raasch, 3 years ago

inifor integrated in build mechanism, some bugfixes in inifor to avoid compile time errors, batch_scp for sending back the job protocol file is called via login-node if a login-node has been set in the config-file, ssh-calls rearranged to avoid output of system/user-profile messages

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