Home

Context Navigation

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

Last change on this file since 3716 was 3716, checked in by eckhard, 5 years ago
inifor: bugfix: removed dependency on soilmoisture input files; added netcdf preprocessor flag
Property svn:keywords set to `Id`
File size: 49.8 KB

Rev	Line
[3447]	1	!> @file src/inifor_io.f90
[2696]	2	!------------------------------------------------------------------------------!
[2718]	3	! This file is part of the PALM model system.
[2696]	4	!
[2718]	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
[2696]	8	! version.
	9	!
[2718]	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.
[2696]	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	!
[2718]	17	! Copyright 2017-2018 Leibniz Universitaet Hannover
	18	! Copyright 2017-2018 Deutscher Wetterdienst Offenbach
[2696]	19	!------------------------------------------------------------------------------!
	20	!
	21	! Current revisions:
	22	! -----------------
[3183]	23	!
	24	!
	25	! Former revisions:
	26	! -----------------
	27	! $Id: inifor_io.f90 3716 2019-02-05 17:02:38Z eckhard $
[3716]	28	! Removed dependency on soilmoisture input files
	29	!
	30	!
	31	! 3680 2019-01-18 14:54:12Z knoop
[3678]	32	! Moved get_input_file_list() here from grid module, added check for presence of
	33	! input files
	34	!
	35	!
	36	!
	37	!
	38	! 3618 2018-12-10 13:25:22Z eckhard
[3618]	39	! Prefixed all INIFOR modules with inifor_
	40	!
	41	!
	42	! 3615 2018-12-10 07:21:03Z raasch
[3615]	43	! bugfix: abort replaced by inifor_abort
	44	!
	45	! 3557 2018-11-22 16:01:22Z eckhard
[3557]	46	! Updated documentation, removed unused subroutine write_netcdf_variable_2d()
	47	!
	48	!
	49	! 3537 2018-11-20 10:53:14Z eckhard
[3537]	50	! New routine get_netcdf_dim_vector()
	51	!
	52	!
	53	! 3534 2018-11-19 15:35:16Z raasch
[3534]	54	! bugfix: INTENT attribute changed
	55	!
	56	! 3456 2018-10-30 14:29:54Z eckhard
[3456]	57	! NetCDf output of internal arrays only with --debug option
	58	!
	59	!
	60	! 3447 2018-10-29 15:52:54Z eckhard
[3447]	61	! Removed INCLUDE statement for get_netcdf_variable()
	62	! Renamed source files for compatibilty with PALM build system
	63	!
	64	!
	65	! 3395 2018-10-22 17:32:49Z eckhard
[3395]	66	! Added command-line options for configuring the computation of geostrophic
	67	! winds (--averagin-mode, --averaging-angle)
	68	! Added command-line option --input-prefix for setting input file prefixes all
	69	! at once
	70	! Added --debug option for more verbose terminal output
	71	! Added option-specific *_is_set LOGICALs to indicate invocation from the
	72	! command-line
	73	! Improved error messages in case of empty file-name strings
	74	! Improved routine naming
[3262]	75	!
	76	! 3183 2018-07-27 14:25:55Z suehring
[3182]	77	! Introduced new PALM grid stretching
	78	! Updated variable names and metadata for PIDS v1.9 compatibility
	79	! Improved handling of the start date string
	80	! Better compatibility with older Intel compilers:
	81	! - avoiding implicit array allocation with new get_netcdf_variable()
	82	! subroutine instead of function
	83	! Improved command line interface:
	84	! - Added configuration validation
	85	! - New options to configure input file prefixes
	86	! - GNU-style short and long option names
	87	! - Added version and copyright output
[2696]	88	!
[3182]	89	!
[3183]	90	! 3182 2018-07-27 13:36:03Z suehring
[2696]	91	! Initial revision
	92	!
	93	!
	94	!
	95	! Authors:
	96	! --------
[3557]	97	!> @author Eckhard Kadasch (Deutscher Wetterdienst, Offenbach)
[2696]	98	!
	99	! Description:
	100	! ------------
	101	!> The io module contains the functions needed to read and write netCDF data in
	102	!> INIFOR.
	103	!------------------------------------------------------------------------------!
[3680]	104	#if defined ( __netcdf )
[3618]	105	MODULE inifor_io
[2696]	106
[3618]	107	USE inifor_control
	108	USE inifor_defs, &
[3182]	109	ONLY: DATE, SNAME, PATH, PI, dp, hp, TO_RADIANS, TO_DEGREES, VERSION
[3618]	110	USE inifor_types
	111	USE inifor_util, &
[3678]	112	ONLY: add_hours_to, reverse, str, real_to_str
[2696]	113	USE netcdf
	114
	115	IMPLICIT NONE
	116
[3557]	117	!------------------------------------------------------------------------------!
	118	! Description:
	119	! ------------
	120	!> get_netcdf_variable() reads the netCDF data and metadate for the netCDF
	121	!> variable 'in_var % name' from the file 'in_file'. The netCDF data array is
	122	!> stored in the 'buffer' array and metadata added to the respective members of
	123	!> the given 'in_var'.
	124	!------------------------------------------------------------------------------!
[3182]	125	INTERFACE get_netcdf_variable
	126	MODULE PROCEDURE get_netcdf_variable_int
	127	MODULE PROCEDURE get_netcdf_variable_real
	128	END INTERFACE get_netcdf_variable
	129
	130	PRIVATE :: get_netcdf_variable_int, get_netcdf_variable_real
	131
[2696]	132	CONTAINS
	133
[3557]	134	!------------------------------------------------------------------------------!
	135	! Description:
	136	! ------------
	137	!> get_netcdf_variable_int() implements the integer variant for the
	138	!> get_netcdf_variable interface.
	139	!------------------------------------------------------------------------------!
[3182]	140	SUBROUTINE get_netcdf_variable_int(in_file, in_var, buffer)
	141
	142	CHARACTER(LEN=PATH), INTENT(IN) :: in_file
	143	TYPE(nc_var), INTENT(INOUT) :: in_var
	144	INTEGER(hp), ALLOCATABLE, INTENT(INOUT) :: buffer(:,:,:)
	145
[3447]	146	INTEGER :: ncid
	147	INTEGER, DIMENSION(3) :: start, count
[3182]	148
[3447]	149	IF ( nf90_open( TRIM(in_file), NF90_NOWRITE, ncid ) .EQ. NF90_NOERR .AND. &
	150	nf90_inq_varid( ncid, in_var % name, in_var % varid ) .EQ. NF90_NOERR ) THEN
	151
	152	CALL get_input_dimensions(in_var, ncid)
	153
	154	CALL get_netcdf_start_and_count(in_var, start, count)
	155	CALL run_control('time', 'read')
	156
	157	ALLOCATE( buffer( count(1), count(2), count(3) ) )
	158	CALL run_control('time', 'alloc')
	159
	160	CALL check(nf90_get_var( ncid, in_var % varid, buffer, &
	161	start = start, &
	162	count = count ))
	163
	164	ELSE
	165
	166	message = "Failed to read '" // TRIM(in_var % name) // &
	167	"' from file '" // TRIM(in_file) // "'."
[3615]	168	CALL inifor_abort('get_netcdf_variable', message)
[3447]	169
	170	END IF
	171
	172	CALL check(nf90_close(ncid))
	173	CALL run_control('time', 'read')
	174
[3182]	175	END SUBROUTINE get_netcdf_variable_int
	176
	177
[3557]	178	!------------------------------------------------------------------------------!
	179	! Description:
	180	! ------------
	181	!> get_netcdf_variable_real() implements the real variant for the
	182	!> get_netcdf_variable interface.
	183	!------------------------------------------------------------------------------!
