Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |