source: palm/trunk/UTIL/inifor/src/io.f90 @ 3182

Last change on this file since 3182 was 3182, checked in by suehring, 3 years ago

New Inifor features: grid stretching, improved command-interface, support start dates in different formats in both YYYYMMDD and YYYYMMDDHH, Ability to manually control input file prefixes (--radiation-prefix, --soil-preifx, --flow-prefix, --soilmoisture-prefix) for compatiblity with DWD forcast naming scheme; GNU-style short and long option; Prepared output of large-scale forcing profiles (no computation yet); Added preprocessor flag netcdf4 to switch output format between netCDF 3 and 4; Updated netCDF variable names and attributes to comply with PIDS v1.9; Inifor bugfixes: Improved compatibility with older Intel Intel compilers by avoiding implicit array allocation; Added origin_lon/_lat values and correct reference time in dynamic driver global attributes; corresponding PALM changes: adjustments to revised Inifor; variables names in dynamic driver adjusted; enable geostrophic forcing also in offline nested mode; variable names in LES-LES and COSMO offline nesting changed; lateral boundary flags for nesting, in- and outflow conditions renamed

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