[3182]	184	SUBROUTINE get_netcdf_variable_real(in_file, in_var, buffer)
	185
	186	CHARACTER(LEN=PATH), INTENT(IN) :: in_file
	187	TYPE(nc_var), INTENT(INOUT) :: in_var
	188	REAL(dp), ALLOCATABLE, INTENT(INOUT) :: buffer(:,:,:)
	189
[3447]	190	INTEGER :: ncid
	191	INTEGER, DIMENSION(3) :: start, count
[3182]	192
[3447]	193	IF ( nf90_open( TRIM(in_file), NF90_NOWRITE, ncid ) .EQ. NF90_NOERR .AND. &
	194	nf90_inq_varid( ncid, in_var % name, in_var % varid ) .EQ. NF90_NOERR ) THEN
	195
	196	CALL get_input_dimensions(in_var, ncid)
	197
	198	CALL get_netcdf_start_and_count(in_var, start, count)
	199	CALL run_control('time', 'read')
	200
	201	ALLOCATE( buffer( count(1), count(2), count(3) ) )
	202	CALL run_control('time', 'alloc')
	203
	204	CALL check(nf90_get_var( ncid, in_var % varid, buffer, &
	205	start = start, &
	206	count = count ))
	207
	208	ELSE
	209
	210	message = "Failed to read '" // TRIM(in_var % name) // &
	211	"' from file '" // TRIM(in_file) // "'."
[3615]	212	CALL inifor_abort('get_netcdf_variable', message)
[3447]	213
	214	END IF
	215
	216	CALL check(nf90_close(ncid))
	217	CALL run_control('time', 'read')
	218
[3182]	219	END SUBROUTINE get_netcdf_variable_real
	220
	221
[3557]	222	!------------------------------------------------------------------------------!
	223	! Description:
	224	! ------------
	225	!> get_netcdf_dim_vector() reads the coordinate array 'coordname' from the
	226	!> netCDF file 'filename'.
	227	!------------------------------------------------------------------------------!
	228	SUBROUTINE get_netcdf_dim_vector(filename, coordname, coords)
[3537]	229
	230	CHARACTER(LEN=*), INTENT(IN) :: filename
[3557]	231	CHARACTER(LEN=*), INTENT(IN) :: coordname
	232	REAL(dp), ALLOCATABLE, INTENT(INOUT) :: coords(:)
[3537]	233
	234	INTEGER :: ncid, varid, dimlen
	235	INTEGER :: dimids(NF90_MAX_VAR_DIMS)
	236
	237	IF ( nf90_open( TRIM(filename), NF90_NOWRITE, ncid ) .EQ. NF90_NOERR .AND. &
[3557]	238	nf90_inq_varid( ncid, coordname, varid ) .EQ. NF90_NOERR ) THEN
[3537]	239
	240	CALL check(nf90_inquire_variable( ncid, varid, dimids = dimids ))
	241	CALL check(nf90_inquire_dimension( ncid, dimids(1), len = dimlen ))
	242
[3557]	243	ALLOCATE(coords(dimlen))
	244	CALL check(nf90_get_var( ncid, varid, coords))
[3537]	245
	246	ELSE
	247
[3557]	248	message = "Failed to read '" // TRIM(coordname) // &
[3537]	249	"' from file '" // TRIM(filename) // "'."
[3615]	250	CALL inifor_abort('get_netcdf_dim_vector', message)
[3537]	251
	252	END IF
	253
	254	END SUBROUTINE get_netcdf_dim_vector
	255
	256
[3557]	257	!------------------------------------------------------------------------------!
	258	! Description:
	259	! ------------
	260	!> get_input_dimensions() reads dimensions metadata of the netCDF variable given
	261	!> by 'in_var % name' into 'in_var' data structure.
	262	!------------------------------------------------------------------------------!
[3447]	263	SUBROUTINE get_input_dimensions(in_var, ncid)
	264
[3534]	265	TYPE(nc_var), INTENT(INOUT) :: in_var
	266	INTEGER, INTENT(IN) :: ncid
[3447]	267
	268	INTEGER :: i
	269
	270	CALL check(nf90_get_att( ncid, in_var % varid, "long_name", &
	271	in_var % long_name))
	272
	273	CALL check(nf90_get_att( ncid, in_var % varid, "units", &
	274	in_var % units ))
	275
	276	CALL check(nf90_inquire_variable( ncid, in_var % varid, &
	277	ndims = in_var % ndim, &
	278	dimids = in_var % dimids ))
	279
	280	DO i = 1, in_var % ndim
	281	CALL check(nf90_inquire_dimension( ncid, in_var % dimids(i), &
	282	name = in_var % dimname(i), &
	283	len = in_var % dimlen(i) ))
	284	END DO
	285
	286	END SUBROUTINE get_input_dimensions
	287
	288
[3557]	289	!------------------------------------------------------------------------------!
	290	! Description:
	291	! ------------
	292	!> get_netcdf_start_and_count() gets the start position and element counts for
	293	!> the given netCDF file. This information is used in get_netcdf_variable_int()
	294	!> and _real() for reading input variables..
	295	!------------------------------------------------------------------------------!
[3447]	296	SUBROUTINE get_netcdf_start_and_count(in_var, start, count)
	297
	298	TYPE(nc_var), INTENT(INOUT) :: in_var
	299	INTEGER, DIMENSION(3), INTENT(OUT) :: start, count
	300
	301	INTEGER :: ndim
	302
	303	IF ( in_var % ndim .LT. 2 .OR. in_var % ndim .GT. 4 ) THEN
	304
	305	message = "Failed reading NetCDF variable " // &
[3537]	306	TRIM(in_var % name) // " with " // TRIM(str(in_var % ndim)) // &
[3447]	307	" dimensions because only two- and and three-dimensional" // &
	308	" variables are supported."
[3615]	309	CALL inifor_abort('get_netcdf_start_and_count', message)
[3447]	310
	311	END IF
	312
	313	start = (/ 1, 1, 1 /)
	314	IF ( TRIM(in_var % name) .EQ. 'T_SO' ) THEN
[3557]	315	!
	316	!-- Skip depth = 0.0 for T_SO and reduce number of depths from 9 to 8
[3447]	317	in_var % dimlen(3) = in_var % dimlen(3) - 1
	318
[3557]	319	!
	320	!-- Start reading from second level, e.g. depth = 0.005 instead of 0.0
[3447]	321	start(3) = 2
	322	END IF
	323
	324	IF (in_var % ndim .EQ. 2) THEN
	325	in_var % dimlen(3) = 1
	326	ENDIF
	327
	328	ndim = MIN(in_var % ndim, 3)
	329	count = (/ 1, 1, 1 /)
	330	count(1:ndim) = in_var % dimlen(1:ndim)
	331
	332	END SUBROUTINE get_netcdf_start_and_count
	333
	334
[3557]	335	!------------------------------------------------------------------------------!
	336	! Description:
	337	! ------------
	338	!> Routine for defining netCDF variables in the dynamic driver, INIFOR's netCDF
	339	!> output.
	340	!------------------------------------------------------------------------------!
[2696]	341	SUBROUTINE netcdf_define_variable(var, ncid)
	342
	343	TYPE(nc_var), INTENT(INOUT) :: var
	344	INTEGER, INTENT(IN) :: ncid
	345
	346	CALL check(nf90_def_var(ncid, var % name, NF90_FLOAT, var % dimids(1:var % ndim), var % varid))
	347	CALL check(nf90_put_att(ncid, var % varid, "long_name", var % long_name))
	348	CALL check(nf90_put_att(ncid, var % varid, "units", var % units))
[3182]	349	IF ( var % lod .GE. 0 ) THEN
	350	CALL check(nf90_put_att(ncid, var % varid, "lod", var % lod))
	351	END IF
