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

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

last commit documented

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