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

Last change on this file since 4793 was 4790, checked in by eckhard, 4 years ago

inifor: Validate netCDF dimensions of COSMO input files

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