[2696]	352	CALL check(nf90_put_att(ncid, var % varid, "source", var % source))
	353	CALL check(nf90_put_att(ncid, var % varid, "_FillValue", NF90_FILL_REAL))
	354
	355	END SUBROUTINE netcdf_define_variable
	356
	357
[3557]	358	!------------------------------------------------------------------------------!
	359	! Description:
	360	! ------------
	361	!> netcdf_get_dimensions() reads in all dimensions and their lengths and stores
	362	!> them in the given the 'var' data structure. This information is used later
	363	!> for writing output variables in update_output().
	364	!------------------------------------------------------------------------------!
[2696]	365	SUBROUTINE netcdf_get_dimensions(var, ncid)
	366
	367	TYPE(nc_var), INTENT(INOUT) :: var
	368	INTEGER, INTENT(IN) :: ncid
	369	INTEGER :: i
	370	CHARACTER(SNAME) :: null
	371
	372	DO i = 1, var % ndim
	373	CALL check(nf90_inquire_dimension(ncid, var % dimids(i), &
	374	name = null, &
	375	len = var % dimlen(i) ) )
	376	END DO
	377
	378	END SUBROUTINE netcdf_get_dimensions
	379
	380
	381	!------------------------------------------------------------------------------!
	382	! Description:
	383	! ------------
[3557]	384	!> This routine parses and interpretes the command-line options and stores the
	385	!> resulting settings in the 'cfg' data structure.
[2696]	386	!------------------------------------------------------------------------------!
[3182]	387	SUBROUTINE parse_command_line_arguments( cfg )
[2696]	388
[3182]	389	TYPE(inifor_config), INTENT(INOUT) :: cfg
[2696]	390
[3182]	391	CHARACTER(LEN=PATH) :: option, arg
	392	INTEGER :: arg_count, i
[2696]	393
[3395]	394	cfg % p0_is_set = .FALSE.
	395	cfg % ug_is_set = .FALSE.
	396	cfg % vg_is_set = .FALSE.
	397	cfg % flow_prefix_is_set = .FALSE.
	398	cfg % input_prefix_is_set = .FALSE.
	399	cfg % radiation_prefix_is_set = .FALSE.
	400	cfg % soil_prefix_is_set = .FALSE.
	401	cfg % soilmoisture_prefix_is_set = .FALSE.
	402
[2696]	403	arg_count = COMMAND_ARGUMENT_COUNT()
	404	IF (arg_count .GT. 0) THEN
	405
	406	message = "The -clon and -clat command line options are depricated. " // &
	407	"Please remove them form your inifor command and specify the " // &
	408	"location of the PALM-4U origin either" // NEW_LINE(' ') // &
[3182]	409	" - by setting the namelist parameters 'longitude' and 'latitude', or" // NEW_LINE(' ') // &
[2696]	410	" - by providing a static driver netCDF file via the -static command-line option."
	411
[3182]	412	i = 1
	413	DO WHILE (i .LE. arg_count)
[2696]	414
	415	CALL GET_COMMAND_ARGUMENT( i, option )
	416
	417	SELECT CASE( TRIM(option) )
	418
[3395]	419	CASE( '--averaging-mode' )
	420	CALL get_option_argument( i, arg )
	421	cfg % averaging_mode = TRIM(arg)
	422
[3182]	423	CASE( '-date', '-d', '--date' )
	424	CALL get_option_argument( i, arg )
	425	cfg % start_date = TRIM(arg)
[2696]	426
[3395]	427	CASE( '--debug' )
	428	cfg % debug = .TRUE.
	429
[3182]	430	CASE( '-z0', '-z', '--elevation' )
	431	CALL get_option_argument( i, arg )
	432	READ(arg, *) cfg % z0
[2696]	433
[3182]	434	CASE( '-p0', '-r', '--surface-pressure' )
[3395]	435	cfg % p0_is_set = .TRUE.
[3182]	436	CALL get_option_argument( i, arg )
	437	READ(arg, *) cfg % p0
[2696]	438
[3182]	439	CASE( '-ug', '-u', '--geostrophic-u' )
[3395]	440	cfg % ug_is_set = .TRUE.
[3182]	441	CALL get_option_argument( i, arg )
	442	READ(arg, *) cfg % ug
[2696]	443
[3182]	444	CASE( '-vg', '-v', '--geostrophic-v' )
[3395]	445	cfg % vg_is_set = .TRUE.
[3182]	446	CALL get_option_argument( i, arg )
	447	READ(arg, *) cfg % vg
[2696]	448
[3182]	449	CASE( '-clon', '-clat' )
[3615]	450	CALL inifor_abort('parse_command_line_arguments', message)
[2696]	451
[3182]	452	CASE( '-path', '-p', '--path' )
	453	CALL get_option_argument( i, arg )
	454	cfg % input_path = TRIM(arg)
[2696]	455
[3182]	456	CASE( '-hhl', '-l', '--hhl-file' )
	457	CALL get_option_argument( i, arg )
[3395]	458	cfg % hhl_file = TRIM(arg)
[2696]	459
[3395]	460	CASE( '--input-prefix')
	461	CALL get_option_argument( i, arg )
	462	cfg % input_prefix = TRIM(arg)
	463	cfg % input_prefix_is_set = .TRUE.
	464
	465	CASE( '-a', '--averaging-angle' )
	466	CALL get_option_argument( i, arg )
	467	READ(arg, *) cfg % averaging_angle
	468
[3182]	469	CASE( '-static', '-t', '--static-driver' )
	470	CALL get_option_argument( i, arg )
[3395]	471	cfg % static_driver_file = TRIM(arg)
[2696]	472
[3182]	473	CASE( '-soil', '-s', '--soil-file')
	474	CALL get_option_argument( i, arg )
[3395]	475	cfg % soiltyp_file = TRIM(arg)
[2696]	476
[3182]	477	CASE( '--flow-prefix')
	478	CALL get_option_argument( i, arg )
[3395]	479	cfg % flow_prefix = TRIM(arg)
	480	cfg % flow_prefix_is_set = .TRUE.
	481
[3182]	482	CASE( '--radiation-prefix')
	483	CALL get_option_argument( i, arg )
[3395]	484	cfg % radiation_prefix = TRIM(arg)
	485	cfg % radiation_prefix_is_set = .TRUE.
	486
[3182]	487	CASE( '--soil-prefix')
	488	CALL get_option_argument( i, arg )
[3395]	489	cfg % soil_prefix = TRIM(arg)
	490	cfg % soil_prefix_is_set = .TRUE.
	491
[3182]	492	CASE( '--soilmoisture-prefix')
	493	CALL get_option_argument( i, arg )
[3395]	494	cfg % soilmoisture_prefix = TRIM(arg)
	495	cfg % soilmoisture_prefix_is_set = .TRUE.
[2696]	496
[3182]	497	CASE( '-o', '--output' )
	498	CALL get_option_argument( i, arg )
	499	cfg % output_file = TRIM(arg)
[2696]	500
[3182]	501	CASE( '-n', '--namelist' )
	502	CALL get_option_argument( i, arg )
	503	cfg % namelist_file = TRIM(arg)
[2696]	504
[3182]	505	CASE( '-mode', '-i', '--init-mode' )
	506	CALL get_option_argument( i, arg )
	507	cfg % ic_mode = TRIM(arg)
	508
	509	CASE( '-f', '--forcing-mode' )
	510	CALL get_option_argument( i, arg )
	511	cfg % bc_mode = TRIM(arg)
	512
	513	CASE( '--version' )
	514	CALL print_version()
	515	STOP
	516
	517	CASE( '--help' )
	518	CALL print_version()
	519	PRINT *, ""
	520	PRINT *, "For a list of command-line options have a look at the README file."
	521	STOP
	522
