!> @file src/inifor_io.f90
!------------------------------------------------------------------------------!
! This file is part of the PALM model system.
!
! PALM is free software: you can redistribute it and/or modify it under the
! terms of the GNU General Public License as published by the Free Software
! Foundation, either version 3 of the License, or (at your option) any later
! version.
!
! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
! A PARTICULAR PURPOSE. See the GNU General Public License for more details.
!
! You should have received a copy of the GNU General Public License along with
! PALM. If not, see .
!
! Copyright 2017-2021 Leibniz Universitaet Hannover
! Copyright 2017-2021 Deutscher Wetterdienst Offenbach
!------------------------------------------------------------------------------!
!
! Current revisions:
! -----------------
!
!
! Former revisions:
! -----------------
! $Id: inifor_io.f90 4843 2021-01-15 15:22:11Z raasch $
! Validate mesoscale input data by the number of dimensions and their lengths
!
!
! 4675 2020-09-11 10:00:26Z eckhard
! New command-line options for soil profile initialization and homogeneous
! boundary conditions
!
!
! 4659 2020-08-31 11:21:17Z eckhard
! New command-line option '--precipitation' enables output of precipitation
! surface forcing
! Only define netCDF variables in enabled IO groups
! Improved code formatting
!
!
! 4569 2020-06-19 12:46:17Z eckhard
! Removed unused variable
!
!
! 4568 2020-06-19 11:56:30Z eckhard
! Handle COSMO soil data with and without additional surface temperature
!
!
! 4553 2020-06-03 16:34:15Z eckhard
! Option --help now points the user to INIFOR's wiki page
! Remove deprecated command-line options -clon and -clat
!
!
! 4538 2020-05-18 13:45:35Z eckhard
! Make setting the vertical PALM origin mandatory
!
!
! 4523 2020-05-07 15:58:16Z eckhard
! respect integer working precision (iwp) specified in inifor_defs.f90
!
!
! 4481 2020-03-31 18:55:54Z maronga
! Pass hhl_file directly instead of entire INIFOR configuration
!
!
! 3997 2019-05-23 12:35:57Z eckhard
! Added boolean indicator for --elevation option invocation
! Stop INIFOR if no command-line options given
!
!
! 3866 2019-04-05 14:25:01Z eckhard
! Use PALM's working precision
! Improved coding style
!
!
! 3801 2019-03-15 17:14:25Z eckhard
! Added routine get_cosmo_grid() to read in COSMO rotated pole from COSMO domain
! Moved get_soil_layer_thickness() here from inifor_grid
!
! 3785 2019-03-06 10:41:14Z eckhard
! Temporariliy disabled height-based geostrophic wind averaging
! Improved variable naming
!
!
! 3764 2019-02-26 13:42:09Z eckhard
! Removed dependency on radiation input files
!
!
! 3716 2019-02-05 17:02:38Z eckhard
! Removed dependency on soilmoisture input files
!
!
! 3680 2019-01-18 14:54:12Z knoop
! Moved get_input_file_list() here from grid module, added check for presence of
! input files
!
!
! 3618 2018-12-10 13:25:22Z eckhard
! Prefixed all INIFOR modules with inifor_
!
!
! 3615 2018-12-10 07:21:03Z raasch
! bugfix: abort replaced by inifor_abort
!
! 3557 2018-11-22 16:01:22Z eckhard
! Updated documentation, removed unused subroutine write_netcdf_variable_2d()
!
!
! 3537 2018-11-20 10:53:14Z eckhard
! New routine get_netcdf_dim_vector()
!
!
! 3534 2018-11-19 15:35:16Z raasch
! bugfix: INTENT attribute changed
!
! 3456 2018-10-30 14:29:54Z eckhard
! NetCDf output of internal arrays only with --debug option
!
!
! 3447 2018-10-29 15:52:54Z eckhard
! Removed INCLUDE statement for get_netcdf_variable()
! Renamed source files for compatibilty with PALM build system
!
!
! 3395 2018-10-22 17:32:49Z eckhard
! Added command-line options for configuring the computation of geostrophic
! winds (--averagin-mode, --averaging-angle)
! Added command-line option --input-prefix for setting input file prefixes all
! at once
! Added --debug option for more verbose terminal output
! Added option-specific *_is_set LOGICALs to indicate invocation from the
! command-line
! Improved error messages in case of empty file-name strings
! Improved routine naming
!
! 3183 2018-07-27 14:25:55Z suehring
! Introduced new PALM grid stretching
! Updated variable names and metadata for PIDS v1.9 compatibility
! Improved handling of the start date string
! Better compatibility with older Intel compilers:
! - avoiding implicit array allocation with new get_netcdf_variable()
! subroutine instead of function
! Improved command line interface:
! - Added configuration validation
! - New options to configure input file prefixes
! - GNU-style short and long option names
! - Added version and copyright output
!
!
! 3182 2018-07-27 13:36:03Z suehring
! Initial revision
!
!
!
! Authors:
! --------
!> @author Eckhard Kadasch (Deutscher Wetterdienst, Offenbach)
!
! Description:
! ------------
!> The io module contains the functions needed to read and write netCDF data in
!> INIFOR.
!------------------------------------------------------------------------------!
#if defined ( __netcdf )
MODULE inifor_io
USE inifor_control
USE inifor_defs, &
ONLY: CFG_INIT_PROFILE, CFG_INIT_VOLUME, &
CFG_INIT_SOIL_PROFILE, CFG_INIT_SOIL_VOLUME, &
CFG_FORCING_HETERO, CFG_FORCING_HOMO, CFG_FORCING_NUDGING, &
DATE, SNAME, PATH, PI, TO_RADIANS, TO_DEGREES, VERSION, &
NC_DEPTH_DIM_IDX, NC_DEPTH_NAME, NC_HHL_NAME, NC_RLAT_NAME, &
NC_RLON_NAME, NC_ROTATED_POLE_NAME, NC_POLE_LATITUDE_NAME, &
NC_POLE_LONGITUDE_NAME, RHO_L, iwp, wp, &
PIDS_ORIGIN_LON, PIDS_ORIGIN_LAT, PIDS_ORIGIN_Z
USE inifor_types
USE inifor_util, &
ONLY: add_hours_to, nearly_equal, reverse, str, real_to_str
USE netcdf
IMPLICIT NONE
!------------------------------------------------------------------------------!
! Description:
! ------------
!> get_netcdf_variable() reads the netCDF data and metadate for the netCDF
!> variable 'in_var%name' from the file 'in_file'. The netCDF data array is
!> stored in the 'buffer' array and metadata added to the respective members of
!> the given 'in_var'.
!------------------------------------------------------------------------------!
INTERFACE get_netcdf_variable
MODULE PROCEDURE get_netcdf_variable_int
MODULE PROCEDURE get_netcdf_variable_real
END INTERFACE get_netcdf_variable
PRIVATE :: get_netcdf_variable_int, get_netcdf_variable_real
CONTAINS
!------------------------------------------------------------------------------!
! Description:
! ------------
!> get_netcdf_variable_int() implements the integer variant for the
!> get_netcdf_variable interface.
!------------------------------------------------------------------------------!
SUBROUTINE get_netcdf_variable_int(in_file, in_var, buffer)
CHARACTER(LEN=PATH), INTENT(IN) :: in_file
TYPE(nc_var), INTENT(INOUT) :: in_var
INTEGER(iwp), ALLOCATABLE, INTENT(INOUT) :: buffer(:,:,:)
INTEGER :: ncid
INTEGER, DIMENSION(3) :: start, count
IF ( nf90_open( TRIM(in_file), NF90_NOWRITE, ncid ) .EQ. NF90_NOERR .AND. &
nf90_inq_varid( ncid, in_var%name, in_var%varid ) .EQ. NF90_NOERR ) THEN
CALL get_input_dimensions(in_var, ncid)
CALL get_netcdf_start_and_count(in_var, start, count)
CALL log_runtime('time', 'read')
ALLOCATE( buffer( count(1), count(2), count(3) ) )
CALL log_runtime('time', 'alloc')
CALL check(nf90_get_var( ncid, in_var%varid, buffer, &
start = start, &
count = count ))
ELSE
message = "Failed to read '" // TRIM(in_var%name) // &
"' from file '" // TRIM(in_file) // "'."
CALL inifor_abort('get_netcdf_variable', message)
ENDIF
CALL check(nf90_close(ncid))
CALL log_runtime('time', 'read')
END SUBROUTINE get_netcdf_variable_int
!------------------------------------------------------------------------------!
! Description:
! ------------
!> get_netcdf_variable_real() implements the real variant for the
!> get_netcdf_variable interface.
!------------------------------------------------------------------------------!
SUBROUTINE get_netcdf_variable_real(in_file, in_var, buffer)
CHARACTER(LEN=PATH), INTENT(IN) :: in_file
TYPE(nc_var), INTENT(INOUT) :: in_var
REAL(wp), ALLOCATABLE, INTENT(INOUT) :: buffer(:,:,:)
INTEGER :: ncid
INTEGER, DIMENSION(3) :: start, count
IF ( nf90_open( TRIM(in_file), NF90_NOWRITE, ncid ) .EQ. NF90_NOERR .AND. &
nf90_inq_varid( ncid, in_var%name, in_var%varid ) .EQ. NF90_NOERR ) THEN
CALL get_input_dimensions(in_var, ncid)
CALL get_netcdf_start_and_count(in_var, start, count)
CALL log_runtime('time', 'read')
ALLOCATE( buffer( count(1), count(2), count(3) ) )
CALL log_runtime('time', 'alloc')
CALL check(nf90_get_var( ncid, in_var%varid, buffer, &
start = start, &
count = count ))
ELSE
message = "Failed to read '" // TRIM(in_var%name) // &
"' from file '" // TRIM(in_file) // "'."
CALL inifor_abort('get_netcdf_variable', message)
ENDIF
CALL check(nf90_close(ncid))
CALL log_runtime('time', 'read')
END SUBROUTINE get_netcdf_variable_real
!------------------------------------------------------------------------------!
! Description:
! ------------
!> get_netcdf_dim_vector() reads the coordinate array 'coordname' from the
!> netCDF file 'filename'.
!------------------------------------------------------------------------------!
SUBROUTINE get_netcdf_dim_vector(filename, coordname, coords)
CHARACTER(LEN=*), INTENT(IN) :: filename
CHARACTER(LEN=*), INTENT(IN) :: coordname
REAL(wp), ALLOCATABLE, INTENT(INOUT) :: coords(:)
INTEGER :: ncid, varid, dimlen
INTEGER :: dimids(NF90_MAX_VAR_DIMS)
IF ( nf90_open( TRIM(filename), NF90_NOWRITE, ncid ) .EQ. NF90_NOERR .AND. &
nf90_inq_varid( ncid, coordname, varid ) .EQ. NF90_NOERR ) THEN
CALL check(nf90_inquire_variable( ncid, varid, dimids = dimids ))
CALL check(nf90_inquire_dimension( ncid, dimids(1), len = dimlen ))
ALLOCATE(coords(dimlen))
CALL check(nf90_get_var( ncid, varid, coords))
ELSE
message = "Failed to read '" // TRIM(coordname) // &
"' from file '" // TRIM(filename) // "'."
CALL inifor_abort('get_netcdf_dim_vector', message)
ENDIF
END SUBROUTINE get_netcdf_dim_vector
!------------------------------------------------------------------------------!
! Description:
! ------------
!> get_input_dimensions() reads dimensions metadata of the netCDF variable given
!> by 'in_var%name' into 'in_var' data structure.
!------------------------------------------------------------------------------!
SUBROUTINE get_input_dimensions(in_var, ncid)
TYPE(nc_var), INTENT(INOUT) :: in_var
INTEGER, INTENT(IN) :: ncid
INTEGER :: i
CALL check(nf90_get_att( ncid, in_var%varid, "long_name", &
in_var%long_name))
CALL check(nf90_get_att( ncid, in_var%varid, "units", &
in_var%units ))
CALL check(nf90_inquire_variable( ncid, in_var%varid, &
ndims = in_var%ndim, &
dimids = in_var%dimids ))
DO i = 1, in_var%ndim
CALL check(nf90_inquire_dimension( ncid, in_var%dimids(i), &
name = in_var%dimname(i), &
len = in_var%dimlen(i) ))
ENDDO
END SUBROUTINE get_input_dimensions
!------------------------------------------------------------------------------!
! Description:
! ------------
!> get_netcdf_start_and_count() gets the start position and element counts for
!> the given netCDF file. This information is used in get_netcdf_variable_int()
!> and _real() for reading input variables..
!------------------------------------------------------------------------------!
SUBROUTINE get_netcdf_start_and_count(in_var, start, count)
TYPE(nc_var), INTENT(INOUT) :: in_var
INTEGER, DIMENSION(3), INTENT(OUT) :: start, count
INTEGER :: ndim
IF ( in_var%ndim .LT. 2 .OR. in_var%ndim .GT. 4 ) THEN
message = "Failed reading NetCDF variable " // &
TRIM(in_var%name) // " with " // &
TRIM(str(INT(in_var%ndim, kind=iwp))) // &
" dimensions because only two- and and three-dimensional" // &
" variables are supported."
CALL inifor_abort('get_netcdf_start_and_count', message)
ENDIF
start = (/ 1, 1, 1 /)
IF ( TRIM(in_var%name) .EQ. 'T_SO' .AND. &
in_var%has_redundant_first_level ) THEN
!
!-- Skip depth = 0.0 for T_SO and reduce number of depths from 9 to 8
in_var%dimlen(3) = in_var%dimlen(3) - 1
!
!-- Start reading from second level, e.g. depth = 0.005 instead of 0.0
start(3) = 2
ENDIF
IF (in_var%ndim .EQ. 2) THEN
in_var%dimlen(3) = 1
ENDIF
ndim = MIN(in_var%ndim, 3)
count = (/ 1, 1, 1 /)
count(1:ndim) = in_var%dimlen(1:ndim)
END SUBROUTINE get_netcdf_start_and_count
!------------------------------------------------------------------------------!
! Description:
! ------------
!> Routine for defining netCDF variables in the dynamic driver, INIFOR's netCDF
!> output.
!------------------------------------------------------------------------------!
SUBROUTINE netcdf_define_variable(var, ncid)
TYPE(nc_var), INTENT(INOUT) :: var
INTEGER, INTENT(IN) :: ncid
CALL check(nf90_def_var(ncid, var%name, NF90_FLOAT, var%dimids(1:var%ndim), var%varid))
CALL check(nf90_put_att(ncid, var%varid, "long_name", var%long_name))
CALL check(nf90_put_att(ncid, var%varid, "units", var%units))
IF ( var%lod .GE. 0 ) THEN
CALL check(nf90_put_att(ncid, var%varid, "lod", var%lod))
ENDIF
CALL check(nf90_put_att(ncid, var%varid, "source", var%source))
CALL check(nf90_put_att(ncid, var%varid, "_FillValue", NF90_FILL_REAL))
END SUBROUTINE netcdf_define_variable
!------------------------------------------------------------------------------!
! Description:
! ------------
!> netcdf_get_dimensions() reads in all dimensions and their lengths and stores
!> them in the given the 'var' data structure. This information is used later
!> for writing output variables in update_output().
!------------------------------------------------------------------------------!
SUBROUTINE netcdf_get_dimensions(var, ncid)
TYPE(nc_var), INTENT(INOUT) :: var
INTEGER, INTENT(IN) :: ncid
INTEGER :: i
CHARACTER(SNAME) :: null
DO i = 1, var%ndim
CALL check(nf90_inquire_dimension(ncid, var%dimids(i), &
name = null, &
len = var%dimlen(i) ) )
ENDDO
END SUBROUTINE netcdf_get_dimensions
!------------------------------------------------------------------------------!
! Description:
! ------------
!> This routine parses and interpretes the command-line options and stores the
!> resulting settings in the 'cfg' data structure.
!------------------------------------------------------------------------------!
SUBROUTINE parse_command_line_arguments( cfg )
TYPE(inifor_config), INTENT(INOUT) :: cfg
CHARACTER(LEN=PATH) :: option, arg
INTEGER :: arg_count, i
cfg%flow_prefix_is_set = .FALSE.
cfg%input_prefix_is_set = .FALSE.
cfg%p0_is_set = .FALSE.
cfg%radiation_prefix_is_set = .FALSE.
cfg%soil_prefix_is_set = .FALSE.
cfg%precipitation_prefix_is_set = .FALSE.
cfg%process_precipitation = .FALSE.
cfg%static_driver_is_set = .FALSE.
cfg%ug_defined_by_user = .FALSE.
cfg%vg_defined_by_user = .FALSE.
cfg%z0_is_set = .FALSE.
arg_count = COMMAND_ARGUMENT_COUNT()
IF (arg_count .GT. 0) THEN
i = 1
DO WHILE (i .LE. arg_count)
CALL GET_COMMAND_ARGUMENT( i, option )
SELECT CASE( TRIM(option) )
CASE( '--averaging-mode' )
CALL get_option_argument( i, arg )
cfg%averaging_mode = TRIM(arg)
CASE( '-date', '-d', '--date' )
CALL get_option_argument( i, arg )
cfg%start_date = TRIM(arg)
CASE( '--debug' )
cfg%debug = .TRUE.
CASE( '-z0', '-z', '--elevation' )
cfg%z0_is_set = .TRUE.
CALL get_option_argument( i, arg )
READ(arg, *) cfg%z0
CASE( '-p0', '-r', '--surface-pressure' )
cfg%p0_is_set = .TRUE.
CALL get_option_argument( i, arg )
READ(arg, *) cfg%p0
CASE( '-ug', '-u', '--geostrophic-u' )
cfg%ug_defined_by_user = .TRUE.
CALL get_option_argument( i, arg )
READ(arg, *) cfg%ug
CASE( '-vg', '-v', '--geostrophic-v' )
cfg%vg_defined_by_user = .TRUE.
CALL get_option_argument( i, arg )
READ(arg, *) cfg%vg
CASE( '-path', '-p', '--path' )
CALL get_option_argument( i, arg )
cfg%input_path = TRIM(arg)
CASE( '-hhl', '-l', '--hhl-file' )
CALL get_option_argument( i, arg )
cfg%hhl_file = TRIM(arg)
CASE( '--input-prefix')
CALL get_option_argument( i, arg )
cfg%input_prefix = TRIM(arg)
cfg%input_prefix_is_set = .TRUE.
CASE( '-a', '--averaging-angle' )
CALL get_option_argument( i, arg )
READ(arg, *) cfg%averaging_angle
CASE( '-static', '-t', '--static-driver' )
cfg%static_driver_is_set = .TRUE.
CALL get_option_argument( i, arg )
cfg%static_driver_file = TRIM(arg)
CASE( '-soil', '-s', '--soil-file')
CALL get_option_argument( i, arg )
cfg%soiltyp_file = TRIM(arg)
CASE( '--flow-prefix')
CALL get_option_argument( i, arg )
cfg%flow_prefix = TRIM(arg)
cfg%flow_prefix_is_set = .TRUE.
CASE( '--radiation-prefix')
CALL get_option_argument( i, arg )
cfg%radiation_prefix = TRIM(arg)
cfg%radiation_prefix_is_set = .TRUE.
CASE( '--soil-prefix')
CALL get_option_argument( i, arg )
cfg%soil_prefix = TRIM(arg)
cfg%soil_prefix_is_set = .TRUE.
CASE( '--precipitation-prefix')
CALL get_option_argument( i, arg )
cfg%precipitation_prefix = TRIM(arg)
cfg%precipitation_prefix_is_set = .TRUE.
CASE( '--precipitation')
cfg%process_precipitation = .TRUE.
CASE( '-o', '--output' )
CALL get_option_argument( i, arg )
cfg%output_file = TRIM(arg)
CASE( '-n', '--namelist' )
CALL get_option_argument( i, arg )
cfg%namelist_file = TRIM(arg)
CASE( '-i', '--init-mode' )
CALL get_option_argument( i, arg )
cfg%ic_mode = TRIM(arg)
CASE( '--soil-init-mode' )
CALL get_option_argument( i, arg )
cfg%isc_mode = TRIM(arg)
CASE( '-f', '--forcing-mode' )
CALL get_option_argument( i, arg )
cfg%bc_mode = TRIM(arg)
CASE( '--version' )
CALL print_version
STOP
CASE( '--help' )
CALL print_version
PRINT *, ""
PRINT *, &
"For documentation and a list of available command-line options " // NEW_LINE( " " ) // &
" please visit https://palm.muk.uni-hannover.de/trac/wiki/doc/app/iofiles/inifor."
STOP
CASE DEFAULT
message = "unknown option '" // TRIM(option) // "'."
CALL inifor_abort('parse_command_line_arguments', message)
END SELECT
i = i + 1
ENDDO
ELSE
CALL print_version
CALL report( 'parse_command_line_arguments', 'No arguments present, exiting.' )
STOP
ENDIF
END SUBROUTINE parse_command_line_arguments
SUBROUTINE get_datetime_file_list( start_date_string, start_hour, end_hour, &
step_hour, input_path, prefix, suffix, &
file_list )
CHARACTER (LEN=DATE), INTENT(IN) :: start_date_string
CHARACTER (LEN=*), INTENT(IN) :: prefix, suffix, input_path
INTEGER(iwp), INTENT(IN) :: start_hour, end_hour, step_hour
CHARACTER(LEN=*), ALLOCATABLE, INTENT(INOUT) :: file_list(:)
INTEGER(iwp) :: number_of_intervals, hour, i
CHARACTER(LEN=DATE) :: date_string
number_of_intervals = CEILING( REAL(end_hour - start_hour) / step_hour )
ALLOCATE( file_list(number_of_intervals + 1) )
DO i = 0, number_of_intervals
hour = start_hour + i * step_hour
date_string = add_hours_to(start_date_string, hour)
file_list(i+1) = TRIM(input_path) // TRIM(prefix) // &
TRIM(date_string) // TRIM(suffix) // '.nc'
ENDDO
END SUBROUTINE get_datetime_file_list
!------------------------------------------------------------------------------!
! Description:
! ------------
!> Establish a list of files based on the given start and end times and file
!> prefixes and suffixes.
!------------------------------------------------------------------------------!
SUBROUTINE get_input_file_list( start_date_string, start_hour, end_hour, &
step_hour, input_path, prefix, suffix, &
file_list )
CHARACTER (LEN=DATE), INTENT(IN) :: start_date_string
CHARACTER (LEN=*), INTENT(IN) :: prefix, suffix, input_path
INTEGER(iwp), INTENT(IN) :: start_hour, end_hour, step_hour
CHARACTER(LEN=*), ALLOCATABLE, INTENT(INOUT) :: file_list(:)
CALL get_datetime_file_list( start_date_string, start_hour, end_hour, &
step_hour, input_path, prefix, suffix, &
file_list )
END SUBROUTINE get_input_file_list
!------------------------------------------------------------------------------!
! Description:
! ------------
!> Abort INIFOR if the given file is not present.
!------------------------------------------------------------------------------!
LOGICAL FUNCTION file_is_present( filename, file_kind, message )
CHARACTER(LEN=*), INTENT(IN) :: filename, file_kind
CHARACTER(LEN=*), INTENT(OUT) :: message
file_is_present = file_present(filename)
IF (.NOT. file_is_present) THEN
IF (LEN( TRIM( filename ) ) == 0) THEN
message = "No name was given for the " // TRIM( file_kind ) // " file."
ELSE
message = "The " // TRIM( file_kind ) // " file '" // &
TRIM( filename ) // "' was not found."
ENDIF
ENDIF
END FUNCTION file_is_present
!------------------------------------------------------------------------------!
! Description:
! ------------
!> Abort INIFOR if the given file is not present.
!------------------------------------------------------------------------------!
SUBROUTINE verify_file( file_name, file_kind )
CHARACTER(LEN=*), INTENT(IN) :: file_name, file_kind
IF (.NOT. file_is_present( file_name, file_kind, message )) THEN
CALL inifor_abort( 'verify_file', message )
ENDIF
message = "Set up input file name '" // TRIM( file_name ) // "'"
CALL report( 'verify_file', message )
END SUBROUTINE verify_file
!------------------------------------------------------------------------------!
! Description:
! ------------
!> Get the argument of the i'th command line option, which is at the location
!> i+1 of the argument list.
!------------------------------------------------------------------------------!
SUBROUTINE get_option_argument( i, arg )
CHARACTER(LEN=PATH), INTENT(INOUT) :: arg
INTEGER, INTENT(INOUT) :: i
i = i + 1
CALL GET_COMMAND_ARGUMENT( i, arg )
END SUBROUTINE
!------------------------------------------------------------------------------!
! Description:
! ------------
!> Checks the INIFOR configuration 'cfg' for plausibility.
!------------------------------------------------------------------------------!
SUBROUTINE validate_config( cfg )
TYPE(inifor_config), INTENT(IN) :: cfg
CALL verify_file( cfg%hhl_file, 'HHL' )
CALL verify_file( cfg%namelist_file, 'NAMELIST' )
CALL verify_file( cfg%soiltyp_file, 'SOILTYP' )
!
!-- Only check optional static driver file name, if it has been given.
IF (TRIM( cfg%static_driver_file ) .NE. '') THEN
CALL verify_file( cfg%static_driver_file, 'static driver' )
ENDIF
SELECT CASE( TRIM( cfg%ic_mode ) )
CASE( CFG_INIT_PROFILE, CFG_INIT_VOLUME )
CASE DEFAULT
message = "Initialization mode '" // TRIM( cfg%ic_mode ) //&
"' is not supported. " //&
"Please select either '" // CFG_INIT_PROFILE //"' or '" // &
CFG_INIT_VOLUME //"', " //&
"or omit the -i/--init-mode option entirely, which corresponds "//&
"to the latter."
CALL inifor_abort( 'validate_config', message )
END SELECT
SELECT CASE( TRIM( cfg%isc_mode ) )
CASE( CFG_INIT_SOIL_PROFILE, CFG_INIT_SOIL_VOLUME )
CASE DEFAULT
message = "Soil initialization mode '" // TRIM( cfg%isc_mode ) //&
"' is not supported. " //&
"Please select either '" // CFG_INIT_SOIL_PROFILE //"' or '" // &
CFG_INIT_SOIL_VOLUME //"', " //&
"or omit the --soil-init-mode option entirely, which corresponds "//&
"to the latter."
CALL inifor_abort( 'validate_config', message )
END SELECT
SELECT CASE( TRIM(cfg%bc_mode) )
CASE( CFG_FORCING_HOMO, CFG_FORCING_HETERO, CFG_FORCING_NUDGING )
CASE DEFAULT
message = "Forcing mode '" // TRIM( cfg%bc_mode ) //&
"' is not supported. " //&
"Please select either '" // CFG_FORCING_NUDGING // &
"', '" // CFG_FORCING_HOMO // "', or '" // &
CFG_FORCING_HETERO // "' " //&
"or omit the -f/--forcing-mode option entirely, which corresponds "//&
"to the latter."
CALL inifor_abort( 'validate_config', message )
END SELECT
SELECT CASE( TRIM( cfg%averaging_mode ) )
CASE( 'level' )
CASE( 'height' )
message = "Averaging mode '" // TRIM( cfg%averaging_mode ) //&
"' is currently not supported. " //&
"Please use level-based averaging by selecting 'level', " //&
"or by omitting the --averaging-mode option entirely."
CALL inifor_abort( 'validate_config', message )
CASE DEFAULT
message = "Averaging mode '" // TRIM( cfg%averaging_mode ) //&
"' is not supported. " //&
! "Please select either 'height' or 'level', " //&
! "or omit the --averaging-mode option entirely, which corresponds "//&
! "to the latter."
"Please use level-based averaging by selecting 'level', " //&
"or by omitting the --averaging-mode option entirely."
CALL inifor_abort( 'validate_config', message )
END SELECT
IF ( cfg%ug_defined_by_user .NEQV. cfg%vg_defined_by_user ) THEN
message = "You specified only one component of the geostrophic " // &
"wind. Please specify either both or none."
CALL inifor_abort( 'validate_config', message )
ENDIF
IF ( .NOT. cfg%static_driver_is_set .AND. .NOT. cfg%z0_is_set ) THEN
message = &
"The vertical origin of the PALM grid has not been defined. " // NEW_LINE( " " ) // &
"Please specify the right value for your setup by either " // NEW_LINE( " " ) // &
" - using the command-line option --elevation , or by" // NEW_LINE( " " ) // &
" - specifying a static driver file using --static in order to use " // NEW_LINE( " " ) // &
" use the value of origin_z (and origin_lon and origin_lat) specifed therein."
CALL inifor_abort( 'validate_config', message )
ENDIF
END SUBROUTINE validate_config
!------------------------------------------------------------------------------!
! Description:
! ------------
!> Checks wheather the COSMO grid matches the shape of the meteorological input
!> data by comparing the number of netCDF dimensions and their lengths in the
!> hhl.nc and the first of the *-flow files.
!------------------------------------------------------------------------------!
SUBROUTINE validate_dataset(flow_files, hhl_file)
CHARACTER(LEN=PATH), INTENT(IN) :: flow_files(:) !< paths to files containing atmospheric variables
CHARACTER(LEN=PATH), INTENT(IN) :: hhl_file !< path to file containing the HHL variable (height of half layers)
CHARACTER(SNAME), PARAMETER :: NC_W_NAME = 'W'
TYPE(nc_var) :: hhl_var, flow_var
INTEGER :: dim_idx, ncid, ndims_flow, ndims_hhl, varid
REAL(wp), ALLOCATABLE :: hhl_dim_vector(:), flow_dim_vector(:)
LOGICAL :: dims_have_same_length
IF ( nf90_open( TRIM( flow_files(1) ), NF90_NOWRITE, ncid ) .EQ. NF90_NOERR ) THEN
CALL check( nf90_inq_varid( ncid, NC_W_NAME, varid=varid ) )
CALL check( nf90_inquire_variable( ncid, varid, ndims=ndims_flow ) )
CALL netcdf_get_dimensions( flow_var, ncid )
CALL check( nf90_close( ncid ) )
ELSE
message = "Failed to read netCDF dimensions'" // &
"' from file '" // TRIM( flow_files(1) ) // "'."
CALL inifor_abort( 'validate_dataset', message )
ENDIF
IF ( nf90_open( TRIM( hhl_file ), NF90_NOWRITE, ncid ) .EQ. NF90_NOERR ) THEN
CALL check( nf90_inq_varid( ncid, NC_HHL_NAME, varid=varid ) )
CALL check( nf90_inquire_variable( ncid, varid, ndims=ndims_hhl ) )
CALL netcdf_get_dimensions( hhl_var, ncid )
CALL check( nf90_close( ncid ) )
ELSE
message = "Failed to read netCDF dimensions'" // &
"' from file '" // TRIM(hhl_file) // "'."
CALL inifor_abort( 'validate_dataset', message )
ENDIF
!
!-- Compare number dimensions of 'HHL' in hhl file and 'W' in first flow file
IF ( .NOT. ndims_flow .EQ. ndims_hhl ) THEN
message = "Mesoscale data inconsistent. Number of dimensions in the " //&
"hhl file does not match with the meteorologial fields " // &
"in the *-flow files (" // &
"HHL: ndims = " // TRIM( str( ndims_hhl ) ) // ", " // &
"W: ndims = " // TRIM( str( ndims_flow ) ) // ")."
CALL inifor_abort( 'validate_dataset', message )
ENDIF
!
!-- Compare lengths of each dimension, ignoring time (dim_idx = 1)
DO dim_idx = 2, ndims_hhl
CALL get_dimension_vector_of_variable( &
NC_HHL_NAME, &
dim_idx = dim_idx, &
filename = hhl_file, &
dim_vector = hhl_dim_vector &
)
CALL get_dimension_vector_of_variable( &
NC_W_NAME, &
dim_idx = dim_idx, &
filename = flow_files(1), &
dim_vector = flow_dim_vector &
)
dims_have_same_length = SIZE( flow_dim_vector ) .EQ. SIZE( hhl_dim_vector )
IF ( .NOT. dims_have_same_length ) THEN
message = &
"Mesoscale data inconsistent. Number of grid points " // &
"in dimension #" // TRIM( str( dim_idx ) ) // &
" do not match in the hhl and *-flow files (" // &
"HHL: len = " // TRIM( str( SIZE( hhl_dim_vector ) ) ) // ", " // &
"W: len = " // TRIM( str( SIZE( flow_dim_vector ) ) )// ")."
CALL inifor_abort( 'validate_dataset', message )
ENDIF
ENDDO
END SUBROUTINE validate_dataset
SUBROUTINE get_cosmo_grid( hhl_file, soil_file, rlon, rlat, hhl, hfl, depths, &
d_depth, d_depth_rho_inv, phi_n, lambda_n, &
phi_equat, &
lonmin_cosmo, lonmax_cosmo, &
latmin_cosmo, latmax_cosmo, &
nlon, nlat, nlev, ndepths )
CHARACTER(LEN=PATH), INTENT(IN) :: hhl_file !< path to file containing the HHL variable (height of half layers)
CHARACTER(LEN=PATH), INTENT(IN) :: soil_file !< path to one of the soil input files for reading soil layer depths
REAL(wp), DIMENSION(:), ALLOCATABLE, INTENT(OUT) :: rlon !< longitudes of COSMO-DE's rotated-pole grid
REAL(wp), DIMENSION(:), ALLOCATABLE, INTENT(OUT) :: rlat !< latitudes of COSMO-DE's rotated-pole grid
REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, INTENT(OUT) :: hhl !< heights of half layers (cell faces) above sea level in COSMO-DE, read in from external file
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
REAL(wp), DIMENSION(:), ALLOCATABLE, INTENT(OUT) :: depths !< COSMO-DE's TERRA-ML soil layer depths
REAL(wp), DIMENSION(:), ALLOCATABLE, INTENT(OUT) :: d_depth
REAL(wp), DIMENSION(:), ALLOCATABLE, INTENT(OUT) :: d_depth_rho_inv
REAL(wp), INTENT(OUT) :: phi_n
REAL(wp), INTENT(OUT) :: phi_equat
REAL(wp), INTENT(OUT) :: lambda_n
REAL(wp), INTENT(OUT) :: lonmin_cosmo !< Minimunm longitude of COSMO-DE's rotated-pole grid [COSMO rotated-pole rad]
REAL(wp), INTENT(OUT) :: lonmax_cosmo !< Maximum longitude of COSMO-DE's rotated-pole grid [COSMO rotated-pole rad]
REAL(wp), INTENT(OUT) :: latmin_cosmo !< Minimunm latitude of COSMO-DE's rotated-pole grid [COSMO rotated-pole rad]
REAL(wp), INTENT(OUT) :: latmax_cosmo !< Maximum latitude of COSMO-DE's rotated-pole grid [COSMO rotated-pole rad]
INTEGER(iwp), INTENT(OUT) :: nlon, nlat, nlev, ndepths
TYPE(nc_var) :: cosmo_var !< COSMO dummy variable, used for reading HHL, rlon, rlat
INTEGER(iwp) :: k
!
!-- Read in COSMO's heights of half layers (vertical cell faces)
cosmo_var%name = NC_HHL_NAME
CALL get_netcdf_variable( hhl_file, cosmo_var, hhl )
CALL get_netcdf_dim_vector( hhl_file, NC_RLON_NAME, rlon )
CALL get_netcdf_dim_vector( hhl_file, NC_RLAT_NAME, rlat )
CALL get_netcdf_dim_vector( soil_file, NC_DEPTH_NAME, depths)
CALL log_runtime( 'time', 'read' )
CALL reverse( hhl )
nlon = SIZE( hhl, 1 )
nlat = SIZE( hhl, 2 )
nlev = SIZE( hhl, 3 )
ndepths = SIZE( depths )
CALL log_runtime( 'time', 'comp' )
ALLOCATE( hfl( nlon, nlat, nlev-1 ) )
ALLOCATE( d_depth( ndepths ), d_depth_rho_inv( ndepths ) )
CALL log_runtime('time', 'alloc')
CALL get_soil_layer_thickness( depths, d_depth )
d_depth_rho_inv = 1.0_wp / ( d_depth * RHO_L )
!
!-- Compute COSMO's heights of full layers (cell centres)
DO k = 1, nlev-1
hfl(:,:,k) = 0.5_wp * ( hhl(:,:,k) + &
hhl(:,:,k+1) )
ENDDO
!
!-- COSMO rotated pole coordinates
phi_n = TO_RADIANS &
* get_netcdf_variable_attribute( hhl_file, &
NC_ROTATED_POLE_NAME, &
NC_POLE_LATITUDE_NAME )
lambda_n = TO_RADIANS &
* get_netcdf_variable_attribute( hhl_file, &
NC_ROTATED_POLE_NAME, &
NC_POLE_LONGITUDE_NAME )
phi_equat = 90.0_wp * TO_RADIANS - phi_n
lonmin_cosmo = MINVAL( rlon ) * TO_RADIANS
lonmax_cosmo = MAXVAL( rlon ) * TO_RADIANS
latmin_cosmo = MINVAL( rlat ) * TO_RADIANS
latmax_cosmo = MAXVAL( rlat ) * TO_RADIANS
CALL log_runtime('time', 'comp')
END SUBROUTINE get_cosmo_grid
!------------------------------------------------------------------------------!
! Description:
! ------------
!> Fills the thickness array of the COSMO soil layers. Since COSMO's (i.e.
!> TERRA_ML's [1]) soil layer boundaries follow the rule
!>
!> depth(0) = 0.0, and
!> depth(k) = 0.01 * 3**(k-1), k in [1,2,3,...,7]
!>
!> and full levels are defined as the midpoints between two layer boundaries,
!> all except the first layer thicknesses equal the depth of the midpoint.
!>
!> [1] A Description of the Nonhydrostatic Regional COSMO Model Part II :
!> Physical Parameterization*, Sect. 11 TERRA_ML.
!> http://www.cosmo-model.org/content/model/documentation/core/cosmoPhysParamtr.pdf)
!>
!> Input parameters:
!> -----------------
!>
!> depths: array of full soil layer depths (cell midpoints)
!>
!>
!> Output parameters:
!> ------------------
!>
!> d_depth: array of soil layer thicknesses
!>
!------------------------------------------------------------------------------!
SUBROUTINE get_soil_layer_thickness( depths, d_depth )
REAL(wp), INTENT(IN) :: depths(:)
REAL(wp), INTENT(OUT) :: d_depth(:)
d_depth(:) = depths(:)
d_depth(1) = 2.0_wp * depths(1)
END SUBROUTINE get_soil_layer_thickness
!------------------------------------------------------------------------------!
! Description:
! ------------
!> Check whether the given file is present on the filesystem.
!------------------------------------------------------------------------------!
LOGICAL FUNCTION file_present( filename )
CHARACTER(LEN=PATH), INTENT(IN) :: filename
INQUIRE( FILE=filename, EXIST=file_present )
END FUNCTION file_present
!------------------------------------------------------------------------------!
! Description:
! ------------
!> This routine initializes the dynamic driver file, i.e. INIFOR's netCDF output
!> file.
!>
!> Besides writing metadata, such as global attributes, coordinates, variables,
!> in the netCDF file, various netCDF IDs are saved for later, when INIFOR
!> writes the actual data.
!------------------------------------------------------------------------------!
SUBROUTINE setup_netcdf_dimensions( output_file, palm_grid, &
start_date_string, origin_lon, origin_lat )
TYPE(nc_file), INTENT(INOUT) :: output_file
TYPE(grid_definition), INTENT(IN) :: palm_grid
CHARACTER (LEN=DATE), INTENT(IN) :: start_date_string
REAL(wp), INTENT(IN) :: origin_lon, origin_lat
CHARACTER (LEN=8) :: date_string
CHARACTER (LEN=10) :: time_string
CHARACTER (LEN=5) :: zone_string
CHARACTER (LEN=SNAME) :: history_string
INTEGER :: nx, ny, nz, nt
INTEGER :: ncid, dimids(3), dimvarids(3)
REAL(wp) :: z0
message = "Initializing PALM-4U dynamic driver file '" // &
TRIM(output_file%name) // "' and setting up dimensions."
CALL report('setup_netcdf_dimensions', message)
!
!-- Create the netCDF file as in netCDF-4/HDF5 format if __netcdf4 preprocessor flag is given
#if defined( __netcdf4 )
CALL check(nf90_create(TRIM(output_file%name), OR(NF90_CLOBBER, NF90_HDF5), ncid))
#else
CALL check(nf90_create(TRIM(output_file%name), NF90_CLOBBER, ncid))
#endif
!------------------------------------------------------------------------------
!- Section 1: Define NetCDF dimensions and coordinates
!------------------------------------------------------------------------------
nt = SIZE(output_file%time)
nx = palm_grid%nx
ny = palm_grid%ny
nz = palm_grid%nz
z0 = palm_grid%z0
!
!------------------------------------------------------------------------------
!- Section 2: Write global NetCDF attributes
!------------------------------------------------------------------------------
CALL date_and_time(DATE=date_string, TIME=time_string, ZONE=zone_string)
history_string = &
'Created on '// date_string // &
' at ' // time_string(1:2) // ':' // time_string(3:4) // &
' (UTC' // zone_string // ')'
CALL check(nf90_put_att(ncid, NF90_GLOBAL, 'title', 'PALM input file for scenario ...'))
CALL check(nf90_put_att(ncid, NF90_GLOBAL, 'institution', 'Deutscher Wetterdienst, Offenbach'))
CALL check(nf90_put_att(ncid, NF90_GLOBAL, 'author', 'Eckhard Kadasch, eckhard.kadasch@dwd.de'))
CALL check(nf90_put_att(ncid, NF90_GLOBAL, 'history', TRIM(history_string)))
CALL check(nf90_put_att(ncid, NF90_GLOBAL, 'references', '--'))
CALL check(nf90_put_att(ncid, NF90_GLOBAL, 'comment', '--'))
CALL check(nf90_put_att(ncid, NF90_GLOBAL, 'origin_lat', TRIM(real_to_str(origin_lat*TO_DEGREES, '(F18.13)'))))
CALL check(nf90_put_att(ncid, NF90_GLOBAL, 'origin_lon', TRIM(real_to_str(origin_lon*TO_DEGREES, '(F18.13)'))))
!
!-- FIXME: This is the elevation relative to COSMO-DE/D2 sea level and does
!-- FIXME: not necessarily comply with DHHN2016 (c.f. PALM Input Data
!-- FIXME: Standard v1.9., origin_z)
CALL check(nf90_put_att(ncid, NF90_GLOBAL, 'origin_z', TRIM(real_to_str(z0, '(F18.13)'))))
CALL check(nf90_put_att(ncid, NF90_GLOBAL, 'inifor_version', TRIM(VERSION)))
CALL check(nf90_put_att(ncid, NF90_GLOBAL, 'palm_version', '--'))
!
!
!------------------------------------------------------------------------------
!- Section 2a: Define dimensions for cell centers (scalars in soil and atmosph.)
!------------------------------------------------------------------------------
!
!-- reset dimids first
dimids = (/0, 0, 0/)
CALL check( nf90_def_dim(ncid, "x", nx+1, dimids(1)) )
CALL check( nf90_def_dim(ncid, "y", ny+1, dimids(2)) )
CALL check( nf90_def_dim(ncid, "z", nz, dimids(3)) )
!
!-- save dimids for later
output_file%dimids_scl = dimids
!
!-- reset dimvarids first
dimvarids = (/0, 0, 0/)
CALL check(nf90_def_var(ncid, "x", NF90_FLOAT, dimids(1), dimvarids(1)))
CALL check(nf90_put_att(ncid, dimvarids(1), "standard_name", "x coordinate of cell centers"))
CALL check(nf90_put_att(ncid, dimvarids(1), "units", "m"))
CALL check(nf90_def_var(ncid, "y", NF90_FLOAT, dimids(2), dimvarids(2)))
CALL check(nf90_put_att(ncid, dimvarids(2), "standard_name", "y coordinate of cell centers"))
CALL check(nf90_put_att(ncid, dimvarids(2), "units", "m"))
CALL check(nf90_def_var(ncid, "z", NF90_FLOAT, dimids(3), dimvarids(3)))
CALL check(nf90_put_att(ncid, dimvarids(3), "standard_name", "z coordinate of cell centers"))
CALL check(nf90_put_att(ncid, dimvarids(3), "units", "m"))
!
!-- save dimvarids for later
output_file%dimvarids_scl = dimvarids
!
!-- overwrite third dimid with the one of depth
CALL check(nf90_def_dim(ncid, "zsoil", SIZE(palm_grid%depths), dimids(3)) )
!
!-- save dimids for later
output_file%dimids_soil = dimids
!
!-- overwrite third dimvarid with the one of depth
CALL check(nf90_def_var(ncid, "zsoil", NF90_FLOAT, output_file%dimids_soil(3), dimvarids(3)))
CALL check(nf90_put_att(ncid, dimvarids(3), "standard_name", "depth_below_land"))
CALL check(nf90_put_att(ncid, dimvarids(3), "positive", "down"))
CALL check(nf90_put_att(ncid, dimvarids(3), "units", "m"))
!
!-- save dimvarids for later
output_file%dimvarids_soil = dimvarids
!
!------------------------------------------------------------------------------
!- Section 2b: Define dimensions for cell faces/velocities
!------------------------------------------------------------------------------
!
!-- reset dimids first
dimids = (/0, 0, 0/)
CALL check(nf90_def_dim(ncid, "xu", nx, dimids(1)) )
CALL check(nf90_def_dim(ncid, "yv", ny, dimids(2)) )
CALL check(nf90_def_dim(ncid, "zw", nz-1, dimids(3)) )
!
!-- save dimids for later
output_file%dimids_vel = dimids
!
!-- reset dimvarids first
dimvarids = (/0, 0, 0/)
CALL check(nf90_def_var(ncid, "xu", NF90_FLOAT, dimids(1), dimvarids(1)))
CALL check(nf90_put_att(ncid, dimvarids(1), "standard_name", "x coordinate of cell faces"))
CALL check(nf90_put_att(ncid, dimvarids(1), "units", "m"))
CALL check(nf90_def_var(ncid, "yv", NF90_FLOAT, dimids(2), dimvarids(2)))
CALL check(nf90_put_att(ncid, dimvarids(2), "standard_name", "y coordinate of cell faces"))
CALL check(nf90_put_att(ncid, dimvarids(2), "units", "m"))
CALL check(nf90_def_var(ncid, "zw", NF90_FLOAT, dimids(3), dimvarids(3)))
CALL check(nf90_put_att(ncid, dimvarids(3), "standard_name", "z coordinate of cell faces"))
CALL check(nf90_put_att(ncid, dimvarids(3), "units", "m"))
!
!-- save dimvarids for later
output_file%dimvarids_vel = dimvarids
!
!------------------------------------------------------------------------------
!- Section 2c: Define time dimension
!------------------------------------------------------------------------------
CALL check(nf90_def_dim(ncid, "time", nt, output_file%dimid_time) )
CALL check(nf90_def_var(ncid, "time", NF90_FLOAT, &
output_file%dimid_time, &
output_file%dimvarid_time))
CALL check(nf90_put_att(ncid, output_file%dimvarid_time, "standard_name", "time"))
CALL check(nf90_put_att(ncid, output_file%dimvarid_time, "long_name", "time"))
CALL check(nf90_put_att(ncid, output_file%dimvarid_time, "units", &
"seconds since " // start_date_string // " UTC"))
CALL check(nf90_enddef(ncid))
!
!------------------------------------------------------------------------------
!- Section 3: Write grid coordinates
!------------------------------------------------------------------------------
CALL check(nf90_put_var(ncid, output_file%dimvarids_scl(1), palm_grid%x))
CALL check(nf90_put_var(ncid, output_file%dimvarids_scl(2), palm_grid%y))
CALL check(nf90_put_var(ncid, output_file%dimvarids_scl(3), palm_grid%z))
CALL check(nf90_put_var(ncid, output_file%dimvarids_vel(1), palm_grid%xu))
CALL check(nf90_put_var(ncid, output_file%dimvarids_vel(2), palm_grid%yv))
CALL check(nf90_put_var(ncid, output_file%dimvarids_vel(3), palm_grid%zw))
!
!-- TODO Read in soil depths from input file before this.
CALL check(nf90_put_var(ncid, output_file%dimvarids_soil(3), palm_grid%depths))
!
!-- Write time vector
CALL check(nf90_put_var(ncid, output_file%dimvarid_time, output_file%time))
!
!-- Close the file
CALL check(nf90_close(ncid))
END SUBROUTINE setup_netcdf_dimensions
!------------------------------------------------------------------------------!
! Description:
! ------------
!> Defines the netCDF variables to be written to the dynamic driver file
!------------------------------------------------------------------------------!
SUBROUTINE setup_netcdf_variables(filename, io_group_list)
CHARACTER (LEN=*), INTENT(IN) :: filename
TYPE(io_group), INTENT(IN), TARGET :: io_group_list(:)
TYPE(io_group), POINTER :: group
TYPE(nc_var), POINTER :: var
INTEGER(iwp) :: group_idx, var_idx, n_var
INTEGER :: ncid
LOGICAL :: to_be_written
message = "Defining variables in dynamic driver '" // TRIM(filename) // "'."
CALL report('setup_netcdf_variables', message)
CALL check(nf90_open(TRIM(filename), NF90_WRITE, ncid))
CALL check(nf90_redef(ncid))
n_var = 0
DO group_idx = 1, SIZE( io_group_list )
group => io_group_list(group_idx)
DO var_idx = 1, SIZE( group%out_vars )
to_be_written = .FALSE.
IF (ALLOCATED( group%out_vars )) THEN
var => group%out_vars(var_idx)
n_var = n_var + 1
to_be_written = ( &
group%to_be_processed .AND. var%to_be_processed &
.AND. .NOT. var%is_internal &
)
ENDIF
IF ( to_be_written ) THEN
message = " variable #" // TRIM(str(n_var)) // " '" // TRIM(var%name) // "'."
CALL report('setup_netcdf_variables', message)
CALL netcdf_define_variable(var, ncid)
CALL netcdf_get_dimensions(var, ncid)
ENDIF
ENDDO
ENDDO
CALL check(nf90_enddef(ncid))
CALL check(nf90_close(ncid))
message = "Dynamic driver '" // TRIM(filename) // "' initialized successfully."
CALL report('setup_netcdf_variables', message)
END SUBROUTINE setup_netcdf_variables
!------------------------------------------------------------------------------!
! Description:
! ------------
!> This routine reads and returns all input variables of the given IO group
!> It accomodates the data by allocating a container variable that represents a
!> list of arrays of the same length as the groups variable list. (This list
!> will typically contain one or two items.) After the container, its members
!> are allocated one by one with the appropriate, possibly different,
!> dimensions.
!>
!> The 'group' is an INTENT(INOUT) variable so that 'get_netcdf_variable()' can
!> record netCDF IDs in the 'in_var_list()' member variable.
!------------------------------------------------------------------------------!
SUBROUTINE read_input_variables(group, iter, buffer)
TYPE(io_group), INTENT(INOUT), TARGET :: group
INTEGER(iwp), INTENT(IN) :: iter
TYPE(container), ALLOCATABLE, INTENT(INOUT) :: buffer(:)
INTEGER(iwp) :: hour, buf_id
TYPE(nc_var), POINTER :: input_var
CHARACTER(LEN=PATH), POINTER :: input_file
INTEGER(iwp) :: ivar, nbuffers
message = "Reading data for I/O group '" // TRIM(group%in_var_list(1)%name) // "'."
CALL report('read_input_variables', message)
input_file => group%in_files(iter)
!
!------------------------------------------------------------------------------
!- Section 1: Load input buffers for accumulated variables
!------------------------------------------------------------------------------
!
!-- radiation budgets, precipitation
IF ( group%kind == 'running average' .OR. &
group%kind == 'accumulated' ) THEN
IF ( SIZE( group%in_var_list ) .GT. 1 ) THEN
message = "I/O groups may not contain more than one " // &
"accumulated variable. Group '" // TRIM(group%kind) //&
"' contains " // &
TRIM( str( SIZE( group%in_var_list, kind=iwp ) ) ) // "."
CALL inifor_abort( 'read_input_variables | accumulation', message )
ENDIF
!
!-- use two buffer arrays
nbuffers = 2
IF ( .NOT. ALLOCATED( buffer ) ) ALLOCATE( buffer(nbuffers) )
!
!-- hour of the day
hour = iter - 1
!
!-- chose correct buffer array
buf_id = select_buffer(hour)
CALL log_runtime('time', 'read')
IF ( ALLOCATED(buffer(buf_id)%array) ) DEALLOCATE(buffer(buf_id)%array)
CALL log_runtime('time', 'alloc')
input_var => group%in_var_list(1)
CALL get_netcdf_variable(input_file, input_var, buffer(buf_id)%array)
CALL report('read_input_variables', "Read accumulated " // TRIM(group%in_var_list(1)%name))
IF ( input_var%is_upside_down ) CALL reverse(buffer(buf_id)%array)
CALL log_runtime('time', 'comp')
!------------------------------------------------------------------------------
!- Section 2: Load input buffers for normal I/O groups
!------------------------------------------------------------------------------
ELSE
!
!-- Allocate one input buffer per output quantity. If more quantities
!-- have to be computed than input variables exist in this group,
!-- allocate more buffers. For instance, in the thermodynamics group,
!-- there are three input variabels (PP, T, Qv) and four quantities
!-- necessart (P, Theta, Rho, qv) for the corresponding output fields
!-- (p0, Theta, qv, ug, and vg)
ALLOCATE( buffer(group%n_output_quantities) )
CALL log_runtime('time', 'alloc')
!
!-- Read in all input variables, leave extra buffers-if any-untouched.
DO ivar = 1, group%n_inputs
input_var => group%in_var_list(ivar)
IF ( input_var%to_be_processed ) THEN
! Check wheather P or PP is present in input file
IF (input_var%name == 'P') THEN
input_var%name = TRIM( get_pressure_varname(input_file) )
CALL log_runtime('time', 'read')
ENDIF
CALL get_netcdf_variable(input_file, input_var, buffer(ivar)%array)
IF ( input_var%is_upside_down ) CALL reverse(buffer(ivar)%array)
CALL log_runtime('time', 'comp')
ENDIF
ENDDO
ENDIF
END SUBROUTINE read_input_variables
!------------------------------------------------------------------------------!
! Description:
! ------------
!> Select the appropriate buffer ID for accumulated COSMO input variables
!> depending on the current hour.
!------------------------------------------------------------------------------!
INTEGER(iwp) FUNCTION select_buffer(hour)
INTEGER(iwp), INTENT(IN) :: hour
INTEGER(iwp) :: step
select_buffer = 0_iwp
step = MODULO(hour, 3_iwp) + 1_iwp
SELECT CASE(step)
CASE(1, 3)
select_buffer = 1_iwp
CASE(2)
select_buffer = 2_iwp
CASE DEFAULT
message = "Invalid step '" // TRIM(str(step))
CALL inifor_abort('select_buffer', message)
END SELECT
END FUNCTION select_buffer
!------------------------------------------------------------------------------!
! Description:
! ------------
!> Checks if the input_file contains the absolute pressure, 'P', or the pressure
!> perturbation, 'PP', and returns the appropriate string.
!------------------------------------------------------------------------------!
CHARACTER(LEN=2) FUNCTION get_pressure_varname(input_file) RESULT(var)
CHARACTER(LEN=*) :: input_file
INTEGER :: ncid, varid
CALL check(nf90_open( TRIM(input_file), NF90_NOWRITE, ncid ))
IF ( nf90_inq_varid( ncid, 'P', varid ) .EQ. NF90_NOERR ) THEN
var = 'P'
ELSE IF ( nf90_inq_varid( ncid, 'PP', varid ) .EQ. NF90_NOERR ) THEN
var = 'PP'
CALL report('get_pressure_var', 'Using PP instead of P')
ELSE
message = "Failed to read '" // TRIM(var) // &
"' from file '" // TRIM(input_file) // "'."
CALL inifor_abort('get_pressure_var', message)
ENDIF
CALL check(nf90_close(ncid))
END FUNCTION get_pressure_varname
!------------------------------------------------------------------------------!
! Description:
! ------------
!> Read the given global attribute form the given netCDF file.
!------------------------------------------------------------------------------!
FUNCTION get_netcdf_attribute(filename, attribute) RESULT(attribute_value)
CHARACTER(LEN=*), INTENT(IN) :: filename, attribute
REAL(wp) :: attribute_value
INTEGER :: ncid
IF ( nf90_open( TRIM(filename), NF90_NOWRITE, ncid ) == NF90_NOERR ) THEN
CALL check(nf90_get_att(ncid, NF90_GLOBAL, TRIM(attribute), attribute_value))
CALL check(nf90_close(ncid))
ELSE
message = "Failed to read '" // TRIM(attribute) // &
"' from file '" // TRIM(filename) // "'."
CALL inifor_abort('get_netcdf_attribute', message)
ENDIF
END FUNCTION get_netcdf_attribute
!------------------------------------------------------------------------------!
! Description:
! ------------
!> Read the attribute of the given variable form the given netCDF file.
!------------------------------------------------------------------------------!
FUNCTION get_netcdf_variable_attribute(filename, varname, attribute) &
RESULT(attribute_value)
CHARACTER(LEN=*), INTENT(IN) :: filename, varname, attribute
REAL(wp) :: attribute_value
INTEGER :: ncid, varid
IF ( nf90_open( TRIM(filename), NF90_NOWRITE, ncid ) == NF90_NOERR ) THEN
CALL check( nf90_inq_varid( ncid, TRIM( varname ), varid ) )
CALL check( nf90_get_att( ncid, varid, TRIM( attribute ), &
attribute_value ) )
CALL check( nf90_close( ncid ) )
ELSE
message = "Failed to read '" // TRIM( varname ) // ":" // &
TRIM( attribute ) // "' from file '" // TRIM(filename) // "'."
CALL inifor_abort('get_netcdf_variable_attribute', message)
ENDIF
END FUNCTION get_netcdf_variable_attribute
!------------------------------------------------------------------------------!
! Description:
! ------------
!> Updates the dynamic driver with the interpolated field of the current
!> variable at the current time step.
!------------------------------------------------------------------------------!
SUBROUTINE update_output(var, array, iter, output_file, cfg)
TYPE(nc_var), INTENT(IN) :: var
REAL(wp), INTENT(IN) :: array(:,:,:)
INTEGER(iwp), INTENT(IN) :: iter
TYPE(nc_file), INTENT(IN) :: output_file
TYPE(inifor_config) :: cfg
INTEGER :: ncid, ndim, start(4), count(4)
LOGICAL :: var_is_time_dependent
var_is_time_dependent = ( &
var%dimids( var%ndim ) == output_file%dimid_time &
)
!
!-- Skip time dimension for output
ndim = var%ndim
IF ( var_is_time_dependent ) ndim = var%ndim - 1
start(:) = (/1,1,1,1/)
start(ndim+1) = iter
count(1:ndim) = var%dimlen(1:ndim)
CALL check(nf90_open(output_file%name, NF90_WRITE, ncid))
!
!-- Reduce dimension of output array according to variable kind
SELECT CASE (TRIM(var%kind))
CASE( 'init scalar profile', 'init u profile', 'init v profile', &
'init w profile', 'init soil profile' )
CALL check(nf90_put_var( ncid, var%varid, array(1,1,:) ) )
CASE( 'init soil', 'init scalar', 'init u', 'init v', 'init w' )
CALL check(nf90_put_var( ncid, var%varid, array(:,:,:) ) )
CASE( 'left scalar', 'right scalar', 'left w', 'right w' )
CALL check(nf90_put_var( ncid, var%varid, array(1,:,:), &
start=start(1:ndim+1), &
count=count(1:ndim) ) )
IF (.NOT. SIZE(array, 2) .EQ. var%dimlen(1)) THEN
PRINT *, "inifor: update_output: Dimension ", 1, " of variable ", &
TRIM(var%name), " (", var%dimlen(1), &
") does not match the dimension of the output array (", &
SIZE(array, 2), ")."
STOP
ENDIF
CASE( 'north scalar', 'south scalar', 'north w', 'south w' )
CALL check(nf90_put_var( ncid, var%varid, array(:,1,:), &
start=start(1:ndim+1), &
count=count(1:ndim) ) )
CASE( 'surface forcing', 'top scalar', 'top w' )
CALL check(nf90_put_var( ncid, var%varid, array(:,:,1), &
start=start(1:ndim+1), &
count=count(1:ndim) ) )
CASE( 'left u', 'right u', 'left v', 'right v' )
CALL check(nf90_put_var( ncid, var%varid, array(1,:,:), &
start=start(1:ndim+1), &
count=count(1:ndim) ) )
CASE( 'north u', 'south u', 'north v', 'south v' )
CALL check(nf90_put_var( ncid, var%varid, array(:,1,:), &
start=start(1:ndim+1), &
count=count(1:ndim) ) )
CASE( 'top u', 'top v' )
CALL check(nf90_put_var( ncid, var%varid, array(:,:,1), &
start=start(1:ndim+1), &
count=count(1:ndim) ) )
CASE( 'time series' )
CALL check(nf90_put_var( ncid, var%varid, array(1,1,1), &
start=start(1:ndim+1) ) )
CASE( 'constant scalar profile', 'geostrophic', &
'left scalar profile', 'right scalar profile', &
'north scalar profile', 'south scalar profile', &
'left u profile', 'right u profile', &
'north u profile', 'south u profile', &
'left v profile', 'right v profile', &
'north v profile', 'south v profile', &
'top scalar profile', 'top u profile', 'top v profile', &
'top w profile', &
'left w profile', 'right w profile', &
'north w profile', 'south w profile' )
CALL check(nf90_put_var( ncid, var%varid, array(1,1,:), &
start=start(1:ndim+1), &
count=count(1:ndim) ) )
CASE ( 'internal profile' )
IF ( cfg%debug ) THEN
CALL check(nf90_put_var( ncid, var%varid, array(1,1,:), &
start=start(1:ndim+1), &
count=count(1:ndim) ) )
ENDIF
CASE ( 'large-scale scalar forcing', 'large-scale w forcing' )
message = "Doing nothing in terms of writing large-scale forings."
CALL report('update_output', message)
CASE DEFAULT
message = "Variable kind '" // TRIM(var%kind) // &
"' not recognized."
CALL inifor_abort('update_output', message)
END SELECT
CALL check(nf90_close(ncid))
END SUBROUTINE update_output
!------------------------------------------------------------------------------!
! Description:
! ------------
!> Checks the status of a netCDF API call and aborts if an error occured
!------------------------------------------------------------------------------!
SUBROUTINE check(status)
INTEGER, INTENT(IN) :: status
IF (status /= nf90_noerr) THEN
message = "NetCDF API call failed with error: " // &
TRIM( nf90_strerror(status) )
CALL inifor_abort('io.check', message)
ENDIF
END SUBROUTINE check
!------------------------------------------------------------------------------!
! Description:
! ------------
!> Setup the origin of the PALM coordinate system based on what is given in the
!> INIFOR namelist file and via an optional static driver.
!------------------------------------------------------------------------------!
SUBROUTINE set_palm_origin( cfg, namelist_longitude, namelist_latitude, &
origin_lon, origin_lat, z0 )
TYPE(inifor_config), INTENT(IN) :: cfg
REAL(wp), INTENT(IN) :: namelist_longitude, namelist_latitude
REAL(wp), INTENT(OUT) :: origin_lon, origin_lat, z0
message = 'Reading PALM-4U origin from'
IF ( TRIM( cfg%static_driver_file ) .NE. '' ) THEN
origin_lon = get_netcdf_attribute( cfg%static_driver_file, TRIM( PIDS_ORIGIN_LON ) )
origin_lat = get_netcdf_attribute( cfg%static_driver_file, TRIM( PIDS_ORIGIN_LAT ) )
z0 = get_netcdf_attribute( cfg%static_driver_file, TRIM( PIDS_ORIGIN_Z ) )
message = TRIM(message) // " static driver file '" &
// TRIM( cfg%static_driver_file ) // "'"
ELSE
origin_lon = namelist_longitude
origin_lat = namelist_latitude
message = TRIM( message ) // " namlist file '" &
// TRIM( cfg%namelist_file ) // "'"
ENDIF
origin_lon = origin_lon * TO_RADIANS
origin_lat = origin_lat * TO_RADIANS
CALL report('set_palm_origin', message)
IF ( cfg%z0_is_set ) THEN
z0 = cfg%z0
IF ( TRIM( cfg%static_driver_file ) .NE. '' ) THEN
message = 'You specified both --static-driver/-t and --elevation/-z0. ' // &
'Using the command line value (' // TRIM( real_to_str_f( cfg%z0 ) ) // &
') instead of static driver value (' // TRIM( real_to_str_f( z0 ) ) // ').'
CALL warn( 'set_palm_origin', message )
ENDIF
ENDIF
END SUBROUTINE set_palm_origin
!------------------------------------------------------------------------------!
! Description:
! ------------
! This function is meant to check weather a COSMO soil variable has an
! additional and redunant surface value at depth = 0.0. For instance operational
! DWD COSMO output contains the surface temperature in T_SO as a copy of the
! values in the first soil layer.
!------------------------------------------------------------------------------!
LOGICAL FUNCTION has_surface_value( soil_var, filename )
TYPE(nc_var), INTENT(IN) :: soil_var
CHARACTER(LEN=*), INTENT(IN) :: filename
REAL(wp), ALLOCATABLE :: depths(:)
CALL get_dimension_vector_of_variable( &
soil_var%name, &
dim_idx = NC_DEPTH_DIM_IDX, &
filename = filename, &
dim_vector = depths &
)
has_surface_value = nearly_equal( depths(1), 0.0_wp, 10 * EPSILON(1.0_wp) )
END FUNCTION has_surface_value
!------------------------------------------------------------------------------!
! Description:
! ------------
! This routine reads the dim_idx-th dimension vector of the variable varname
! from netCDF file filename. It is used for finding the depth coordinate vector
! of COSMO soil variables without knowing its name.
!------------------------------------------------------------------------------!
SUBROUTINE get_dimension_vector_of_variable( varname, dim_idx, filename, dim_vector )
CHARACTER(LEN=*), INTENT(IN) :: varname, filename
INTEGER, INTENT(IN) :: dim_idx
REAL(wp), INTENT(OUT), ALLOCATABLE :: dim_vector(:)
INTEGER :: dimids(NF90_MAX_VAR_DIMS)
INTEGER :: varid
CHARACTER(LEN=NF90_MAX_NAME) :: dimname
INTEGER :: ncid
IF ( nf90_open( TRIM( filename ), NF90_NOWRITE, ncid ) .EQ. NF90_NOERR ) THEN
!
!-- get id of variable varname
CALL check( nf90_inq_varid( ncid, TRIM( varname ), varid ) )
!
!-- get dimension ids of variable with varid
CALL check( nf90_inquire_variable( ncid, varid, dimids = dimids ) )
!
!-- get name of dim_idx-th dimension variable
CALL check( nf90_inquire_dimension( ncid, dimids(dim_idx), name = dimname ) )
CALL check( nf90_close( ncid ) )
ELSE
message = "Failed to open file '" // TRIM(filename) // "'."
CALL inifor_abort('get_netcdf_variable', message)
ENDIF
! get dimension vector with dimname
CALL get_netcdf_dim_vector( filename, dimname, dim_vector )
END SUBROUTINE get_dimension_vector_of_variable
LOGICAL FUNCTION netcdf_variable_present_in_file( varname, filename )
CHARACTER(LEN=SNAME), INTENT(IN) :: varname
CHARACTER(LEN=PATH), INTENT(IN) :: filename
INTEGER :: ncid, varid
netcdf_variable_present_in_file = ( &
nf90_open( TRIM( filename ), NF90_NOWRITE, ncid ) .EQ. NF90_NOERR .AND. &
nf90_inq_varid( ncid, varname, varid ) .EQ. NF90_NOERR .AND. &
nf90_close( ncid ) .EQ. NF90_NOERR &
)
END FUNCTION netcdf_variable_present_in_file
END MODULE inifor_io
#endif