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

Last change on this file since 4659 was 4659, checked in by eckhard, 2 years ago

inifor: Support for COSMO cloud water and precipitation

  • Property svn:keywords set to Id
File size: 63.7 KB
Line 
1!> @file src/inifor_io.f90
2!------------------------------------------------------------------------------!
3! This file is part of the PALM model system.
4!
5! PALM is free software: you can redistribute it and/or modify it under the
6! terms of the GNU General Public License as published by the Free Software
7! Foundation, either version 3 of the License, or (at your option) any later
8! version.
9!
10! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
11! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
12! A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
13!
14! You should have received a copy of the GNU General Public License along with
15! PALM. If not, see <http://www.gnu.org/licenses/>.
16!
17! Copyright 2017-2020 Leibniz Universitaet Hannover
18! Copyright 2017-2020 Deutscher Wetterdienst Offenbach
19!------------------------------------------------------------------------------!
20!
21! Current revisions:
22! -----------------
23!
24!
25! Former revisions:
26! -----------------
27! $Id: inifor_io.f90 4659 2020-08-31 11:21:17Z eckhard $
28! New command-line option '--precipitation' enables output of precipitation
29!    surface forcing
30! Only define netCDF variables in enabled IO groups
31! Improved code formatting
32!
33!
34! 4569 2020-06-19 12:46:17Z eckhard
35! Removed unused variable
36!
37!
38! 4568 2020-06-19 11:56:30Z eckhard
39! Handle COSMO soil data with and without additional surface temperature
40!
41!
42! 4553 2020-06-03 16:34:15Z eckhard
43! Option --help now points the user to INIFOR's wiki page
44! Remove deprecated command-line options -clon and -clat
45!
46!
47! 4538 2020-05-18 13:45:35Z eckhard
48! Make setting the vertical PALM origin mandatory
49!
50!
51! 4523 2020-05-07 15:58:16Z eckhard
52! respect integer working precision (iwp) specified in inifor_defs.f90
53!
54!
55! 4481 2020-03-31 18:55:54Z maronga
56! Pass hhl_file directly instead of entire INIFOR configuration
57!
58!
59! 3997 2019-05-23 12:35:57Z eckhard
60! Added boolean indicator for --elevation option invocation
61! Stop INIFOR if no command-line options given
62!
63!
64! 3866 2019-04-05 14:25:01Z eckhard
65! Use PALM's working precision
66! Improved coding style
67!
68!
69! 3801 2019-03-15 17:14:25Z eckhard
70! Added routine get_cosmo_grid() to read in COSMO rotated pole from COSMO domain
71! Moved get_soil_layer_thickness() here from inifor_grid
72!
73! 3785 2019-03-06 10:41:14Z eckhard
74! Temporariliy disabled height-based geostrophic wind averaging
75! Improved variable naming
76!
77!
78! 3764 2019-02-26 13:42:09Z eckhard
79! Removed dependency on radiation input files
80!
81!
82! 3716 2019-02-05 17:02:38Z eckhard
83! Removed dependency on soilmoisture input files
84!
85!
86! 3680 2019-01-18 14:54:12Z knoop
87! Moved get_input_file_list() here from grid module, added check for presence of
88!    input files
89!
90!
91! 3618 2018-12-10 13:25:22Z eckhard
92! Prefixed all INIFOR modules with inifor_
93!
94!
95! 3615 2018-12-10 07:21:03Z raasch
96! bugfix: abort replaced by inifor_abort
97!
98! 3557 2018-11-22 16:01:22Z eckhard
99! Updated documentation, removed unused subroutine write_netcdf_variable_2d()
100!
101!
102! 3537 2018-11-20 10:53:14Z eckhard
103! New routine get_netcdf_dim_vector()
104!
105!
106! 3534 2018-11-19 15:35:16Z raasch
107! bugfix: INTENT attribute changed
108!
109! 3456 2018-10-30 14:29:54Z eckhard
110! NetCDf output of internal arrays only with --debug option
111!
112!
113! 3447 2018-10-29 15:52:54Z eckhard
114! Removed INCLUDE statement for get_netcdf_variable()
115! Renamed source files for compatibilty with PALM build system
116!
117!
118! 3395 2018-10-22 17:32:49Z eckhard
119! Added command-line options for configuring the computation of geostrophic
120!     winds (--averagin-mode, --averaging-angle)
121! Added command-line option --input-prefix for setting input file prefixes all
122!     at once
123! Added --debug option for more verbose terminal output
124! Added option-specific *_is_set LOGICALs to indicate invocation from the
125!     command-line
126! Improved error messages in case of empty file-name strings
127! Improved routine naming
128!
129! 3183 2018-07-27 14:25:55Z suehring
130! Introduced new PALM grid stretching
131! Updated variable names and metadata for PIDS v1.9 compatibility
132! Improved handling of the start date string
133! Better compatibility with older Intel compilers:
134! - avoiding implicit array allocation with new get_netcdf_variable()
135!   subroutine instead of function
136! Improved command line interface:
137! - Added configuration validation
138! - New options to configure input file prefixes
139! - GNU-style short and long option names
140! - Added version and copyright output
141!
142!
143! 3182 2018-07-27 13:36:03Z suehring
144! Initial revision
145!
146!
147!
148! Authors:
149! --------
150!> @author Eckhard Kadasch (Deutscher Wetterdienst, Offenbach)
151!
152! Description:
153! ------------
154!> The io module contains the functions needed to read and write netCDF data in
155!> INIFOR.
156!------------------------------------------------------------------------------!
157#if defined ( __netcdf )
158 MODULE inifor_io
159
160    USE inifor_control
161    USE inifor_defs,                                                           &
162        ONLY:  DATE, SNAME, PATH, PI, TO_RADIANS, TO_DEGREES, VERSION,         &
163               NC_DEPTH_DIM_IDX, NC_DEPTH_NAME, NC_HHL_NAME, NC_RLAT_NAME,     &
164               NC_RLON_NAME, NC_ROTATED_POLE_NAME, NC_POLE_LATITUDE_NAME,      &
165               NC_POLE_LONGITUDE_NAME, RHO_L, iwp, wp,                         &
166               PIDS_ORIGIN_LON, PIDS_ORIGIN_LAT, PIDS_ORIGIN_Z
167    USE inifor_types
168    USE inifor_util,                                                           &
169        ONLY:  add_hours_to, nearly_equal, reverse, str, real_to_str
170    USE netcdf
171
172    IMPLICIT NONE
173
174!------------------------------------------------------------------------------!
175! Description:
176! ------------
177!> get_netcdf_variable() reads the netCDF data and metadate for the netCDF
178!> variable 'in_var%name' from the file 'in_file'. The netCDF data array is
179!> stored in the 'buffer' array and metadata added to the respective members of
180!> the given 'in_var'.
181!------------------------------------------------------------------------------!
182    INTERFACE get_netcdf_variable
183       MODULE PROCEDURE get_netcdf_variable_int
184       MODULE PROCEDURE get_netcdf_variable_real
185    END INTERFACE get_netcdf_variable
186
187    PRIVATE ::  get_netcdf_variable_int, get_netcdf_variable_real
188
189 CONTAINS
190
191!------------------------------------------------------------------------------!
192! Description:
193! ------------
194!> get_netcdf_variable_int() implements the integer variant for the
195!> get_netcdf_variable interface.
196!------------------------------------------------------------------------------!
197 SUBROUTINE get_netcdf_variable_int(in_file, in_var, buffer)
198
199    CHARACTER(LEN=PATH), INTENT(IN)          ::  in_file
200    TYPE(nc_var), INTENT(INOUT)              ::  in_var
201    INTEGER(iwp), ALLOCATABLE, INTENT(INOUT) ::  buffer(:,:,:)
202
203    INTEGER               ::  ncid
204    INTEGER, DIMENSION(3) ::  start, count
205
206    IF ( nf90_open( TRIM(in_file), NF90_NOWRITE, ncid ) .EQ. NF90_NOERR .AND. &
207         nf90_inq_varid( ncid, in_var%name, in_var%varid ) .EQ. NF90_NOERR )  THEN
208
209       CALL get_input_dimensions(in_var, ncid)
210
211       CALL get_netcdf_start_and_count(in_var, start, count)
212       CALL log_runtime('time', 'read')
213
214       ALLOCATE( buffer( count(1), count(2), count(3) ) )
215       CALL log_runtime('time', 'alloc')
216
217       CALL check(nf90_get_var( ncid, in_var%varid, buffer,                  &
218                                start = start,                                 &
219                                count = count ))
220
221    ELSE
222
223       message = "Failed to read '" // TRIM(in_var%name) // &
224          "' from file '" // TRIM(in_file) // "'."
225       CALL inifor_abort('get_netcdf_variable', message)
226
227    ENDIF
228
229    CALL check(nf90_close(ncid))
230    CALL log_runtime('time', 'read')
231
232 END SUBROUTINE get_netcdf_variable_int
233
234
235!------------------------------------------------------------------------------!
236! Description:
237! ------------
238!> get_netcdf_variable_real() implements the real variant for the
239!> get_netcdf_variable interface.
240!------------------------------------------------------------------------------!
241 SUBROUTINE get_netcdf_variable_real(in_file, in_var, buffer)
242
243    CHARACTER(LEN=PATH), INTENT(IN)      ::  in_file
244    TYPE(nc_var), INTENT(INOUT)          ::  in_var
245    REAL(wp), ALLOCATABLE, INTENT(INOUT) ::  buffer(:,:,:)
246
247    INTEGER               ::  ncid
248    INTEGER, DIMENSION(3) ::  start, count
249
250    IF ( nf90_open( TRIM(in_file), NF90_NOWRITE, ncid ) .EQ. NF90_NOERR .AND. &
251         nf90_inq_varid( ncid, in_var%name, in_var%varid ) .EQ. NF90_NOERR )  THEN
252
253       CALL get_input_dimensions(in_var, ncid)
254
255       CALL get_netcdf_start_and_count(in_var, start, count)
256       CALL log_runtime('time', 'read')
257
258       ALLOCATE( buffer( count(1), count(2), count(3) ) )
259       CALL log_runtime('time', 'alloc')
260
261       CALL check(nf90_get_var( ncid, in_var%varid, buffer,                  &
262                                start = start,                                 &
263                                count = count ))
264
265    ELSE
266
267       message = "Failed to read '" // TRIM(in_var%name) // &
268          "' from file '" // TRIM(in_file) // "'."
269       CALL inifor_abort('get_netcdf_variable', message)
270
271    ENDIF
272
273    CALL check(nf90_close(ncid))
274    CALL log_runtime('time', 'read')
275
276 END SUBROUTINE get_netcdf_variable_real
277
278
279!------------------------------------------------------------------------------!
280! Description:
281! ------------
282!> get_netcdf_dim_vector() reads the coordinate array 'coordname' from the
283!> netCDF file 'filename'.
284!------------------------------------------------------------------------------!
285 SUBROUTINE get_netcdf_dim_vector(filename, coordname, coords)
286
287    CHARACTER(LEN=*), INTENT(IN)         ::  filename
288    CHARACTER(LEN=*), INTENT(IN)         ::  coordname
289    REAL(wp), ALLOCATABLE, INTENT(INOUT) ::  coords(:)
290
291    INTEGER ::  ncid, varid, dimlen
292    INTEGER ::  dimids(NF90_MAX_VAR_DIMS)
293
294    IF ( nf90_open( TRIM(filename), NF90_NOWRITE, ncid ) .EQ. NF90_NOERR .AND. &
295         nf90_inq_varid( ncid, coordname, varid ) .EQ. NF90_NOERR )  THEN
296
297       CALL check(nf90_inquire_variable( ncid, varid, dimids = dimids ))
298       CALL check(nf90_inquire_dimension( ncid, dimids(1), len = dimlen ))
299
300       ALLOCATE(coords(dimlen))
301       CALL check(nf90_get_var( ncid, varid, coords))
302
303    ELSE
304
305       message = "Failed to read '" // TRIM(coordname) // &
306          "' from file '" // TRIM(filename) // "'."
307       CALL inifor_abort('get_netcdf_dim_vector', message)
308
309    ENDIF
310
311 END SUBROUTINE get_netcdf_dim_vector
312
313
314!------------------------------------------------------------------------------!
315! Description:
316! ------------
317!> get_input_dimensions() reads dimensions metadata of the netCDF variable given
318!> by 'in_var%name' into 'in_var' data structure.
319!------------------------------------------------------------------------------!
320 SUBROUTINE get_input_dimensions(in_var, ncid)
321
322    TYPE(nc_var), INTENT(INOUT) ::  in_var
323    INTEGER, INTENT(IN)         ::  ncid
324
325    INTEGER ::  i
326
327    CALL check(nf90_get_att( ncid, in_var%varid, "long_name",             &
328                             in_var%long_name))
329
330    CALL check(nf90_get_att( ncid, in_var%varid, "units",                 &
331                             in_var%units ))
332
333    CALL check(nf90_inquire_variable( ncid, in_var%varid,                 &
334                                      ndims  = in_var%ndim,               &
335                                      dimids = in_var%dimids ))
336
337    DO  i = 1, in_var%ndim
338       CALL check(nf90_inquire_dimension( ncid, in_var%dimids(i),         &
339                                          name = in_var%dimname(i),       &
340                                          len  = in_var%dimlen(i) ))
341    ENDDO
342
343 END SUBROUTINE get_input_dimensions
344
345
346!------------------------------------------------------------------------------!
347! Description:
348! ------------
349!> get_netcdf_start_and_count() gets the start position and element counts for
350!> the given netCDF file. This information is used in get_netcdf_variable_int()
351!> and _real() for reading input variables..
352!------------------------------------------------------------------------------!
353 SUBROUTINE get_netcdf_start_and_count(in_var, start, count)
354
355    TYPE(nc_var), INTENT(INOUT)        ::  in_var
356    INTEGER, DIMENSION(3), INTENT(OUT) ::  start, count
357
358    INTEGER ::  ndim
359
360    IF ( in_var%ndim .LT. 2  .OR.  in_var%ndim .GT. 4 )  THEN
361
362       message = "Failed reading NetCDF variable " //                       &
363          TRIM(in_var%name) // " with " //                                  &
364          TRIM(str(INT(in_var%ndim, kind=iwp))) //                          &
365          " dimensions because only two- and and three-dimensional" //      &
366          " variables are supported."
367       CALL inifor_abort('get_netcdf_start_and_count', message)
368
369    ENDIF
370
371    start = (/ 1, 1, 1 /)
372    IF ( TRIM(in_var%name) .EQ. 'T_SO' .AND.                                &
373         in_var%has_redundant_first_level )  THEN
374       
375!
376!--    Skip depth = 0.0 for T_SO and reduce number of depths from 9 to 8
377       in_var%dimlen(3) = in_var%dimlen(3) - 1
378
379!
380!--    Start reading from second level, e.g. depth = 0.005 instead of 0.0
381       start(3) = 2
382    ENDIF
383
384    IF (in_var%ndim .EQ. 2)  THEN
385       in_var%dimlen(3) = 1
386    ENDIF
387
388    ndim = MIN(in_var%ndim, 3)
389    count = (/ 1, 1, 1 /)
390    count(1:ndim) = in_var%dimlen(1:ndim)
391
392 END SUBROUTINE get_netcdf_start_and_count
393
394
395!------------------------------------------------------------------------------!
396! Description:
397! ------------
398!> Routine for defining netCDF variables in the dynamic driver, INIFOR's netCDF
399!> output.
400!------------------------------------------------------------------------------!
401 SUBROUTINE netcdf_define_variable(var, ncid)
402
403     TYPE(nc_var), INTENT(INOUT) ::  var
404     INTEGER, INTENT(IN)         ::  ncid
405
406     CALL check(nf90_def_var(ncid, var%name, NF90_FLOAT,       var%dimids(1:var%ndim), var%varid))
407     CALL check(nf90_put_att(ncid, var%varid, "long_name",     var%long_name))
408     CALL check(nf90_put_att(ncid, var%varid, "units",         var%units))
409     IF ( var%lod .GE. 0 )  THEN
410        CALL check(nf90_put_att(ncid, var%varid, "lod",           var%lod))
411     ENDIF
412     CALL check(nf90_put_att(ncid, var%varid, "source",        var%source))
413     CALL check(nf90_put_att(ncid, var%varid, "_FillValue",    NF90_FILL_REAL))
414
415 END SUBROUTINE netcdf_define_variable
416   
417
418!------------------------------------------------------------------------------!
419! Description:
420! ------------
421!> netcdf_get_dimensions() reads in all dimensions and their lengths and stores
422!> them in the given the 'var' data structure. This information is used later
423!> for writing output variables in update_output().
424!------------------------------------------------------------------------------!
425 SUBROUTINE netcdf_get_dimensions(var, ncid)
426
427     TYPE(nc_var), INTENT(INOUT) ::  var
428     INTEGER, INTENT(IN)         ::  ncid
429     INTEGER                     ::  i
430     CHARACTER(SNAME)            ::  null
431
432     DO  i = 1, var%ndim
433        CALL check(nf90_inquire_dimension(ncid, var%dimids(i), &
434                                          name = null, &
435                                          len  = var%dimlen(i)  ) )
436     ENDDO
437
438 END SUBROUTINE netcdf_get_dimensions
439
440
441!------------------------------------------------------------------------------!
442! Description:
443! ------------
444!> This routine parses and interpretes the command-line options and stores the
445!> resulting settings in the 'cfg' data structure.
446!------------------------------------------------------------------------------!
447 SUBROUTINE parse_command_line_arguments( cfg )
448
449    TYPE(inifor_config), INTENT(INOUT) ::  cfg
450
451    CHARACTER(LEN=PATH)                ::  option, arg
452    INTEGER                            ::  arg_count, i
453
454    cfg%flow_prefix_is_set = .FALSE.
455    cfg%input_prefix_is_set = .FALSE.
456    cfg%p0_is_set = .FALSE.
457    cfg%radiation_prefix_is_set = .FALSE.
458    cfg%soil_prefix_is_set = .FALSE.
459    cfg%precipitation_prefix_is_set = .FALSE.
460    cfg%process_precipitation = .FALSE.
461    cfg%static_driver_is_set = .FALSE.
462    cfg%ug_defined_by_user = .FALSE.
463    cfg%vg_defined_by_user = .FALSE.
464    cfg%z0_is_set = .FALSE.
465
466    arg_count = COMMAND_ARGUMENT_COUNT()
467    IF (arg_count .GT. 0)  THEN
468
469       i = 1
470       DO  WHILE (i .LE. arg_count)
471
472          CALL GET_COMMAND_ARGUMENT( i, option )
473
474          SELECT CASE( TRIM(option) )
475
476             CASE( '--averaging-mode' )
477                CALL get_option_argument( i, arg )
478                cfg%averaging_mode = TRIM(arg)
479
480             CASE( '-date', '-d', '--date' )
481                CALL get_option_argument( i, arg )
482                cfg%start_date = TRIM(arg)
483
484             CASE( '--debug' )
485                cfg%debug = .TRUE.
486
487             CASE( '-z0', '-z', '--elevation' )
488                cfg%z0_is_set = .TRUE.
489                CALL get_option_argument( i, arg )
490                READ(arg, *) cfg%z0
491
492             CASE( '-p0', '-r', '--surface-pressure' )
493                cfg%p0_is_set = .TRUE.
494                CALL get_option_argument( i, arg )
495                READ(arg, *) cfg%p0
496
497             CASE( '-ug', '-u', '--geostrophic-u' )
498                cfg%ug_defined_by_user = .TRUE.
499                CALL get_option_argument( i, arg )
500                READ(arg, *) cfg%ug
501
502             CASE( '-vg', '-v', '--geostrophic-v' )
503                cfg%vg_defined_by_user = .TRUE.
504                CALL get_option_argument( i, arg )
505                READ(arg, *) cfg%vg
506
507             CASE( '-path', '-p', '--path' )
508                CALL get_option_argument( i, arg )
509                 cfg%input_path = TRIM(arg)
510
511             CASE( '-hhl', '-l', '--hhl-file' )
512                CALL get_option_argument( i, arg )
513                cfg%hhl_file = TRIM(arg)
514
515             CASE( '--input-prefix')
516                CALL get_option_argument( i, arg )
517                cfg%input_prefix = TRIM(arg)
518                cfg%input_prefix_is_set = .TRUE.
519   
520             CASE( '-a', '--averaging-angle' )
521                CALL get_option_argument( i, arg )
522                READ(arg, *) cfg%averaging_angle
523
524             CASE( '-static', '-t', '--static-driver' )
525                cfg%static_driver_is_set = .TRUE.
526                CALL get_option_argument( i, arg )
527                cfg%static_driver_file = TRIM(arg)
528
529             CASE( '-soil', '-s', '--soil-file')
530                CALL get_option_argument( i, arg )
531                cfg%soiltyp_file = TRIM(arg)
532
533             CASE( '--flow-prefix')
534                CALL get_option_argument( i, arg )
535                cfg%flow_prefix = TRIM(arg)
536                cfg%flow_prefix_is_set = .TRUE.
537   
538             CASE( '--radiation-prefix')
539                CALL get_option_argument( i, arg )
540                cfg%radiation_prefix = TRIM(arg)
541                cfg%radiation_prefix_is_set = .TRUE.
542   
543             CASE( '--soil-prefix')
544                CALL get_option_argument( i, arg )
545                cfg%soil_prefix = TRIM(arg)
546                cfg%soil_prefix_is_set = .TRUE.
547   
548             CASE( '--precipitation-prefix')
549                CALL get_option_argument( i, arg )
550                cfg%precipitation_prefix = TRIM(arg)
551                cfg%precipitation_prefix_is_set = .TRUE.
552
553             CASE( '--precipitation')
554                cfg%process_precipitation = .TRUE.
555
556             CASE( '-o', '--output' )
557                CALL get_option_argument( i, arg )
558                cfg%output_file = TRIM(arg)
559
560             CASE( '-n', '--namelist' )
561                CALL get_option_argument( i, arg )
562                cfg%namelist_file = TRIM(arg)
563
564             CASE( '-mode', '-i', '--init-mode' )
565                CALL get_option_argument( i, arg )
566                cfg%ic_mode = TRIM(arg)
567
568             CASE( '-f', '--forcing-mode' )
569                CALL get_option_argument( i, arg )
570                cfg%bc_mode = TRIM(arg)
571
572             CASE( '--version' )
573                CALL print_version
574                STOP
575
576             CASE( '--help' )
577                CALL print_version
578                PRINT *, ""
579                PRINT *, &
580                   "For documentation and a list of available command-line options " // NEW_LINE( " " ) // &
581                   " please visit https://palm.muk.uni-hannover.de/trac/wiki/doc/app/iofiles/inifor."
582                STOP
583
584             CASE DEFAULT
585                message = "unknown option '" // TRIM(option) // "'."
586                CALL inifor_abort('parse_command_line_arguments', message)
587
588          END SELECT
589
590          i = i + 1
591
592       ENDDO
593
594    ELSE
595         
596       CALL print_version
597       CALL report( 'parse_command_line_arguments', 'No arguments present, exiting.' )
598       STOP
599
600    ENDIF
601
602 END SUBROUTINE parse_command_line_arguments
603
604
605
606 SUBROUTINE get_datetime_file_list( start_date_string, start_hour, end_hour, &
607                                    step_hour, input_path, prefix, suffix,   &
608                                    file_list )
609
610    CHARACTER (LEN=DATE), INTENT(IN) ::  start_date_string
611    CHARACTER (LEN=*),    INTENT(IN) ::  prefix, suffix, input_path
612    INTEGER(iwp),         INTENT(IN) ::  start_hour, end_hour, step_hour
613    CHARACTER(LEN=*), ALLOCATABLE, INTENT(INOUT) ::  file_list(:)
614
615    INTEGER(iwp)        ::  number_of_intervals, hour, i
616    CHARACTER(LEN=DATE) ::  date_string
617
618    number_of_intervals = CEILING( REAL(end_hour - start_hour) / step_hour )
619    ALLOCATE( file_list(number_of_intervals + 1) )
620
621    DO  i = 0, number_of_intervals
622
623       hour = start_hour + i * step_hour
624       date_string = add_hours_to(start_date_string, hour)
625
626       file_list(i+1) = TRIM(input_path) // TRIM(prefix) //                  &
627                        TRIM(date_string) // TRIM(suffix) // '.nc'
628
629    ENDDO
630
631 END SUBROUTINE get_datetime_file_list
632
633!------------------------------------------------------------------------------!
634! Description:
635! ------------
636!> Establish a list of files based on the given start and end times and file
637!> prefixes and suffixes.
638!------------------------------------------------------------------------------!
639 SUBROUTINE get_input_file_list( start_date_string, start_hour, end_hour,    &
640                                 step_hour, input_path, prefix, suffix,      &
641                                 file_list )
642
643    CHARACTER (LEN=DATE), INTENT(IN) ::  start_date_string
644    CHARACTER (LEN=*),    INTENT(IN) ::  prefix, suffix, input_path
645    INTEGER(iwp),         INTENT(IN) ::  start_hour, end_hour, step_hour
646    CHARACTER(LEN=*), ALLOCATABLE, INTENT(INOUT) ::  file_list(:)
647
648    CALL get_datetime_file_list( start_date_string, start_hour, end_hour,    &
649                                 step_hour, input_path, prefix, suffix,      &
650                                 file_list )
651
652 END SUBROUTINE get_input_file_list
653
654
655!------------------------------------------------------------------------------!
656! Description:
657! ------------
658!> Abort INIFOR if the given file is not present.
659!------------------------------------------------------------------------------!
660LOGICAL FUNCTION file_is_present( filename, file_kind, message )
661
662    CHARACTER(LEN=*), INTENT(IN)            ::  filename, file_kind
663    CHARACTER(LEN=*), INTENT(OUT)           ::  message
664
665    file_is_present = file_present(filename)
666
667    IF (.NOT. file_is_present)  THEN
668
669       IF (LEN( TRIM( filename ) ) == 0)  THEN
670          message = "No name was given for the " // TRIM( file_kind ) // " file."
671       ELSE
672          message = "The " // TRIM( file_kind ) // " file '" //                  &
673                    TRIM( filename ) // "' was not found."
674       ENDIF
675
676    ENDIF
677
678END FUNCTION file_is_present
679
680
681!------------------------------------------------------------------------------!
682! Description:
683! ------------
684!> Abort INIFOR if the given file is not present.
685!------------------------------------------------------------------------------!
686 SUBROUTINE verify_file( file_name, file_kind )
687
688    CHARACTER(LEN=*), INTENT(IN)           ::  file_name, file_kind
689
690    IF (.NOT. file_is_present( file_name, file_kind, message ))  THEN
691
692       CALL inifor_abort( 'verify_file', message )
693
694    ENDIF
695
696    message = "Set up input file name '" // TRIM( file_name ) // "'"
697    CALL report( 'verify_file', message )
698
699 END SUBROUTINE verify_file
700!------------------------------------------------------------------------------!
701! Description:
702! ------------
703!> Get the argument of the i'th command line option, which is at the location
704!> i+1 of the argument list.
705!------------------------------------------------------------------------------!
706 SUBROUTINE get_option_argument( i, arg )
707    CHARACTER(LEN=PATH), INTENT(INOUT) ::  arg
708    INTEGER, INTENT(INOUT)             ::  i
709
710    i = i + 1
711    CALL GET_COMMAND_ARGUMENT( i, arg )
712
713 END SUBROUTINE
714
715
716!------------------------------------------------------------------------------!
717! Description:
718! ------------
719!> Checks the INIFOR configuration 'cfg' for plausibility.
720!------------------------------------------------------------------------------!
721 SUBROUTINE validate_config( cfg )
722    TYPE(inifor_config), INTENT(IN) ::  cfg
723
724    CALL verify_file( cfg%hhl_file, 'HHL' )
725    CALL verify_file( cfg%namelist_file, 'NAMELIST' )
726    CALL verify_file( cfg%soiltyp_file, 'SOILTYP' )
727
728!
729!-- Only check optional static driver file name, if it has been given.
730    IF (TRIM( cfg%static_driver_file ) .NE. '')  THEN
731       CALL verify_file( cfg%static_driver_file, 'static driver' )
732    ENDIF
733
734    SELECT CASE( TRIM( cfg%ic_mode ) )
735       CASE( 'profile', 'volume' )
736       CASE DEFAULT
737          message = "Initialization mode '" // TRIM( cfg%ic_mode ) //&
738                    "' is not supported. " //&
739                    "Please select either 'profile' or 'volume', " //&
740                    "or omit the -i/--init-mode/-mode option entirely, which corresponds "//&
741                    "to the latter."
742          CALL inifor_abort( 'validate_config', message )
743    END SELECT
744
745    SELECT CASE( TRIM(cfg%bc_mode) )
746       CASE( 'real', 'ideal')
747       CASE DEFAULT
748          message = "Forcing mode '" // TRIM( cfg%bc_mode ) //& 
749                    "' is not supported. " //&
750                    "Please select either 'real' or 'ideal', " //&
751                    "or omit the -f/--forcing-mode option entirely, which corresponds "//&
752                    "to the latter."
753          CALL inifor_abort( 'validate_config', message )
754    END SELECT
755
756    SELECT CASE( TRIM( cfg%averaging_mode ) )
757       CASE( 'level' )
758       CASE( 'height' )
759          message = "Averaging mode '" // TRIM( cfg%averaging_mode ) //&
760                    "' is currently not supported. " //&
761                    "Please use level-based averaging by selecting 'level', " //&
762                    "or by omitting the --averaging-mode option entirely."
763          CALL inifor_abort( 'validate_config', message )
764       CASE DEFAULT
765          message = "Averaging mode '" // TRIM( cfg%averaging_mode ) //&
766                    "' is not supported. " //&
767          !          "Please select either 'height' or 'level', " //&
768          !          "or omit the --averaging-mode option entirely, which corresponds "//&
769          !          "to the latter."
770                    "Please use level-based averaging by selecting 'level', " //&
771                    "or by omitting the --averaging-mode option entirely."
772          CALL inifor_abort( 'validate_config', message )
773    END SELECT
774
775    IF ( cfg%ug_defined_by_user .NEQV. cfg%vg_defined_by_user )  THEN
776       message = "You specified only one component of the geostrophic " // &
777                 "wind. Please specify either both or none."
778       CALL inifor_abort( 'validate_config', message )
779    ENDIF
780
781    IF ( .NOT. cfg%static_driver_is_set .AND. .NOT. cfg%z0_is_set )  THEN
782       message =                                                               &
783          "The vertical origin of the PALM grid has not been defined. " // NEW_LINE( " " ) // &
784          "Please specify the right value for your setup by either " // NEW_LINE( " " ) // &
785          "  - using the command-line option --elevation <height above sea level>, or by" // NEW_LINE( " " ) // &
786          "  - specifying a static driver file using --static <filename> in order to use " // NEW_LINE( " " ) // &
787          "    use the value of origin_z (and origin_lon and origin_lat) specifed therein."
788       CALL inifor_abort( 'validate_config', message )
789    ENDIF
790
791 END SUBROUTINE validate_config
792
793
794 SUBROUTINE get_cosmo_grid( hhl_file, soil_file, rlon, rlat, hhl, hfl, depths, &
795                            d_depth, d_depth_rho_inv, phi_n, lambda_n,       &
796                            phi_equat,                                       &
797                            lonmin_cosmo, lonmax_cosmo,                      &
798                            latmin_cosmo, latmax_cosmo,                      &
799                            nlon, nlat, nlev, ndepths )
800
801    CHARACTER(LEN=PATH), INTENT(IN)                      ::  hhl_file  !< path to file containing the HHL variable (height of half layers)
802    CHARACTER(LEN=PATH), INTENT(IN)                      ::  soil_file !< path to one of the soil input files for reading soil layer depths
803    REAL(wp), DIMENSION(:), ALLOCATABLE, INTENT(OUT)     ::  rlon      !< longitudes of COSMO-DE's rotated-pole grid
804    REAL(wp), DIMENSION(:), ALLOCATABLE, INTENT(OUT)     ::  rlat      !< latitudes of COSMO-DE's rotated-pole grid
805    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, INTENT(OUT) ::  hhl       !< heights of half layers (cell faces) above sea level in COSMO-DE, read in from external file
806    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, INTENT(OUT) ::  hfl       !< heights of full layers (cell centres) above sea level in COSMO-DE, computed as arithmetic average of hhl
807    REAL(wp), DIMENSION(:), ALLOCATABLE, INTENT(OUT)     ::  depths    !< COSMO-DE's TERRA-ML soil layer depths
808    REAL(wp), DIMENSION(:), ALLOCATABLE, INTENT(OUT)     ::  d_depth
809    REAL(wp), DIMENSION(:), ALLOCATABLE, INTENT(OUT)     ::  d_depth_rho_inv
810    REAL(wp), INTENT(OUT)                                ::  phi_n
811    REAL(wp), INTENT(OUT)                                ::  phi_equat
812    REAL(wp), INTENT(OUT)                                ::  lambda_n
813    REAL(wp), INTENT(OUT)                                ::  lonmin_cosmo !< Minimunm longitude of COSMO-DE's rotated-pole grid [COSMO rotated-pole rad]
814    REAL(wp), INTENT(OUT)                                ::  lonmax_cosmo !< Maximum longitude of COSMO-DE's rotated-pole grid [COSMO rotated-pole rad]
815    REAL(wp), INTENT(OUT)                                ::  latmin_cosmo !< Minimunm latitude of COSMO-DE's rotated-pole grid [COSMO rotated-pole rad]
816    REAL(wp), INTENT(OUT)                                ::  latmax_cosmo !< Maximum latitude of COSMO-DE's rotated-pole grid [COSMO rotated-pole rad]
817    INTEGER(iwp), INTENT(OUT)                            ::  nlon, nlat, nlev, ndepths
818
819    TYPE(nc_var) ::  cosmo_var !< COSMO dummy variable, used for reading HHL, rlon, rlat
820    INTEGER(iwp) ::  k
821
822!
823!-- Read in COSMO's heights of half layers (vertical cell faces)
824    cosmo_var%name = NC_HHL_NAME
825    CALL get_netcdf_variable( hhl_file, cosmo_var, hhl )
826    CALL get_netcdf_dim_vector( hhl_file, NC_RLON_NAME, rlon )
827    CALL get_netcdf_dim_vector( hhl_file, NC_RLAT_NAME, rlat )
828    CALL get_netcdf_dim_vector( soil_file, NC_DEPTH_NAME, depths)
829    CALL log_runtime( 'time', 'read' )
830
831    CALL reverse( hhl )
832    nlon = SIZE( hhl, 1 )
833    nlat = SIZE( hhl, 2 )
834    nlev = SIZE( hhl, 3 )
835    ndepths = SIZE( depths )
836
837    CALL log_runtime( 'time', 'comp' )
838
839    ALLOCATE( hfl( nlon, nlat, nlev-1 ) )
840    ALLOCATE( d_depth( ndepths ), d_depth_rho_inv( ndepths ) )
841    CALL log_runtime('time', 'alloc')
842
843    CALL get_soil_layer_thickness( depths, d_depth )
844    d_depth_rho_inv = 1.0_wp / ( d_depth * RHO_L )
845
846!
847!-- Compute COSMO's heights of full layers (cell centres)
848    DO  k = 1, nlev-1
849       hfl(:,:,k) = 0.5_wp * ( hhl(:,:,k) +                                  &
850                               hhl(:,:,k+1) )
851    ENDDO
852!
853!-- COSMO rotated pole coordinates
854    phi_n = TO_RADIANS                                                       &
855          * get_netcdf_variable_attribute( hhl_file,                         &
856                                           NC_ROTATED_POLE_NAME,             &
857                                           NC_POLE_LATITUDE_NAME )
858
859    lambda_n = TO_RADIANS                                                    &
860             * get_netcdf_variable_attribute( hhl_file,                      &
861                                              NC_ROTATED_POLE_NAME,          &
862                                              NC_POLE_LONGITUDE_NAME )
863
864    phi_equat = 90.0_wp * TO_RADIANS - phi_n
865
866    lonmin_cosmo = MINVAL( rlon ) * TO_RADIANS
867    lonmax_cosmo = MAXVAL( rlon ) * TO_RADIANS
868    latmin_cosmo = MINVAL( rlat ) * TO_RADIANS
869    latmax_cosmo = MAXVAL( rlat ) * TO_RADIANS
870    CALL log_runtime('time', 'comp')
871
872 END SUBROUTINE get_cosmo_grid
873
874
875!------------------------------------------------------------------------------!
876! Description:
877! ------------
878!> Fills the thickness array of the COSMO soil layers. Since COSMO's (i.e.
879!> TERRA_ML's [1]) soil layer boundaries follow the rule
880!>
881!>    depth(0) = 0.0, and
882!>    depth(k) = 0.01 * 3**(k-1), k in [1,2,3,...,7]
883!>
884!> and full levels are defined as the midpoints between two layer boundaries,
885!> all except the first layer thicknesses equal the depth of the midpoint.
886!>
887!> [1] A Description of the Nonhydrostatic Regional COSMO Model Part II :
888!>     Physical Parameterization*, Sect. 11 TERRA_ML.
889!>     http://www.cosmo-model.org/content/model/documentation/core/cosmoPhysParamtr.pdf)
890!>
891!> Input parameters:
892!> -----------------
893!>
894!> depths: array of full soil layer depths (cell midpoints)
895!>
896!>
897!> Output parameters:
898!> ------------------
899!>
900!> d_depth: array of soil layer thicknesses
901!>
902!------------------------------------------------------------------------------!
903 SUBROUTINE get_soil_layer_thickness( depths, d_depth )
904
905    REAL(wp), INTENT(IN)  ::  depths(:)
906    REAL(wp), INTENT(OUT) ::  d_depth(:)
907
908    d_depth(:) = depths(:)
909    d_depth(1) = 2.0_wp * depths(1)
910
911 END SUBROUTINE get_soil_layer_thickness
912!------------------------------------------------------------------------------!
913! Description:
914! ------------
915!> Check whether the given file is present on the filesystem.
916!------------------------------------------------------------------------------!
917 LOGICAL FUNCTION file_present( filename )
918    CHARACTER(LEN=PATH), INTENT(IN) ::  filename
919
920    INQUIRE( FILE=filename, EXIST=file_present )
921
922 END FUNCTION file_present
923
924
925!------------------------------------------------------------------------------!
926! Description:
927! ------------
928!> This routine initializes the dynamic driver file, i.e. INIFOR's netCDF output
929!> file.
930!>
931!> Besides writing metadata, such as global attributes, coordinates, variables,
932!> in the netCDF file, various netCDF IDs are saved for later, when INIFOR
933!> writes the actual data.
934!------------------------------------------------------------------------------!
935 SUBROUTINE setup_netcdf_dimensions( output_file, palm_grid,                  &
936                                     start_date_string, origin_lon, origin_lat )
937
938    TYPE(nc_file), INTENT(INOUT)      ::  output_file
939    TYPE(grid_definition), INTENT(IN) ::  palm_grid
940    CHARACTER (LEN=DATE), INTENT(IN)  ::  start_date_string
941    REAL(wp), INTENT(IN)              ::  origin_lon, origin_lat
942
943    CHARACTER (LEN=8)     ::  date_string
944    CHARACTER (LEN=10)    ::  time_string
945    CHARACTER (LEN=5)     ::  zone_string
946    CHARACTER (LEN=SNAME) ::  history_string
947    INTEGER               ::  nx, ny, nz, nt
948    INTEGER               ::  ncid, dimids(3), dimvarids(3)
949    REAL(wp)              ::  z0
950
951    message = "Initializing PALM-4U dynamic driver file '" //               &
952              TRIM(output_file%name) // "' and setting up dimensions."
953    CALL report('setup_netcdf_dimensions', message)
954
955!
956!-- Create the netCDF file as in netCDF-4/HDF5 format if __netcdf4 preprocessor flag is given
957#if defined( __netcdf4 )
958    CALL check(nf90_create(TRIM(output_file%name), OR(NF90_CLOBBER, NF90_HDF5), ncid))
959#else
960    CALL check(nf90_create(TRIM(output_file%name), NF90_CLOBBER, ncid))
961#endif
962
963!------------------------------------------------------------------------------
964!- Section 1: Define NetCDF dimensions and coordinates
965!------------------------------------------------------------------------------
966    nt = SIZE(output_file%time)
967    nx = palm_grid%nx
968    ny = palm_grid%ny
969    nz = palm_grid%nz
970    z0 = palm_grid%z0
971
972
973!
974!------------------------------------------------------------------------------
975!- Section 2: Write global NetCDF attributes
976!------------------------------------------------------------------------------
977    CALL date_and_time(DATE=date_string, TIME=time_string, ZONE=zone_string)
978    history_string =                                                        &
979        'Created on '// date_string      //                                 &
980        ' at '       // time_string(1:2) // ':' // time_string(3:4) //      &
981        ' (UTC'      // zone_string // ')'
982
983    CALL check(nf90_put_att(ncid, NF90_GLOBAL, 'title',          'PALM input file for scenario ...'))
984    CALL check(nf90_put_att(ncid, NF90_GLOBAL, 'institution',    'Deutscher Wetterdienst, Offenbach'))
985    CALL check(nf90_put_att(ncid, NF90_GLOBAL, 'author',         'Eckhard Kadasch, eckhard.kadasch@dwd.de'))
986    CALL check(nf90_put_att(ncid, NF90_GLOBAL, 'history',        TRIM(history_string)))
987    CALL check(nf90_put_att(ncid, NF90_GLOBAL, 'references',     '--'))
988    CALL check(nf90_put_att(ncid, NF90_GLOBAL, 'comment',        '--'))
989    CALL check(nf90_put_att(ncid, NF90_GLOBAL, 'origin_lat',     TRIM(real_to_str(origin_lat*TO_DEGREES, '(F18.13)'))))
990    CALL check(nf90_put_att(ncid, NF90_GLOBAL, 'origin_lon',     TRIM(real_to_str(origin_lon*TO_DEGREES, '(F18.13)'))))
991!
992!-- FIXME: This is the elevation relative to COSMO-DE/D2 sea level and does
993!-- FIXME: not necessarily comply with DHHN2016 (c.f. PALM Input Data
994!-- FIXME: Standard v1.9., origin_z)
995    CALL check(nf90_put_att(ncid, NF90_GLOBAL, 'origin_z',       TRIM(real_to_str(z0, '(F18.13)'))))
996    CALL check(nf90_put_att(ncid, NF90_GLOBAL, 'inifor_version', TRIM(VERSION)))
997    CALL check(nf90_put_att(ncid, NF90_GLOBAL, 'palm_version',   '--'))
998
999!
1000!
1001!------------------------------------------------------------------------------
1002!- Section 2a: Define dimensions for cell centers (scalars in soil and atmosph.)
1003!------------------------------------------------------------------------------
1004!
1005!-- reset dimids first
1006    dimids = (/0, 0, 0/)
1007    CALL check( nf90_def_dim(ncid, "x", nx+1, dimids(1)) )
1008    CALL check( nf90_def_dim(ncid, "y", ny+1, dimids(2)) )
1009    CALL check( nf90_def_dim(ncid, "z", nz, dimids(3)) )
1010!
1011!-- save dimids for later
1012    output_file%dimids_scl = dimids 
1013
1014!
1015!-- reset dimvarids first
1016    dimvarids = (/0, 0, 0/)
1017    CALL check(nf90_def_var(ncid, "x", NF90_FLOAT, dimids(1), dimvarids(1)))
1018    CALL check(nf90_put_att(ncid, dimvarids(1), "standard_name", "x coordinate of cell centers"))
1019    CALL check(nf90_put_att(ncid, dimvarids(1), "units", "m"))
1020
1021    CALL check(nf90_def_var(ncid, "y", NF90_FLOAT, dimids(2), dimvarids(2)))
1022    CALL check(nf90_put_att(ncid, dimvarids(2), "standard_name", "y coordinate of cell centers"))
1023    CALL check(nf90_put_att(ncid, dimvarids(2), "units", "m"))
1024
1025    CALL check(nf90_def_var(ncid, "z", NF90_FLOAT, dimids(3), dimvarids(3)))
1026    CALL check(nf90_put_att(ncid, dimvarids(3), "standard_name", "z coordinate of cell centers"))
1027    CALL check(nf90_put_att(ncid, dimvarids(3), "units", "m"))
1028!
1029!-- save dimvarids for later
1030    output_file%dimvarids_scl = dimvarids
1031
1032!
1033!-- overwrite third dimid with the one of depth
1034    CALL check(nf90_def_dim(ncid, "zsoil", SIZE(palm_grid%depths), dimids(3)) )
1035!
1036!-- save dimids for later
1037    output_file%dimids_soil = dimids
1038
1039!
1040!-- overwrite third dimvarid with the one of depth
1041    CALL check(nf90_def_var(ncid, "zsoil", NF90_FLOAT, output_file%dimids_soil(3), dimvarids(3)))
1042    CALL check(nf90_put_att(ncid, dimvarids(3), "standard_name", "depth_below_land"))
1043    CALL check(nf90_put_att(ncid, dimvarids(3), "positive", "down"))
1044    CALL check(nf90_put_att(ncid, dimvarids(3), "units", "m"))
1045!
1046!-- save dimvarids for later
1047    output_file%dimvarids_soil = dimvarids
1048!
1049!------------------------------------------------------------------------------
1050!- Section 2b: Define dimensions for cell faces/velocities
1051!------------------------------------------------------------------------------
1052!
1053!-- reset dimids first
1054    dimids = (/0, 0, 0/)
1055    CALL check(nf90_def_dim(ncid, "xu", nx, dimids(1)) )
1056    CALL check(nf90_def_dim(ncid, "yv", ny, dimids(2)) )
1057    CALL check(nf90_def_dim(ncid, "zw", nz-1, dimids(3)) )
1058!
1059!-- save dimids for later
1060    output_file%dimids_vel = dimids
1061
1062!
1063!-- reset dimvarids first
1064    dimvarids = (/0, 0, 0/)
1065    CALL check(nf90_def_var(ncid, "xu", NF90_FLOAT, dimids(1), dimvarids(1)))
1066    CALL check(nf90_put_att(ncid, dimvarids(1), "standard_name", "x coordinate of cell faces"))
1067    CALL check(nf90_put_att(ncid, dimvarids(1), "units", "m"))
1068
1069    CALL check(nf90_def_var(ncid, "yv", NF90_FLOAT, dimids(2), dimvarids(2)))
1070    CALL check(nf90_put_att(ncid, dimvarids(2), "standard_name", "y coordinate of cell faces"))
1071    CALL check(nf90_put_att(ncid, dimvarids(2), "units", "m"))
1072
1073    CALL check(nf90_def_var(ncid, "zw", NF90_FLOAT, dimids(3), dimvarids(3)))
1074    CALL check(nf90_put_att(ncid, dimvarids(3), "standard_name", "z coordinate of cell faces"))
1075    CALL check(nf90_put_att(ncid, dimvarids(3), "units", "m"))
1076!
1077!-- save dimvarids for later
1078    output_file%dimvarids_vel = dimvarids
1079
1080!
1081!------------------------------------------------------------------------------
1082!- Section 2c: Define time dimension
1083!------------------------------------------------------------------------------
1084    CALL check(nf90_def_dim(ncid, "time", nt, output_file%dimid_time) )
1085    CALL check(nf90_def_var(ncid, "time", NF90_FLOAT, &
1086                                          output_file%dimid_time, &
1087                                          output_file%dimvarid_time))
1088    CALL check(nf90_put_att(ncid, output_file%dimvarid_time, "standard_name", "time"))
1089    CALL check(nf90_put_att(ncid, output_file%dimvarid_time, "long_name", "time"))
1090    CALL check(nf90_put_att(ncid, output_file%dimvarid_time, "units",     &
1091                            "seconds since " // start_date_string // " UTC"))
1092
1093    CALL check(nf90_enddef(ncid))
1094
1095!
1096!------------------------------------------------------------------------------
1097!- Section 3: Write grid coordinates
1098!------------------------------------------------------------------------------
1099    CALL check(nf90_put_var(ncid, output_file%dimvarids_scl(1), palm_grid%x))
1100    CALL check(nf90_put_var(ncid, output_file%dimvarids_scl(2), palm_grid%y))
1101    CALL check(nf90_put_var(ncid, output_file%dimvarids_scl(3), palm_grid%z))
1102
1103    CALL check(nf90_put_var(ncid, output_file%dimvarids_vel(1), palm_grid%xu))
1104    CALL check(nf90_put_var(ncid, output_file%dimvarids_vel(2), palm_grid%yv))
1105    CALL check(nf90_put_var(ncid, output_file%dimvarids_vel(3), palm_grid%zw))
1106
1107!
1108!-- TODO Read in soil depths from input file before this.
1109    CALL check(nf90_put_var(ncid, output_file%dimvarids_soil(3), palm_grid%depths))
1110
1111!
1112!-- Write time vector
1113    CALL check(nf90_put_var(ncid, output_file%dimvarid_time, output_file%time))
1114
1115!
1116!-- Close the file
1117    CALL check(nf90_close(ncid))
1118
1119 END SUBROUTINE setup_netcdf_dimensions
1120
1121
1122!------------------------------------------------------------------------------!
1123! Description:
1124! ------------
1125!> Defines the netCDF variables to be written to the dynamic driver file
1126!------------------------------------------------------------------------------!
1127 SUBROUTINE setup_netcdf_variables(filename, io_group_list)
1128
1129    CHARACTER (LEN=*), INTENT(IN)       ::  filename
1130    TYPE(io_group), INTENT(IN), TARGET  ::  io_group_list(:)
1131
1132    TYPE(io_group), POINTER             ::  group
1133    TYPE(nc_var), POINTER               ::  var
1134    INTEGER(iwp)                        ::  group_idx, var_idx, n_var
1135    INTEGER                             ::  ncid
1136    LOGICAL                             ::  to_be_written
1137
1138    message = "Defining variables in dynamic driver '" // TRIM(filename) // "'."
1139    CALL report('setup_netcdf_variables', message)
1140
1141    CALL check(nf90_open(TRIM(filename), NF90_WRITE, ncid))
1142    CALL check(nf90_redef(ncid))
1143
1144    n_var = 0
1145    DO  group_idx = 1, SIZE( io_group_list )
1146
1147       group => io_group_list(group_idx)
1148       DO var_idx = 1, SIZE( group%out_vars )
1149
1150          to_be_written = .FALSE.
1151
1152          IF (ALLOCATED( group%out_vars ))  THEN
1153             var => group%out_vars(var_idx)
1154             n_var = n_var + 1
1155
1156             to_be_written = (                                                 &
1157                group%to_be_processed  .AND.  var%to_be_processed              &
1158                .AND. .NOT. var%is_internal                                    &
1159             )
1160          ENDIF
1161
1162          IF ( to_be_written )  THEN
1163             message = "  variable #" // TRIM(str(n_var)) // " '" // TRIM(var%name) // "'."
1164             CALL report('setup_netcdf_variables', message)
1165
1166             CALL netcdf_define_variable(var, ncid)
1167             CALL netcdf_get_dimensions(var, ncid)
1168          ENDIF
1169
1170       ENDDO
1171       
1172    ENDDO
1173
1174    CALL check(nf90_enddef(ncid))
1175    CALL check(nf90_close(ncid))
1176
1177    message = "Dynamic driver '" // TRIM(filename) // "' initialized successfully."
1178    CALL report('setup_netcdf_variables', message)
1179
1180 END SUBROUTINE setup_netcdf_variables
1181
1182
1183!------------------------------------------------------------------------------!
1184! Description:
1185! ------------
1186!> This routine reads and returns all input variables of the given IO group
1187!> It accomodates the data by allocating a container variable that represents a
1188!> list of arrays of the same length as the groups variable list. (This list
1189!> will typically contain one or two items.) After the container, its members
1190!> are allocated one by one with the appropriate, possibly different,
1191!> dimensions.
1192!>
1193!> The 'group' is an INTENT(INOUT) variable so that 'get_netcdf_variable()' can
1194!> record netCDF IDs in the 'in_var_list()' member variable.
1195!------------------------------------------------------------------------------!
1196 SUBROUTINE read_input_variables(group, iter, buffer)
1197    TYPE(io_group), INTENT(INOUT), TARGET       ::  group
1198    INTEGER(iwp), INTENT(IN)                    ::  iter
1199    TYPE(container), ALLOCATABLE, INTENT(INOUT) ::  buffer(:)
1200    INTEGER(iwp)                                ::  hour, buf_id
1201    TYPE(nc_var), POINTER                       ::  input_var
1202    CHARACTER(LEN=PATH), POINTER                ::  input_file
1203    INTEGER(iwp)                                ::  ivar, nbuffers
1204
1205    message = "Reading data for I/O group '" // TRIM(group%in_var_list(1)%name) // "'."
1206    CALL report('read_input_variables', message)
1207
1208    input_file => group%in_files(iter)
1209
1210!
1211!------------------------------------------------------------------------------
1212!- Section 1: Load input buffers for accumulated variables
1213!------------------------------------------------------------------------------
1214!
1215!-- radiation budgets, precipitation
1216    IF ( group%kind == 'running average' .OR.                                  &
1217         group%kind == 'accumulated' )  THEN
1218
1219       IF ( SIZE( group%in_var_list ) .GT. 1 ) THEN
1220          message = "I/O groups may not contain more than one " // & 
1221                    "accumulated variable. Group '" // TRIM(group%kind) //&
1222                    "' contains " //                                           &
1223                    TRIM( str( SIZE( group%in_var_list, kind=iwp ) ) ) // "."
1224          CALL inifor_abort( 'read_input_variables | accumulation', message )
1225       ENDIF
1226
1227!
1228!--    use two buffer arrays
1229       nbuffers = 2
1230       IF ( .NOT. ALLOCATED( buffer ) )  ALLOCATE( buffer(nbuffers) )
1231
1232!
1233!--    hour of the day
1234       hour = iter - 1
1235!
1236!--    chose correct buffer array
1237       buf_id = select_buffer(hour)
1238
1239       CALL log_runtime('time', 'read')
1240       IF ( ALLOCATED(buffer(buf_id)%array) )  DEALLOCATE(buffer(buf_id)%array)
1241       CALL log_runtime('time', 'alloc')
1242
1243       input_var => group%in_var_list(1)
1244       CALL get_netcdf_variable(input_file, input_var, buffer(buf_id)%array)
1245       CALL report('read_input_variables', "Read accumulated " // TRIM(group%in_var_list(1)%name)) 
1246
1247       IF ( input_var%is_upside_down )  CALL reverse(buffer(buf_id)%array)
1248       CALL log_runtime('time', 'comp')
1249         
1250!------------------------------------------------------------------------------
1251!- Section 2: Load input buffers for normal I/O groups
1252!------------------------------------------------------------------------------
1253    ELSE
1254
1255!
1256!--    Allocate one input buffer per output quantity. If more quantities
1257!--    have to be computed than input variables exist in this group,
1258!--    allocate more buffers. For instance, in the thermodynamics group,
1259!--    there are three input variabels (PP, T, Qv) and four quantities
1260!--    necessart (P, Theta, Rho, qv) for the corresponding output fields
1261!--    (p0, Theta, qv, ug, and vg)
1262       ALLOCATE( buffer(group%n_output_quantities) )
1263       CALL log_runtime('time', 'alloc')
1264       
1265!
1266!--    Read in all input variables, leave extra buffers-if any-untouched.
1267       DO  ivar = 1, group%n_inputs
1268
1269          input_var => group%in_var_list(ivar)
1270
1271          IF ( input_var%to_be_processed )  THEN
1272!            Check wheather P or PP is present in input file
1273             IF (input_var%name == 'P')  THEN
1274                input_var%name = TRIM( get_pressure_varname(input_file) )
1275             CALL log_runtime('time', 'read')
1276             ENDIF
1277
1278             CALL get_netcdf_variable(input_file, input_var, buffer(ivar)%array)
1279
1280             IF ( input_var%is_upside_down )  CALL reverse(buffer(ivar)%array)
1281             CALL log_runtime('time', 'comp')
1282          ENDIF
1283
1284       ENDDO
1285    ENDIF
1286
1287 END SUBROUTINE read_input_variables
1288
1289
1290!------------------------------------------------------------------------------!
1291! Description:
1292! ------------
1293!> Select the appropriate buffer ID for accumulated COSMO input variables
1294!> depending on the current hour.
1295!------------------------------------------------------------------------------!
1296 INTEGER(iwp) FUNCTION select_buffer(hour)
1297    INTEGER(iwp), INTENT(IN) ::  hour
1298    INTEGER(iwp)             ::  step
1299
1300    select_buffer = 0_iwp
1301    step = MODULO(hour, 3_iwp) + 1_iwp
1302
1303    SELECT CASE(step)
1304       CASE(1, 3)
1305           select_buffer = 1_iwp
1306       CASE(2)
1307           select_buffer = 2_iwp
1308       CASE DEFAULT
1309           message = "Invalid step '" // TRIM(str(step))
1310           CALL inifor_abort('select_buffer', message)
1311    END SELECT
1312 END FUNCTION select_buffer
1313
1314
1315!------------------------------------------------------------------------------!
1316! Description:
1317! ------------
1318!> Checks if the input_file contains the absolute pressure, 'P', or the pressure
1319!> perturbation, 'PP', and returns the appropriate string.
1320!------------------------------------------------------------------------------!
1321 CHARACTER(LEN=2) FUNCTION get_pressure_varname(input_file) RESULT(var)
1322    CHARACTER(LEN=*) ::  input_file
1323    INTEGER          ::  ncid, varid
1324
1325    CALL check(nf90_open( TRIM(input_file), NF90_NOWRITE, ncid ))
1326    IF ( nf90_inq_varid( ncid, 'P', varid ) .EQ. NF90_NOERR )  THEN
1327
1328       var = 'P'
1329
1330    ELSE IF ( nf90_inq_varid( ncid, 'PP', varid ) .EQ. NF90_NOERR )  THEN
1331
1332       var = 'PP'
1333       CALL report('get_pressure_var', 'Using PP instead of P')
1334
1335    ELSE
1336
1337       message = "Failed to read '" // TRIM(var) // &
1338                 "' from file '" // TRIM(input_file) // "'."
1339       CALL inifor_abort('get_pressure_var', message)
1340
1341    ENDIF
1342
1343    CALL check(nf90_close(ncid))
1344
1345 END FUNCTION get_pressure_varname
1346
1347
1348!------------------------------------------------------------------------------!
1349! Description:
1350! ------------
1351!> Read the given global attribute form the given netCDF file.
1352!------------------------------------------------------------------------------!
1353 FUNCTION get_netcdf_attribute(filename, attribute) RESULT(attribute_value)
1354
1355    CHARACTER(LEN=*), INTENT(IN) ::  filename, attribute
1356    REAL(wp)                     ::  attribute_value
1357
1358    INTEGER                      ::  ncid
1359
1360    IF ( nf90_open( TRIM(filename), NF90_NOWRITE, ncid ) == NF90_NOERR )  THEN
1361
1362       CALL check(nf90_get_att(ncid, NF90_GLOBAL, TRIM(attribute), attribute_value))
1363       CALL check(nf90_close(ncid))
1364
1365    ELSE
1366
1367       message = "Failed to read '" // TRIM(attribute) // &
1368                 "' from file '" // TRIM(filename) // "'."
1369       CALL inifor_abort('get_netcdf_attribute', message)
1370
1371    ENDIF
1372
1373 END FUNCTION get_netcdf_attribute
1374
1375
1376!------------------------------------------------------------------------------!
1377! Description:
1378! ------------
1379!> Read the attribute of the given variable form the given netCDF file.
1380!------------------------------------------------------------------------------!
1381 FUNCTION get_netcdf_variable_attribute(filename, varname, attribute)       &
1382    RESULT(attribute_value)
1383
1384    CHARACTER(LEN=*), INTENT(IN) ::  filename, varname, attribute
1385    REAL(wp)                     ::  attribute_value
1386
1387    INTEGER                      ::  ncid, varid
1388
1389    IF ( nf90_open( TRIM(filename), NF90_NOWRITE, ncid ) == NF90_NOERR )  THEN
1390
1391       CALL check( nf90_inq_varid( ncid, TRIM( varname ), varid ) )
1392       CALL check( nf90_get_att( ncid, varid, TRIM( attribute ),            &
1393                   attribute_value ) )
1394       CALL check( nf90_close( ncid ) )
1395
1396    ELSE
1397
1398       message = "Failed to read '" // TRIM( varname ) // ":" //            &
1399                 TRIM( attribute ) // "' from file '" // TRIM(filename) // "'."
1400       CALL inifor_abort('get_netcdf_variable_attribute', message)
1401
1402    ENDIF
1403
1404 END FUNCTION get_netcdf_variable_attribute
1405
1406!------------------------------------------------------------------------------!
1407! Description:
1408! ------------
1409!> Updates the dynamic driver with the interpolated field of the current
1410!> variable at the current time step.
1411!------------------------------------------------------------------------------!
1412 SUBROUTINE update_output(var, array, iter, output_file, cfg)
1413    TYPE(nc_var), INTENT(IN)  ::  var
1414    REAL(wp), INTENT(IN)      ::  array(:,:,:)
1415    INTEGER(iwp), INTENT(IN)  ::  iter
1416    TYPE(nc_file), INTENT(IN) ::  output_file
1417    TYPE(inifor_config)       ::  cfg
1418
1419    INTEGER      ::  ncid, ndim, start(4), count(4)
1420    LOGICAL      ::  var_is_time_dependent
1421
1422    var_is_time_dependent = (                                               &
1423       var%dimids( var%ndim ) == output_file%dimid_time               &
1424    )
1425
1426!
1427!-- Skip time dimension for output
1428    ndim = var%ndim
1429    IF ( var_is_time_dependent )  ndim = var%ndim - 1
1430
1431    start(:)      = (/1,1,1,1/)
1432    start(ndim+1) = iter
1433    count(1:ndim) = var%dimlen(1:ndim)
1434
1435    CALL check(nf90_open(output_file%name, NF90_WRITE, ncid))
1436
1437!
1438!-- Reduce dimension of output array according to variable kind
1439    SELECT CASE (TRIM(var%kind))
1440       
1441       CASE ( 'init scalar profile', 'init u profile', 'init v profile',       &
1442              'init w profile' )
1443
1444          CALL check(nf90_put_var( ncid, var%varid, array(1,1,:) ) )
1445
1446       CASE ( 'init soil', 'init scalar', 'init u', 'init v', 'init w' )
1447
1448          CALL check(nf90_put_var( ncid, var%varid, array(:,:,:) ) )
1449
1450       CASE( 'left scalar', 'right scalar', 'left w', 'right w' ) 
1451
1452          CALL check(nf90_put_var( ncid, var%varid, array(1,:,:),              &
1453                                   start=start(1:ndim+1),                      &
1454                                   count=count(1:ndim) ) )
1455         
1456
1457          IF (.NOT. SIZE(array, 2) .EQ. var%dimlen(1))  THEN
1458             PRINT *, "inifor: update_output: Dimension ", 1, " of variable ", &
1459                 TRIM(var%name), " (", var%dimlen(1),                      &
1460                 ") does not match the dimension of the output array (",       &
1461                 SIZE(array, 2), ")."
1462             STOP
1463          ENDIF
1464         
1465
1466       CASE( 'north scalar', 'south scalar', 'north w', 'south w' )
1467
1468          CALL check(nf90_put_var( ncid, var%varid, array(:,1,:),              &
1469                                   start=start(1:ndim+1),                      &
1470                                   count=count(1:ndim) ) )
1471         
1472
1473       CASE( 'surface forcing', 'top scalar', 'top w' )
1474
1475          CALL check(nf90_put_var( ncid, var%varid, array(:,:,1),              &
1476                                   start=start(1:ndim+1),                      &
1477                                   count=count(1:ndim) ) )
1478         
1479       CASE ( 'left u', 'right u', 'left v', 'right v' )
1480
1481          CALL check(nf90_put_var( ncid, var%varid, array(1,:,:),              &
1482                                   start=start(1:ndim+1),                      &
1483                                   count=count(1:ndim) ) )
1484
1485       CASE ( 'north u', 'south u', 'north v', 'south v' )
1486
1487          CALL check(nf90_put_var( ncid, var%varid, array(:,1,:),              &
1488                                   start=start(1:ndim+1),                      &
1489                                   count=count(1:ndim) ) )
1490
1491       CASE ( 'top u', 'top v' )
1492
1493          CALL check(nf90_put_var( ncid, var%varid, array(:,:,1),              &
1494                                   start=start(1:ndim+1),                      &
1495                                   count=count(1:ndim) ) )
1496
1497       CASE ( 'time series' )
1498
1499          CALL check(nf90_put_var( ncid, var%varid, array(1,1,1),              &
1500                                   start=start(1:ndim+1) ) )
1501
1502       CASE ( 'constant scalar profile', 'geostrophic' )
1503
1504          CALL check(nf90_put_var( ncid, var%varid, array(1,1,:),              &
1505                                   start=start(1:ndim+1),                      &
1506                                   count=count(1:ndim) ) )
1507
1508       CASE ( 'internal profile' )
1509
1510          IF ( cfg%debug )  THEN
1511             CALL check(nf90_put_var( ncid, var%varid, array(1,1,:),           &
1512                                      start=start(1:ndim+1),                   &
1513                                      count=count(1:ndim) ) )
1514          ENDIF
1515
1516       CASE ( 'large-scale scalar forcing', 'large-scale w forcing' )
1517
1518           message = "Doing nothing in terms of writing large-scale forings."
1519           CALL report('update_output', message)
1520
1521       CASE DEFAULT
1522
1523           message = "Variable kind '" // TRIM(var%kind) //                  &
1524                    "' not recognized."
1525           CALL inifor_abort('update_output', message)
1526
1527    END SELECT
1528
1529    CALL check(nf90_close(ncid))
1530
1531 END SUBROUTINE update_output
1532
1533
1534!------------------------------------------------------------------------------!
1535! Description:
1536! ------------
1537!> Checks the status of a netCDF API call and aborts if an error occured
1538!------------------------------------------------------------------------------!
1539 SUBROUTINE check(status)
1540
1541    INTEGER, INTENT(IN) ::  status
1542
1543    IF (status /= nf90_noerr)  THEN
1544       message = "NetCDF API call failed with error: " //                     &
1545                 TRIM( nf90_strerror(status) )
1546       CALL inifor_abort('io.check', message)
1547    ENDIF
1548
1549 END SUBROUTINE check
1550
1551
1552!------------------------------------------------------------------------------!
1553! Description:
1554! ------------
1555!> Setup the origin of the PALM coordinate system based on what is given in the
1556!> INIFOR namelist file and via an optional static driver.
1557!------------------------------------------------------------------------------!
1558 SUBROUTINE set_palm_origin( cfg, namelist_longitude, namelist_latitude,       &
1559                             origin_lon, origin_lat, z0 )
1560
1561    TYPE(inifor_config), INTENT(IN) ::  cfg
1562    REAL(wp), INTENT(IN)            ::  namelist_longitude, namelist_latitude
1563    REAL(wp), INTENT(OUT)           ::  origin_lon, origin_lat, z0
1564
1565    message = 'Reading PALM-4U origin from'
1566    IF ( TRIM( cfg%static_driver_file ) .NE. '' )  THEN
1567
1568       origin_lon = get_netcdf_attribute( cfg%static_driver_file, TRIM( PIDS_ORIGIN_LON ) )
1569       origin_lat = get_netcdf_attribute( cfg%static_driver_file, TRIM( PIDS_ORIGIN_LAT ) )
1570       z0         = get_netcdf_attribute( cfg%static_driver_file, TRIM( PIDS_ORIGIN_Z ) )
1571
1572       message = TRIM(message) // " static driver file '"                      &
1573                               // TRIM( cfg%static_driver_file ) // "'"
1574
1575
1576    ELSE
1577
1578       origin_lon = namelist_longitude
1579       origin_lat = namelist_latitude
1580
1581       message = TRIM( message ) // " namlist file '"                          &
1582                                 // TRIM( cfg%namelist_file ) // "'"
1583
1584    ENDIF
1585    origin_lon = origin_lon * TO_RADIANS
1586    origin_lat = origin_lat * TO_RADIANS
1587
1588    CALL report('set_palm_origin', message)
1589
1590    IF ( cfg%z0_is_set )  THEN
1591       z0 = cfg%z0
1592       IF ( TRIM( cfg%static_driver_file ) .NE. '' )  THEN
1593          message = 'You specified both --static-driver/-t and --elevation/-z0. ' // &
1594                    'Using the command line value (' // TRIM( real_to_str_f( cfg%z0 ) ) // &
1595                    ') instead of static driver value (' // TRIM( real_to_str_f( z0 ) ) // ').'
1596          CALL warn( 'set_palm_origin', message )
1597       ENDIF
1598    ENDIF
1599
1600 END SUBROUTINE set_palm_origin
1601
1602
1603!------------------------------------------------------------------------------!
1604! Description:
1605! ------------
1606! This function is meant to check weather a COSMO soil variable has an
1607! additional and redunant surface value at depth = 0.0. For instance operational
1608! DWD COSMO output contains the surface temperature in T_SO as a copy of the
1609! values in the first soil layer.
1610!------------------------------------------------------------------------------!
1611 LOGICAL FUNCTION has_surface_value( soil_var, filename ) 
1612
1613    TYPE(nc_var), INTENT(IN)     ::  soil_var
1614    CHARACTER(LEN=*), INTENT(IN) ::  filename
1615   
1616    REAL(wp), ALLOCATABLE        ::  depths(:)
1617
1618    CALL get_dimension_vector_of_variable(                                     &
1619       soil_var%name,                                                          &
1620       dim_idx = NC_DEPTH_DIM_IDX,                                             &
1621       filename = filename,                                                    &
1622       dim_vector = depths                                                     &
1623    )
1624
1625    has_surface_value = nearly_equal( depths(1), 0.0_wp, 10 * EPSILON(1.0_wp) )
1626
1627 END FUNCTION has_surface_value 
1628
1629
1630!------------------------------------------------------------------------------!
1631! Description:
1632! ------------
1633! This routine reads the dim_idx-th dimension vector of the variable varname
1634! from netCDF file filename. It is used for finding the depth coordinate vector
1635! of COSMO soil variables without knowing its name.
1636!------------------------------------------------------------------------------!
1637 SUBROUTINE get_dimension_vector_of_variable( varname, dim_idx, filename, dim_vector )
1638    CHARACTER(LEN=*), INTENT(IN)              ::  varname, filename
1639    INTEGER, INTENT(IN)                       ::  dim_idx
1640
1641    REAL(wp), INTENT(OUT), ALLOCATABLE ::  dim_vector(:)
1642
1643    INTEGER                            ::  dimids(NF90_MAX_VAR_DIMS)
1644    INTEGER                            ::  varid
1645    CHARACTER(LEN=NF90_MAX_NAME)       ::  dimname
1646
1647    INTEGER ::  ncid
1648
1649    IF ( nf90_open( TRIM( filename ), NF90_NOWRITE, ncid ) .EQ. NF90_NOERR )  THEN
1650
1651!
1652!--    get id of variable varname
1653       CALL check( nf90_inq_varid( ncid, TRIM( varname ), varid ) )
1654       
1655!
1656!--    get dimension ids of variable with varid
1657       CALL check( nf90_inquire_variable( ncid, varid, dimids = dimids ) )
1658
1659!
1660!--    get name of dim_idx-th dimension variable
1661       CALL check( nf90_inquire_dimension( ncid, dimids(dim_idx), name = dimname ) )
1662       CALL check( nf90_close( ncid ) )
1663
1664    ELSE
1665
1666       message = "Failed to open file '" // TRIM(filename) // "'."
1667       CALL inifor_abort('get_netcdf_variable', message)
1668
1669    ENDIF
1670
1671    ! get dimension vector with dimname
1672    CALL get_netcdf_dim_vector( filename, dimname, dim_vector )
1673
1674 END SUBROUTINE get_dimension_vector_of_variable
1675
1676
1677 LOGICAL FUNCTION netcdf_variable_present_in_file( varname, filename )
1678    CHARACTER(LEN=SNAME), INTENT(IN) ::  varname
1679    CHARACTER(LEN=PATH), INTENT(IN)  ::  filename
1680
1681    INTEGER ::  ncid, varid
1682   
1683    netcdf_variable_present_in_file = (                                        &
1684       nf90_open( TRIM( filename ), NF90_NOWRITE, ncid ) .EQ. NF90_NOERR .AND. &
1685       nf90_inq_varid( ncid, varname, varid ) .EQ. NF90_NOERR .AND.            &
1686       nf90_close( ncid ) .EQ. NF90_NOERR                                      &
1687    )
1688 
1689 END FUNCTION netcdf_variable_present_in_file
1690
1691
1692 END MODULE inifor_io
1693#endif
Note: See TracBrowser for help on using the repository browser.