[2696]	523	CASE DEFAULT
[3182]	524	message = "unknown option '" // TRIM(option) // "'."
[3615]	525	CALL inifor_abort('parse_command_line_arguments', message)
[2696]	526
	527	END SELECT
	528
[3182]	529	i = i + 1
	530
[2696]	531	END DO
	532
	533	ELSE
	534
	535	message = "No arguments present, using default input and output files"
	536	CALL report('parse_command_line_arguments', message)
	537
	538	END IF
	539
	540	END SUBROUTINE parse_command_line_arguments
	541
[3182]	542
[3678]	543
	544	SUBROUTINE get_datetime_file_list( start_date_string, start_hour, end_hour, &
	545	step_hour, input_path, prefix, suffix, &
	546	file_list )
	547
	548	CHARACTER (LEN=DATE), INTENT(IN) :: start_date_string
	549	CHARACTER (LEN=*), INTENT(IN) :: prefix, suffix, input_path
	550	INTEGER, INTENT(IN) :: start_hour, end_hour, step_hour
	551	CHARACTER(LEN=*), ALLOCATABLE, INTENT(INOUT) :: file_list(:)
	552
	553	INTEGER :: number_of_intervals, hour, i
	554	CHARACTER(LEN=DATE) :: date_string
	555
	556	number_of_intervals = CEILING( REAL(end_hour - start_hour) / step_hour )
	557	ALLOCATE( file_list(number_of_intervals + 1) )
	558
	559	DO i = 0, number_of_intervals
	560
	561	hour = start_hour + i * step_hour
	562	date_string = add_hours_to(start_date_string, hour)
	563
	564	file_list(i+1) = TRIM(input_path) // TRIM(prefix) // &
	565	TRIM(date_string) // TRIM(suffix) // '.nc'
	566
	567	message = "Set up input file name '" // TRIM(file_list(i+1)) // "'"
	568	CALL report('input_file_list', message)
	569
	570	END DO
	571
	572	END SUBROUTINE get_datetime_file_list
	573
[3557]	574	!------------------------------------------------------------------------------!
	575	! Description:
	576	! ------------
[3678]	577	!> Establish a list of files based on the given start and end times and file
	578	!> prefixes and suffixes.
	579	!------------------------------------------------------------------------------!
	580	SUBROUTINE get_input_file_list( start_date_string, start_hour, end_hour, &
	581	step_hour, input_path, prefix, suffix, &
[3716]	582	file_list, nocheck )
[3678]	583
	584	CHARACTER (LEN=DATE), INTENT(IN) :: start_date_string
	585	CHARACTER (LEN=*), INTENT(IN) :: prefix, suffix, input_path
	586	INTEGER, INTENT(IN) :: start_hour, end_hour, step_hour
	587	CHARACTER(LEN=*), ALLOCATABLE, INTENT(INOUT) :: file_list(:)
[3716]	588	LOGICAL, OPTIONAL, INTENT(IN) :: nocheck
[3678]	589
	590	INTEGER :: number_of_intervals, hour, i
	591	CHARACTER(LEN=DATE) :: date_string
	592	CHARACTER(LEN=PATH) :: file_name
[3716]	593	LOGICAL :: check_files
[3678]	594
	595	CALL get_datetime_file_list( start_date_string, start_hour, end_hour, &
	596	step_hour, input_path, prefix, suffix, &
	597	file_list )
	598
[3716]	599	check_files = .TRUE.
	600	IF ( PRESENT ( nocheck ) ) THEN
	601	IF ( nocheck ) check_files = .FALSE.
	602	END IF
[3678]	603
[3716]	604	IF ( check_files ) THEN
[3678]	605
[3716]	606	tip = "Please check if you specified the correct file prefix " // &
	607	"using the options --input-prefix, --flow-prefix, etc."
	608
	609	DO i = 1, SIZE(file_list)
	610	CALL verify_file(file_list(i), 'input', tip)
	611	END DO
	612
	613	END IF
	614
[3678]	615	END SUBROUTINE get_input_file_list
	616
	617
	618	!------------------------------------------------------------------------------!
	619	! Description:
	620	! ------------
	621	!> Abort INIFOR if the given file is not present.
	622	!------------------------------------------------------------------------------!
	623	SUBROUTINE verify_file(file_name, file_kind, tip)
	624
	625	CHARACTER(LEN=*), INTENT(IN) :: file_name, file_kind
	626	CHARACTER(LEN=*), INTENT(IN), OPTIONAL :: tip
	627
	628	IF (.NOT. file_present(file_name)) THEN
	629
	630	IF (LEN(TRIM(file_name)) == 0) THEN
	631
	632	message = "No name was given for the " // TRIM(file_kind) // " file."
	633
	634	ELSE
	635
	636	message = "The " // TRIM(file_kind) // " file '" // &
	637	TRIM(file_name) // "' was not found."
	638
	639	IF (PRESENT(tip)) THEN
	640	message = TRIM(message) // " " // TRIM(tip)
	641	END IF
	642
	643	END IF
	644
	645	CALL inifor_abort('verify_file', message)
	646
	647	END IF
	648
	649	END SUBROUTINE verify_file
	650
	651
	652	!------------------------------------------------------------------------------!
	653	! Description:
	654	! ------------
[3557]	655	!> Get the argument of the i'th command line option, which is at the location
	656	!> i+1 of the argument list.
	657	!------------------------------------------------------------------------------!
[3182]	658	SUBROUTINE get_option_argument(i, arg)
	659	CHARACTER(LEN=PATH), INTENT(INOUT) :: arg
	660	INTEGER, INTENT(INOUT) :: i
[2696]	661
[3182]	662	i = i + 1
	663	CALL GET_COMMAND_ARGUMENT(i, arg)
	664
	665	END SUBROUTINE
	666
	667
[3557]	668	!------------------------------------------------------------------------------!
	669	! Description:
	670	! ------------
	671	!> Checks the INIFOR configuration 'cfg' for plausibility.
	672	!------------------------------------------------------------------------------!
[3182]	673	SUBROUTINE validate_config(cfg)
	674	TYPE(inifor_config), INTENT(IN) :: cfg
	675
[3678]	676	CALL verify_file(cfg % hhl_file, 'HHL')
	677	CALL verify_file(cfg % namelist_file, 'NAMELIST')
	678	CALL verify_file(cfg % soiltyp_file, 'SOILTYP')
[3182]	679
[3557]	680	!
	681	!-- Only check optional static driver file name, if it has been given.
[3182]	682	IF (TRIM(cfg % static_driver_file) .NE. '') THEN
[3678]	683	CALL verify_file(cfg % static_driver_file, 'static driver')
[3182]	684	END IF
	685
	686	SELECT CASE( TRIM(cfg % ic_mode) )
	687	CASE( 'profile', 'volume')
	688	CASE DEFAULT
	689	message = "Initialization mode '" // TRIM(cfg % ic_mode) //&
	690	"' is not supported. " //&
	691	"Please select either 'profile' or 'volume', " //&
	692	"or omit the -i/--init-mode/-mode option entirely, which corresponds "//&
	693	"to the latter."
[3615]	694	CALL inifor_abort( 'validate_config', message )
[3182]	695	END SELECT
	696
	697
	698	SELECT CASE( TRIM(cfg % bc_mode) )
	699	CASE( 'real', 'ideal')
	700	CASE DEFAULT
	701	message = "Forcing mode '" // TRIM(cfg % bc_mode) //&
	702	"' is not supported. " //&
	703	"Please select either 'real' or 'ideal', " //&
	704	"or omit the -f/--forcing-mode option entirely, which corresponds "//&
	705	"to the latter."
