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

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

inifor: Added computation of geostrophic winds from COSMO input

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