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

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

Removed unnecessary file check

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