[3615]	706	CALL inifor_abort( 'validate_config', message )
[3182]	707	END SELECT
	708
[3395]	709	SELECT CASE( TRIM(cfg % averaging_mode) )
	710	CASE( 'level', 'height')
	711	CASE DEFAULT
	712	message = "Averaging mode '" // TRIM(cfg % averaging_mode) //&
	713	"' is not supported. " //&
	714	"Please select either 'height' or 'level', " //&
	715	"or omit the --averaging-mode option entirely, which corresponds "//&
	716	"to the latter."
[3615]	717	CALL inifor_abort( 'validate_config', message )
[3395]	718	END SELECT
[3182]	719
[3395]	720	IF ( cfg % ug_is_set .NEQV. cfg % vg_is_set ) THEN
	721	message = "You specified only one component of the geostrophic " // &
	722	"wind. Please specify either both or none."
[3615]	723	CALL inifor_abort( 'validate_config', message )
[3395]	724	END IF
	725
[3182]	726	END SUBROUTINE validate_config
	727
	728
[3557]	729	!------------------------------------------------------------------------------!
	730	! Description:
	731	! ------------
	732	!> Check whether the given file is present on the filesystem.
	733	!------------------------------------------------------------------------------!
[3678]	734	LOGICAL FUNCTION file_present(filename)
[3182]	735	CHARACTER(LEN=PATH), INTENT(IN) :: filename
	736
[3678]	737	INQUIRE(FILE=filename, EXIST=file_present)
[3182]	738
	739	END FUNCTION file_present
	740
	741
[2696]	742	!------------------------------------------------------------------------------!
	743	! Description:
	744	! ------------
[3557]	745	!> This routine initializes the dynamic driver file, i.e. INIFOR's netCDF output
	746	!> file.
[2696]	747	!>
	748	!> Besides writing metadata, such as global attributes, coordinates, variables,
[3557]	749	!> in the netCDF file, various netCDF IDs are saved for later, when INIFOR
[2696]	750	!> writes the actual data.
	751	!------------------------------------------------------------------------------!
[3182]	752	SUBROUTINE setup_netcdf_dimensions(output_file, palm_grid, &
	753	start_date_string, origin_lon, origin_lat)
[2696]	754
	755	TYPE(nc_file), INTENT(INOUT) :: output_file
	756	TYPE(grid_definition), INTENT(IN) :: palm_grid
[3182]	757	CHARACTER (LEN=DATE), INTENT(IN) :: start_date_string
	758	REAL(dp), INTENT(IN) :: origin_lon, origin_lat
[2696]	759
[3182]	760	CHARACTER (LEN=8) :: date_string
	761	CHARACTER (LEN=10) :: time_string
	762	CHARACTER (LEN=5) :: zone_string
	763	CHARACTER (LEN=SNAME) :: history_string
[2696]	764	INTEGER :: ncid, nx, ny, nz, nt, dimids(3), dimvarids(3)
	765	REAL(dp) :: z0
	766
[3182]	767	message = "Initializing PALM-4U dynamic driver file '" // &
	768	TRIM(output_file % name) // "' and setting up dimensions."
	769	CALL report('setup_netcdf_dimensions', message)
	770
[3557]	771	!
	772	!-- Create the netCDF file as in netCDF-4/HDF5 format if __netcdf4 preprocessor flag is given
[3182]	773	#if defined( __netcdf4 )
[2696]	774	CALL check(nf90_create(TRIM(output_file % name), OR(NF90_CLOBBER, NF90_HDF5), ncid))
[3182]	775	#else
	776	CALL check(nf90_create(TRIM(output_file % name), NF90_CLOBBER, ncid))
	777	#endif
[2696]	778
[3395]	779	!------------------------------------------------------------------------------
	780	!- Section 1: Define NetCDF dimensions and coordinates
	781	!------------------------------------------------------------------------------
	782	nt = SIZE(output_file % time)
	783	nx = palm_grid % nx
	784	ny = palm_grid % ny
	785	nz = palm_grid % nz
	786	z0 = palm_grid % z0
	787
	788
[2696]	789	!
	790	!------------------------------------------------------------------------------
[3395]	791	!- Section 2: Write global NetCDF attributes
[2696]	792	!------------------------------------------------------------------------------
[3182]	793	CALL date_and_time(DATE=date_string, TIME=time_string, ZONE=zone_string)
	794	history_string = &
	795	'Created on '// date_string // &
	796	' at ' // time_string(1:2) // ':' // time_string(3:4) // &
	797	' (UTC' // zone_string // ')'
	798
[2696]	799	CALL check(nf90_put_att(ncid, NF90_GLOBAL, 'title', 'PALM input file for scenario ...'))
	800	CALL check(nf90_put_att(ncid, NF90_GLOBAL, 'institution', 'Deutscher Wetterdienst, Offenbach'))
	801	CALL check(nf90_put_att(ncid, NF90_GLOBAL, 'author', 'Eckhard Kadasch, eckhard.kadasch@dwd.de'))
[3182]	802	CALL check(nf90_put_att(ncid, NF90_GLOBAL, 'history', TRIM(history_string)))
[2696]	803	CALL check(nf90_put_att(ncid, NF90_GLOBAL, 'references', '--'))
	804	CALL check(nf90_put_att(ncid, NF90_GLOBAL, 'comment', '--'))
[3182]	805	CALL check(nf90_put_att(ncid, NF90_GLOBAL, 'origin_lat', TRIM(real_to_str(origin_lat*TO_DEGREES, '(F18.13)'))))
	806	CALL check(nf90_put_att(ncid, NF90_GLOBAL, 'origin_lon', TRIM(real_to_str(origin_lon*TO_DEGREES, '(F18.13)'))))
[3557]	807	!
	808	!-- FIXME: This is the elevation relative to COSMO-DE/D2 sea level and does
	809	!-- FIXME: not necessarily comply with DHHN2016 (c.f. PALM Input Data
	810	!-- FIXME: Standard v1.9., origin_z)
[3395]	811	CALL check(nf90_put_att(ncid, NF90_GLOBAL, 'origin_z', TRIM(real_to_str(z0, '(F18.13)'))))
[2696]	812	CALL check(nf90_put_att(ncid, NF90_GLOBAL, 'inifor_version', TRIM(VERSION)))
	813	CALL check(nf90_put_att(ncid, NF90_GLOBAL, 'palm_version', '--'))
	814
	815	!
	816	!
	817	!------------------------------------------------------------------------------
	818	!- Section 2a: Define dimensions for cell centers (scalars in soil and atmosph.)
	819	!------------------------------------------------------------------------------
[3557]	820	!
	821	!-- reset dimids first
	822	dimids = (/0, 0, 0/)
	823	CALL check( nf90_def_dim(ncid, "x", nx+1, dimids(1)) )
	824	CALL check( nf90_def_dim(ncid, "y", ny+1, dimids(2)) )
	825	CALL check( nf90_def_dim(ncid, "z", nz, dimids(3)) )
	826	!
	827	!-- save dimids for later
	828	output_file % dimids_scl = dimids
[2696]	829
[3557]	830	!
	831	!-- reset dimvarids first
	832	dimvarids = (/0, 0, 0/)
	833	CALL check(nf90_def_var(ncid, "x", NF90_FLOAT, dimids(1), dimvarids(1)))
	834	CALL check(nf90_put_att(ncid, dimvarids(1), "standard_name", "x coordinate of cell centers"))
	835	CALL check(nf90_put_att(ncid, dimvarids(1), "units", "m"))
[2696]	836
[3557]	837	CALL check(nf90_def_var(ncid, "y", NF90_FLOAT, dimids(2), dimvarids(2)))
	838	CALL check(nf90_put_att(ncid, dimvarids(2), "standard_name", "y coordinate of cell centers"))
	839	CALL check(nf90_put_att(ncid, dimvarids(2), "units", "m"))
