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

Last change on this file since 3447 was 3447, checked in by eckhard, 3 years ago

inifor: Renamed source files for compatibilty with PALM build system

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