[2696]	840
[3557]	841	CALL check(nf90_def_var(ncid, "z", NF90_FLOAT, dimids(3), dimvarids(3)))
	842	CALL check(nf90_put_att(ncid, dimvarids(3), "standard_name", "z coordinate of cell centers"))
	843	CALL check(nf90_put_att(ncid, dimvarids(3), "units", "m"))
	844	!
	845	!-- save dimvarids for later
	846	output_file % dimvarids_scl = dimvarids
[2696]	847
[3557]	848	!
	849	!-- overwrite third dimid with the one of depth
[3182]	850	CALL check(nf90_def_dim(ncid, "zsoil", SIZE(palm_grid % depths), dimids(3)) )
[3557]	851	!
	852	!-- save dimids for later
	853	output_file % dimids_soil = dimids
[2696]	854
[3557]	855	!
	856	!-- overwrite third dimvarid with the one of depth
[3182]	857	CALL check(nf90_def_var(ncid, "zsoil", NF90_FLOAT, output_file % dimids_soil(3), dimvarids(3)))
[2696]	858	CALL check(nf90_put_att(ncid, dimvarids(3), "standard_name", "depth_below_land"))
	859	CALL check(nf90_put_att(ncid, dimvarids(3), "positive", "down"))
	860	CALL check(nf90_put_att(ncid, dimvarids(3), "units", "m"))
	861	!
[3557]	862	!-- save dimvarids for later
	863	output_file % dimvarids_soil = dimvarids
	864	!
[2696]	865	!------------------------------------------------------------------------------
	866	!- Section 2b: Define dimensions for cell faces/velocities
	867	!------------------------------------------------------------------------------
[3557]	868	!
	869	!-- reset dimids first
	870	dimids = (/0, 0, 0/)
	871	CALL check(nf90_def_dim(ncid, "xu", nx, dimids(1)) )
	872	CALL check(nf90_def_dim(ncid, "yv", ny, dimids(2)) )
	873	CALL check(nf90_def_dim(ncid, "zw", nz-1, dimids(3)) )
	874	!
	875	!-- save dimids for later
	876	output_file % dimids_vel = dimids
[2696]	877
[3557]	878	!
	879	!-- reset dimvarids first
	880	dimvarids = (/0, 0, 0/)
	881	CALL check(nf90_def_var(ncid, "xu", NF90_FLOAT, dimids(1), dimvarids(1)))
	882	CALL check(nf90_put_att(ncid, dimvarids(1), "standard_name", "x coordinate of cell faces"))
	883	CALL check(nf90_put_att(ncid, dimvarids(1), "units", "m"))
[2696]	884
[3557]	885	CALL check(nf90_def_var(ncid, "yv", NF90_FLOAT, dimids(2), dimvarids(2)))
	886	CALL check(nf90_put_att(ncid, dimvarids(2), "standard_name", "y coordinate of cell faces"))
	887	CALL check(nf90_put_att(ncid, dimvarids(2), "units", "m"))
[2696]	888
[3557]	889	CALL check(nf90_def_var(ncid, "zw", NF90_FLOAT, dimids(3), dimvarids(3)))
	890	CALL check(nf90_put_att(ncid, dimvarids(3), "standard_name", "z coordinate of cell faces"))
	891	CALL check(nf90_put_att(ncid, dimvarids(3), "units", "m"))
	892	!
	893	!-- save dimvarids for later
	894	output_file % dimvarids_vel = dimvarids
[2696]	895
	896	!
	897	!------------------------------------------------------------------------------
	898	!- Section 2c: Define time dimension
	899	!------------------------------------------------------------------------------
	900	CALL check(nf90_def_dim(ncid, "time", nt, output_file % dimid_time) )
	901	CALL check(nf90_def_var(ncid, "time", NF90_FLOAT, &
	902	output_file % dimid_time, &
	903	output_file % dimvarid_time))
	904	CALL check(nf90_put_att(ncid, output_file % dimvarid_time, "standard_name", "time"))
	905	CALL check(nf90_put_att(ncid, output_file % dimvarid_time, "long_name", "time"))
[3182]	906	CALL check(nf90_put_att(ncid, output_file % dimvarid_time, "units", &
	907	"seconds since " // start_date_string // " UTC"))
[2696]	908
	909	CALL check(nf90_enddef(ncid))
	910
	911	!
	912	!------------------------------------------------------------------------------
	913	!- Section 3: Write grid coordinates
	914	!------------------------------------------------------------------------------
	915	CALL check(nf90_put_var(ncid, output_file % dimvarids_scl(1), palm_grid%x))
	916	CALL check(nf90_put_var(ncid, output_file % dimvarids_scl(2), palm_grid%y))
	917	CALL check(nf90_put_var(ncid, output_file % dimvarids_scl(3), palm_grid%z))
	918
	919	CALL check(nf90_put_var(ncid, output_file % dimvarids_vel(1), palm_grid%xu))
	920	CALL check(nf90_put_var(ncid, output_file % dimvarids_vel(2), palm_grid%yv))
	921	CALL check(nf90_put_var(ncid, output_file % dimvarids_vel(3), palm_grid%zw))
	922
[3557]	923	!
	924	!-- TODO Read in soil depths from input file before this.
[2696]	925	CALL check(nf90_put_var(ncid, output_file % dimvarids_soil(3), palm_grid%depths))
	926
[3557]	927	!
	928	!-- Write time vector
[2696]	929	CALL check(nf90_put_var(ncid, output_file % dimvarid_time, output_file % time))
	930
[3557]	931	!
	932	!-- Close the file
[2696]	933	CALL check(nf90_close(ncid))
	934
	935	END SUBROUTINE setup_netcdf_dimensions
	936
	937
[3557]	938	!------------------------------------------------------------------------------!
	939	! Description:
	940	! ------------
	941	!> Defines the netCDF variables to be written to the dynamic driver file
	942	!------------------------------------------------------------------------------!
[3456]	943	SUBROUTINE setup_netcdf_variables(filename, output_variable_table, debug)
[2696]	944
	945	CHARACTER (LEN=*), INTENT(IN) :: filename
	946	TYPE(nc_var), INTENT(INOUT), TARGET :: output_variable_table(:)
[3456]	947	LOGICAL, INTENT(IN) :: debug
	948
[2696]	949	TYPE(nc_var), POINTER :: var
	950	INTEGER :: i, ncid
[3456]	951	LOGICAL :: to_be_written
[2696]	952
[3182]	953	message = "Defining variables in dynamic driver '" // TRIM(filename) // "'."
[2696]	954	CALL report('setup_netcdf_variables', message)
	955
	956	CALL check(nf90_open(TRIM(filename), NF90_WRITE, ncid))
	957	CALL check(nf90_redef(ncid))
	958
	959	DO i = 1, SIZE(output_variable_table)
	960
	961	var => output_variable_table(i)
	962
[3456]	963	to_be_written = ( var % to_be_processed .AND. .NOT. var % is_internal) .OR. &
	964	( var % is_internal .AND. debug )
	965
	966	IF ( to_be_written ) THEN
[3182]	967	message = " variable #" // TRIM(str(i)) // " '" // TRIM(var%name) // "'."
[2696]	968	CALL report('setup_netcdf_variables', message)
	969
	970	CALL netcdf_define_variable(var, ncid)
	971	CALL netcdf_get_dimensions(var, ncid)
	972	END IF
	973
	974	END DO
	975
	976	CALL check(nf90_enddef(ncid))
	977	CALL check(nf90_close(ncid))
	978
[3182]	979	message = "Dynamic driver '" // TRIM(filename) // "' initialized successfully."
[2696]	980	CALL report('setup_netcdf_variables', message)
	981
	982	END SUBROUTINE setup_netcdf_variables
	983
	984
	985	!------------------------------------------------------------------------------!
	986	! Description:
	987	! ------------
	988	!> This routine reads and returns all input variables of the given IO group
	989	!> It accomodates the data by allocating a container variable that represents a
	990	!> list of arrays of the same length as the groups variable list. (This list
	991	!> will typically contain one or two items.) After the container, its members
	992	!> are allocated one by one with the appropriate, possibly different,
	993	!> dimensions.
	994	!>
	995	!> The 'group' is an INTENT(INOUT) variable so that 'get_netcdf_variable()' can
	996	!> record netCDF IDs in the 'in_var_list()' member variable.
	997	!------------------------------------------------------------------------------!
	998	SUBROUTINE read_input_variables(group, iter, buffer)
	999	TYPE(io_group), INTENT(INOUT), TARGET :: group
	1000	INTEGER, INTENT(IN) :: iter
	1001	TYPE(container), ALLOCATABLE, INTENT(INOUT) :: buffer(:)
	1002	INTEGER :: hour, buf_id
	1003	TYPE(nc_var), POINTER :: input_var
	1004	CHARACTER(LEN=PATH), POINTER :: input_file
	1005	INTEGER :: ivar, nbuffers
	1006
	1007	message = "Reading data for I/O group '" // TRIM(group % in_var_list(1) % name) // "'."
	1008	CALL report('read_input_variables', message)
	1009
	1010	input_file => group % in_files(iter)
	1011
	1012	!
	1013	!------------------------------------------------------------------------------
	1014	!- Section 1: Load input buffers for accumulated variables
	1015	!------------------------------------------------------------------------------
[3557]	1016	!
	1017	!-- radiation budgets, precipitation
[2696]	1018	IF (group % kind == 'running average' .OR. &
[3557]	1019	group % kind == 'accumulated') THEN
[2696]	1020
	1021	IF (SIZE(group % in_var_list) .GT. 1 ) THEN
	1022	message = "I/O groups may not contain more than one " // &
	1023	"accumulated variable. Group '" // TRIM(group % kind) //&
	1024	"' contains " // &
	1025	TRIM( str(SIZE(group % in_var_list)) ) // "."
[3615]	1026	CALL inifor_abort('read_input_variables \| accumulation', message)
[2696]	1027	END IF
	1028
[3557]	1029	!
	1030	!-- use two buffer arrays
[2696]	1031	nbuffers = 2
	1032	IF ( .NOT. ALLOCATED( buffer ) ) ALLOCATE( buffer(nbuffers) )
	1033
[3557]	1034	!
	1035	!-- hour of the day
	1036	hour = iter - 1
	1037	!
	1038	!-- chose correct buffer array
[2696]	1039	buf_id = select_buffer(hour)
	1040
	1041	CALL run_control('time', 'read')
	1042	IF ( ALLOCATED(buffer(buf_id) % array) ) DEALLOCATE(buffer(buf_id) % array)
	1043	CALL run_control('time', 'alloc')
	1044
	1045	input_var => group % in_var_list(1)
[3182]	1046	CALL get_netcdf_variable(input_file, input_var, buffer(buf_id) % array)
[2696]	1047	CALL report('read_input_variables', "Read accumulated " // TRIM(group % in_var_list(1) % name))
	1048
	1049	IF ( input_var % is_upside_down ) CALL reverse(buffer(buf_id) % array)
	1050	CALL run_control('time', 'comp')
	1051
	1052	!------------------------------------------------------------------------------
	1053	!- Section 2: Load input buffers for normal I/O groups
	1054	!------------------------------------------------------------------------------
	1055	ELSE
	1056
[3557]	1057	!
	1058	!-- Allocate one input buffer per input_variable. If more quantities
	1059	!-- have to be computed than input variables exist in this group,
	1060	!-- allocate more buffers. For instance, in the thermodynamics group,
	1061	!-- there are three input variabels (PP, T, Qv) and four quantities
	1062	!-- necessart (P, Theta, Rho, qv) for the corresponding output fields
	1063	!-- (p0, Theta, qv, ug, and vg)
[3395]	1064	nbuffers = MAX( group % n_inputs, group % n_output_quantities )
[2696]	1065	ALLOCATE( buffer(nbuffers) )
	1066	CALL run_control('time', 'alloc')
	1067
[3557]	1068	!
	1069	!-- Read in all input variables, leave extra buffers-if any-untouched.
[3395]	1070	DO ivar = 1, group % n_inputs
[2696]	1071
	1072	input_var => group % in_var_list(ivar)
	1073
[3557]	1074	!
	1075	! Check wheather P or PP is present in input file
[2696]	1076	IF (input_var % name == 'P') THEN
[3395]	1077	input_var % name = TRIM( get_pressure_varname(input_file) )
[2696]	1078	CALL run_control('time', 'read')
	1079	END IF
	1080
[3182]	1081	CALL get_netcdf_variable(input_file, input_var, buffer(ivar) % array)
[2696]	1082
	1083	IF ( input_var % is_upside_down ) CALL reverse(buffer(ivar) % array)
	1084	CALL run_control('time', 'comp')
	1085
	1086	END DO
	1087	END IF
	1088
	1089	END SUBROUTINE read_input_variables
	1090
	1091
[3557]	1092	!------------------------------------------------------------------------------!
	1093	! Description:
	1094	! ------------
	1095	!> Select the appropriate buffer ID for accumulated COSMO input variables
	1096	!> depending on the current hour.
	1097	!------------------------------------------------------------------------------!
[2696]	1098	INTEGER FUNCTION select_buffer(hour)
	1099	INTEGER, INTENT(IN) :: hour
	1100	INTEGER :: step
	1101
	1102	select_buffer = 0
	1103	step = MODULO(hour, 3) + 1
	1104
	1105	SELECT CASE(step)
	1106	CASE(1, 3)
	1107	select_buffer = 1
	1108	CASE(2)
	1109	select_buffer = 2
	1110	CASE DEFAULT
	1111	message = "Invalid step '" // TRIM(str(step))
[3615]	1112	CALL inifor_abort('select_buffer', message)
[2696]	1113	END SELECT
	1114	END FUNCTION select_buffer
	1115
	1116
	1117	!------------------------------------------------------------------------------!
	1118	! Description:
	1119	! ------------
	1120	!> Checks if the input_file contains the absolute pressure, 'P', or the pressure
	1121	!> perturbation, 'PP', and returns the appropriate string.
	1122	!------------------------------------------------------------------------------!
[3395]	1123	CHARACTER(LEN=2) FUNCTION get_pressure_varname(input_file) RESULT(var)
[2696]	1124	CHARACTER(LEN=*) :: input_file
	1125	INTEGER :: ncid, varid
	1126
	1127	CALL check(nf90_open( TRIM(input_file), NF90_NOWRITE, ncid ))
	1128	IF ( nf90_inq_varid( ncid, 'P', varid ) .EQ. NF90_NOERR ) THEN
	1129
	1130	var = 'P'
	1131
	1132	ELSE IF ( nf90_inq_varid( ncid, 'PP', varid ) .EQ. NF90_NOERR ) THEN
	1133
	1134	var = 'PP'
	1135	CALL report('get_pressure_var', 'Using PP instead of P')
	1136
	1137	ELSE
	1138
	1139	message = "Failed to read '" // TRIM(var) // &
	1140	"' from file '" // TRIM(input_file) // "'."
[3615]	1141	CALL inifor_abort('get_pressure_var', message)
[2696]	1142
	1143	END IF
	1144
	1145	CALL check(nf90_close(ncid))
	1146
[3395]	1147	END FUNCTION get_pressure_varname
[2696]	1148
	1149
[3557]	1150	!------------------------------------------------------------------------------!
	1151	! Description:
	1152	! ------------
	1153	!> Read the given global attribute form the given netCDF file.
	1154	!------------------------------------------------------------------------------!
[2696]	1155	FUNCTION get_netcdf_attribute(filename, attribute) RESULT(attribute_value)
	1156
	1157	CHARACTER(LEN=*), INTENT(IN) :: filename, attribute
	1158	REAL(dp) :: attribute_value
	1159
[3557]	1160	INTEGER :: ncid
[2696]	1161
	1162	IF ( nf90_open( TRIM(filename), NF90_NOWRITE, ncid ) == NF90_NOERR ) THEN
	1163
	1164	CALL check(nf90_get_att(ncid, NF90_GLOBAL, TRIM(attribute), attribute_value))
[3182]	1165	CALL check(nf90_close(ncid))
[2696]	1166
	1167	ELSE
	1168
	1169	message = "Failed to read '" // TRIM(attribute) // &
	1170	"' from file '" // TRIM(filename) // "'."
[3615]	1171	CALL inifor_abort('get_netcdf_attribute', message)
[2696]	1172
	1173	END IF
	1174
	1175	END FUNCTION get_netcdf_attribute
	1176
	1177
[3557]	1178
	1179	!------------------------------------------------------------------------------!
	1180	! Description:
	1181	! ------------
	1182	!> Updates the dynamic driver with the interpolated field of the current
	1183	!> variable at the current time step.
	1184	!------------------------------------------------------------------------------!
[3456]	1185	SUBROUTINE update_output(var, array, iter, output_file, cfg)
[2696]	1186	TYPE(nc_var), INTENT(IN) :: var
	1187	REAL(dp), INTENT(IN) :: array(:,:,:)
	1188	INTEGER, INTENT(IN) :: iter
	1189	TYPE(nc_file), INTENT(IN) :: output_file
[3456]	1190	TYPE(inifor_config) :: cfg
[2696]	1191
	1192	INTEGER :: ncid, ndim, start(4), count(4)
	1193	LOGICAL :: var_is_time_dependent
	1194
	1195	var_is_time_dependent = ( &
	1196	var % dimids( var % ndim ) == output_file % dimid_time &
	1197	)
	1198
[3557]	1199	!
	1200	!-- Skip time dimension for output
[3182]	1201	ndim = var % ndim
	1202	IF ( var_is_time_dependent ) ndim = var % ndim - 1
[2696]	1203
	1204	start(:) = (/1,1,1,1/)
	1205	start(ndim+1) = iter
	1206	count(1:ndim) = var%dimlen(1:ndim)
	1207
	1208	CALL check(nf90_open(output_file % name, NF90_WRITE, ncid))
	1209
[3557]	1210	!
	1211	!-- Reduce dimension of output array according to variable kind
[2696]	1212	SELECT CASE (TRIM(var % kind))
	1213
	1214	CASE ( 'init scalar profile', 'init u profile', 'init v profile', &
	1215	'init w profile' )
	1216
	1217	CALL check(nf90_put_var( ncid, var%varid, array(1,1,:) ) )
	1218
	1219	CASE ( 'init soil', 'init scalar', 'init u', 'init v', 'init w' )
	1220
	1221	CALL check(nf90_put_var( ncid, var%varid, array(:,:,:) ) )
	1222
	1223	CASE( 'left scalar', 'right scalar', 'left w', 'right w' )
	1224
	1225	CALL check(nf90_put_var( ncid, var%varid, array(1,:,:), &
	1226	start=start(1:ndim+1), &
	1227	count=count(1:ndim) ) )
	1228
	1229
	1230	IF (.NOT. SIZE(array, 2) .EQ. var % dimlen(1)) THEN
	1231	PRINT *, "inifor: update_output: Dimension ", 1, " of variable ", &
	1232	TRIM(var % name), " (", var % dimlen(1), &
	1233	") does not match the dimension of the output array (", &
	1234	SIZE(array, 2), ")."
	1235	STOP
	1236	END IF
	1237
	1238
	1239	CASE( 'north scalar', 'south scalar', 'north w', 'south w' )
	1240
	1241	CALL check(nf90_put_var( ncid, var%varid, array(:,1,:), &
	1242	start=start(1:ndim+1), &
	1243	count=count(1:ndim) ) )
	1244
	1245
	1246	CASE( 'surface forcing', 'top scalar', 'top w' )
	1247
	1248	CALL check(nf90_put_var( ncid, var%varid, array(:,:,1), &
	1249	start=start(1:ndim+1), &
	1250	count=count(1:ndim) ) )
	1251
	1252	CASE ( 'left u', 'right u', 'left v', 'right v' )
	1253
	1254	CALL check(nf90_put_var( ncid, var%varid, array(1,:,:), &
	1255	start=start(1:ndim+1), &
	1256	count=count(1:ndim) ) )
	1257
	1258	CASE ( 'north u', 'south u', 'north v', 'south v' )
	1259
	1260	CALL check(nf90_put_var( ncid, var%varid, array(:,1,:), &
	1261	start=start(1:ndim+1), &
	1262	count=count(1:ndim) ) )
	1263
	1264	CASE ( 'top u', 'top v' )
	1265
	1266	CALL check(nf90_put_var( ncid, var%varid, array(:,:,1), &
	1267	start=start(1:ndim+1), &
	1268	count=count(1:ndim) ) )
	1269
	1270	CASE ( 'time series' )
	1271
	1272	CALL check(nf90_put_var( ncid, var%varid, array(1,1,1), &
	1273	start=start(1:ndim+1) ) )
	1274
[3456]	1275	CASE ( 'constant scalar profile', 'geostrophic' )
[2696]	1276
	1277	CALL check(nf90_put_var( ncid, var%varid, array(1,1,:), &
	1278	start=start(1:ndim+1), &
	1279	count=count(1:ndim) ) )
	1280
[3456]	1281	CASE ( 'internal profile' )
	1282
	1283	IF ( cfg % debug ) THEN
	1284	CALL check(nf90_put_var( ncid, var%varid, array(1,1,:), &
	1285	start=start(1:ndim+1), &
	1286	count=count(1:ndim) ) )
	1287	END IF
	1288
[3182]	1289	CASE ( 'large-scale scalar forcing', 'large-scale w forcing' )
	1290
	1291	message = "Doing nothing in terms of writing large-scale forings."
	1292	CALL report('update_output', message)
	1293
[2696]	1294	CASE DEFAULT
	1295
	1296	message = "Variable kind '" // TRIM(var % kind) // &
	1297	"' not recognized."
[3615]	1298	CALL inifor_abort('update_output', message)
[2696]	1299
	1300	END SELECT
	1301
	1302	CALL check(nf90_close(ncid))
	1303
	1304	END SUBROUTINE update_output
	1305
	1306
[3557]	1307	!------------------------------------------------------------------------------!
	1308	! Description:
	1309	! ------------
	1310	!> Checks the status of a netCDF API call and aborts if an error occured
	1311	!------------------------------------------------------------------------------!
[2696]	1312	SUBROUTINE check(status)
	1313
	1314	INTEGER, INTENT(IN) :: status
	1315
	1316	IF (status /= nf90_noerr) THEN
	1317	message = "NetCDF API call failed with error: " // &
	1318	TRIM( nf90_strerror(status) )
[3615]	1319	CALL inifor_abort('io.check', message)
[2696]	1320	END IF
	1321
	1322	END SUBROUTINE check
	1323
[3618]	1324	END MODULE inifor_io
[3680]	1325	#endif

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |