Home

Context Navigation

source: palm/trunk/UTIL/inifor/src/io.f90 @ 3359

Last change on this file since 3359 was 3262, checked in by eckhard, 6 years ago
Removed unnecessary file check
Property svn:keywords set to `Id`
File size: 34.1 KB

Rev	Line
[2696]	1	!> @file src/io.f90
	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: io.f90 3262 2018-09-19 15:47:15Z knoop $
[3262]	28	! Removed unnecessary check for output file
	29	!
	30	! 3183 2018-07-27 14:25:55Z suehring
[3182]	31	! Introduced new PALM grid stretching
	32	! Updated variable names and metadata for PIDS v1.9 compatibility
	33	! Improved handling of the start date string
	34	! Better compatibility with older Intel compilers:
	35	! - avoiding implicit array allocation with new get_netcdf_variable()
	36	! subroutine instead of function
	37	! Improved command line interface:
	38	! - Added configuration validation
	39	! - New options to configure input file prefixes
	40	! - GNU-style short and long option names
	41	! - Added version and copyright output
[2696]	42	!
[3182]	43	!
[3183]	44	! 3182 2018-07-27 13:36:03Z suehring
[2696]	45	! Initial revision
	46	!
	47	!
	48	!
	49	! Authors:
	50	! --------
	51	! @author Eckhard Kadasch
	52	!
	53	! Description:
	54	! ------------
	55	!> The io module contains the functions needed to read and write netCDF data in
	56	!> INIFOR.
	57	!------------------------------------------------------------------------------!
	58	MODULE io
	59
	60	USE control
	61	USE defs, &
[3182]	62	ONLY: DATE, SNAME, PATH, PI, dp, hp, TO_RADIANS, TO_DEGREES, VERSION
[2696]	63	USE netcdf
	64	USE types
	65	USE util, &
[3182]	66	ONLY: reverse, str, real_to_str
[2696]	67
	68	IMPLICIT NONE
	69
[3182]	70	INTERFACE get_netcdf_variable
	71	MODULE PROCEDURE get_netcdf_variable_int
	72	MODULE PROCEDURE get_netcdf_variable_real
	73	END INTERFACE get_netcdf_variable
	74
	75	PRIVATE :: get_netcdf_variable_int, get_netcdf_variable_real
	76
[2696]	77	CONTAINS
	78
[3182]	79	SUBROUTINE get_netcdf_variable_int(in_file, in_var, buffer)
	80
	81	CHARACTER(LEN=PATH), INTENT(IN) :: in_file
	82	TYPE(nc_var), INTENT(INOUT) :: in_var
	83	INTEGER(hp), ALLOCATABLE, INTENT(INOUT) :: buffer(:,:,:)
	84
	85	INCLUDE 'get_netcdf_variable.inc'
	86
	87	END SUBROUTINE get_netcdf_variable_int
	88
	89
	90	SUBROUTINE get_netcdf_variable_real(in_file, in_var, buffer)
	91
	92	CHARACTER(LEN=PATH), INTENT(IN) :: in_file
	93	TYPE(nc_var), INTENT(INOUT) :: in_var
	94	REAL(dp), ALLOCATABLE, INTENT(INOUT) :: buffer(:,:,:)
	95
	96	INCLUDE 'get_netcdf_variable.inc'
	97
	98	END SUBROUTINE get_netcdf_variable_real
	99
	100
[2696]	101	SUBROUTINE netcdf_define_variable(var, ncid)
	102
	103	TYPE(nc_var), INTENT(INOUT) :: var
	104	INTEGER, INTENT(IN) :: ncid
	105
	106	CALL check(nf90_def_var(ncid, var % name, NF90_FLOAT, var % dimids(1:var % ndim), var % varid))
	107	CALL check(nf90_put_att(ncid, var % varid, "long_name", var % long_name))
	108	CALL check(nf90_put_att(ncid, var % varid, "units", var % units))
[3182]	109	IF ( var % lod .GE. 0 ) THEN
	110	CALL check(nf90_put_att(ncid, var % varid, "lod", var % lod))
	111	END IF
[2696]	112	CALL check(nf90_put_att(ncid, var % varid, "source", var % source))
	113	CALL check(nf90_put_att(ncid, var % varid, "_FillValue", NF90_FILL_REAL))
	114
	115	END SUBROUTINE netcdf_define_variable
	116
	117
	118	SUBROUTINE netcdf_get_dimensions(var, ncid)
	119
	120	TYPE(nc_var), INTENT(INOUT) :: var
	121	INTEGER, INTENT(IN) :: ncid
	122	INTEGER :: i
	123	CHARACTER(SNAME) :: null
	124
	125	! Remember dimension lenghts for NetCDF writing routine
	126	DO i = 1, var % ndim
	127	CALL check(nf90_inquire_dimension(ncid, var % dimids(i), &
	128	name = null, &
	129	len = var % dimlen(i) ) )
	130	END DO
	131
	132	END SUBROUTINE netcdf_get_dimensions
	133
	134
	135	!------------------------------------------------------------------------------!
	136	! Description:
	137	! ------------
	138	!> This routine initializes Inifor. This includes parsing command-line
	139	!> arguments, setting the names of the input and output file names as well as
	140	!> the name of the input namelist and, subsequently, reading in and setting grid
	141	!> parameters for the PALM-4U computational grid.
	142	!------------------------------------------------------------------------------!
[3182]	143	SUBROUTINE parse_command_line_arguments( cfg )
[2696]	144
[3182]	145	TYPE(inifor_config), INTENT(INOUT) :: cfg
[2696]	146
[3182]	147	CHARACTER(LEN=PATH) :: option, arg
	148	INTEGER :: arg_count, i
[2696]	149
	150	arg_count = COMMAND_ARGUMENT_COUNT()
	151	IF (arg_count .GT. 0) THEN
	152
	153	! Every option should have an argument.
[3182]	154	!IF ( MOD(arg_count, 2) .NE. 0 ) THEN
	155	! message = "Syntax error in command line."
	156	! CALL abort('parse_command_line_arguments', message)
	157	!END IF
[2696]	158
	159	message = "The -clon and -clat command line options are depricated. " // &
	160	"Please remove them form your inifor command and specify the " // &
	161	"location of the PALM-4U origin either" // NEW_LINE(' ') // &
[3182]	162	" - by setting the namelist parameters 'longitude' and 'latitude', or" // NEW_LINE(' ') // &
[2696]	163	" - by providing a static driver netCDF file via the -static command-line option."
	164
[3182]	165	i = 1
	166	DO WHILE (i .LE. arg_count)
[2696]	167
	168	CALL GET_COMMAND_ARGUMENT( i, option )
	169
	170	SELECT CASE( TRIM(option) )
	171
[3182]	172	CASE( '-date', '-d', '--date' )
	173	CALL get_option_argument( i, arg )
	174	cfg % start_date = TRIM(arg)
[2696]	175
	176	! Elevation of the PALM-4U domain above sea level
[3182]	177	CASE( '-z0', '-z', '--elevation' )
	178	CALL get_option_argument( i, arg )
	179	READ(arg, *) cfg % z0
[2696]	180
	181	! surface pressure, at z0
[3182]	182	CASE( '-p0', '-r', '--surface-pressure' )
	183	CALL get_option_argument( i, arg )
	184	READ(arg, *) cfg % p0
[2696]	185
[3182]	186	! geostrophic wind in x direction
	187	CASE( '-ug', '-u', '--geostrophic-u' )
	188	CALL get_option_argument( i, arg )
	189	READ(arg, *) cfg % ug
[2696]	190
[3182]	191	! geostrophic wind in y direction
	192	CASE( '-vg', '-v', '--geostrophic-v' )
	193	CALL get_option_argument( i, arg )
	194	READ(arg, *) cfg % vg
[2696]	195
[3182]	196	! domain centre geographical longitude and latitude
	197	CASE( '-clon', '-clat' )
[2696]	198	CALL abort('parse_command_line_arguments', message)
	199	!READ(arg, *) lambda_cg
	200	!lambda_cg = lambda_cg * TO_RADIANS
	201	!READ(arg, *) phi_cg
	202	!phi_cg = phi_cg * TO_RADIANS
	203
[3182]	204	CASE( '-path', '-p', '--path' )
	205	CALL get_option_argument( i, arg )
	206	cfg % input_path = TRIM(arg)
[2696]	207
[3182]	208	CASE( '-hhl', '-l', '--hhl-file' )
	209	CALL get_option_argument( i, arg )
	210	cfg % hhl_file = TRIM(arg)
[2696]	211
[3182]	212	CASE( '-static', '-t', '--static-driver' )
	213	CALL get_option_argument( i, arg )
	214	cfg % static_driver_file = TRIM(arg)
[2696]	215
[3182]	216	CASE( '-soil', '-s', '--soil-file')
	217	CALL get_option_argument( i, arg )
	218	cfg % soiltyp_file = TRIM(arg)
[2696]	219
[3182]	220	CASE( '--flow-prefix')
	221	CALL get_option_argument( i, arg )
	222	cfg % flow_prefix = TRIM(arg)
[2696]	223
[3182]	224	CASE( '--radiation-prefix')
	225	CALL get_option_argument( i, arg )
	226	cfg % radiation_prefix = TRIM(arg)
[2696]	227
[3182]	228	CASE( '--soil-prefix')
	229	CALL get_option_argument( i, arg )
	230	cfg % soil_prefix = TRIM(arg)
[2696]	231
[3182]	232	CASE( '--soilmoisture-prefix')
	233	CALL get_option_argument( i, arg )
	234	cfg % soilmoisture_prefix = TRIM(arg)
[2696]	235
[3182]	236	CASE( '-o', '--output' )
	237	CALL get_option_argument( i, arg )
	238	cfg % output_file = TRIM(arg)
[2696]	239
[3182]	240	CASE( '-n', '--namelist' )
	241	CALL get_option_argument( i, arg )
	242	cfg % namelist_file = TRIM(arg)
[2696]	243
[3182]	244	! initial condition mode: 'profile' / 'volume'
	245	CASE( '-mode', '-i', '--init-mode' )
	246	CALL get_option_argument( i, arg )
	247	cfg % ic_mode = TRIM(arg)
	248
	249	! boundary conditions / forcing mode: 'ideal' / 'real'
	250	CASE( '-f', '--forcing-mode' )
	251	CALL get_option_argument( i, arg )
	252	cfg % bc_mode = TRIM(arg)
	253
	254	CASE( '--version' )
	255	CALL print_version()
	256	STOP
	257
	258	CASE( '--help' )
	259	CALL print_version()
	260	PRINT *, ""
	261	PRINT *, "For a list of command-line options have a look at the README file."
	262	STOP
	263
[2696]	264	CASE DEFAULT
[3182]	265	message = "unknown option '" // TRIM(option) // "'."
[2696]	266	CALL abort('parse_command_line_arguments', message)
	267
	268	END SELECT
	269
[3182]	270	i = i + 1
	271
[2696]	272	END DO
	273
	274	ELSE
	275
	276	message = "No arguments present, using default input and output files"
	277	CALL report('parse_command_line_arguments', message)
	278
	279	END IF
	280
	281	END SUBROUTINE parse_command_line_arguments
	282
[3182]	283
	284	SUBROUTINE get_option_argument(i, arg)
	285	CHARACTER(LEN=PATH), INTENT(INOUT) :: arg
	286	INTEGER, INTENT(INOUT) :: i
[2696]	287
[3182]	288	i = i + 1
	289	CALL GET_COMMAND_ARGUMENT(i, arg)
	290
	291	END SUBROUTINE
	292
	293
	294	SUBROUTINE validate_config(cfg)
	295	TYPE(inifor_config), INTENT(IN) :: cfg
	296	LOGICAL :: all_files_present
	297
	298	all_files_present = .TRUE.
	299	all_files_present = all_files_present .AND. file_present(cfg % hhl_file)
	300	all_files_present = all_files_present .AND. file_present(cfg % namelist_file)
	301	all_files_present = all_files_present .AND. file_present(cfg % soiltyp_file)
	302
	303	! Only check optional static driver file name, if it has been given.
	304	IF (TRIM(cfg % static_driver_file) .NE. '') THEN
	305	all_files_present = all_files_present .AND. file_present(cfg % static_driver_file)
	306	END IF
	307
	308	IF (.NOT. all_files_present) THEN
	309	message = "INIFOR configuration invalid; some input files are missing."
	310	CALL abort( 'validate_config', message )
	311	END IF
	312
	313
	314	SELECT CASE( TRIM(cfg % ic_mode) )
	315	CASE( 'profile', 'volume')
	316	CASE DEFAULT
	317	message = "Initialization mode '" // TRIM(cfg % ic_mode) //&
	318	"' is not supported. " //&
	319	"Please select either 'profile' or 'volume', " //&
	320	"or omit the -i/--init-mode/-mode option entirely, which corresponds "//&
	321	"to the latter."
	322	CALL abort( 'validate_config', message )
	323	END SELECT
	324
	325
	326	SELECT CASE( TRIM(cfg % bc_mode) )
	327	CASE( 'real', 'ideal')
	328	CASE DEFAULT
	329	message = "Forcing mode '" // TRIM(cfg % bc_mode) //&
	330	"' is not supported. " //&
	331	"Please select either 'real' or 'ideal', " //&
	332	"or omit the -f/--forcing-mode option entirely, which corresponds "//&
	333	"to the latter."
	334	CALL abort( 'validate_config', message )
	335	END SELECT
	336
	337
	338	END SUBROUTINE validate_config
	339
	340
	341	LOGICAL FUNCTION file_present(filename)
	342	CHARACTER(LEN=PATH), INTENT(IN) :: filename
	343
	344	INQUIRE(FILE=filename, EXIST=file_present)
	345
	346	IF (.NOT. file_present) THEN
	347	message = "The given file '" // "' does not exist."
	348	CALL report('file_present', message)
	349	END IF
	350
	351	END FUNCTION file_present
	352
	353
[2696]	354	!------------------------------------------------------------------------------!
	355	! Description:
	356	! ------------
	357	!> This routine initializes the Inifor output file, i.e. the PALM-4U
	358	!> initializing and forcing data as a NetCDF file.
	359	!>
	360	!> Besides writing metadata, such as global attributes, coordinates, variables,
	361	!> in the NetCDF file, various NetCDF IDs are saved for later, when Inifor
	362	!> writes the actual data.
	363	!------------------------------------------------------------------------------!
[3182]	364	SUBROUTINE setup_netcdf_dimensions(output_file, palm_grid, &
	365	start_date_string, origin_lon, origin_lat)
[2696]	366
	367	TYPE(nc_file), INTENT(INOUT) :: output_file
	368	TYPE(grid_definition), INTENT(IN) :: palm_grid
[3182]	369	CHARACTER (LEN=DATE), INTENT(IN) :: start_date_string
	370	REAL(dp), INTENT(IN) :: origin_lon, origin_lat
[2696]	371
[3182]	372	CHARACTER (LEN=8) :: date_string
	373	CHARACTER (LEN=10) :: time_string
	374	CHARACTER (LEN=5) :: zone_string
	375	CHARACTER (LEN=SNAME) :: history_string
[2696]	376	INTEGER :: ncid, nx, ny, nz, nt, dimids(3), dimvarids(3)
	377	REAL(dp) :: z0
	378
[3182]	379	message = "Initializing PALM-4U dynamic driver file '" // &
	380	TRIM(output_file % name) // "' and setting up dimensions."
	381	CALL report('setup_netcdf_dimensions', message)
	382
[2696]	383	! Create the NetCDF file. NF90_CLOBBER selects overwrite mode.
[3182]	384	#if defined( __netcdf4 )
[2696]	385	CALL check(nf90_create(TRIM(output_file % name), OR(NF90_CLOBBER, NF90_HDF5), ncid))
[3182]	386	#else
	387	CALL check(nf90_create(TRIM(output_file % name), NF90_CLOBBER, ncid))
	388	#endif
[2696]	389
	390	!
	391	!------------------------------------------------------------------------------
	392	!- Section 1: Write global NetCDF attributes
	393	!------------------------------------------------------------------------------
[3182]	394	CALL date_and_time(DATE=date_string, TIME=time_string, ZONE=zone_string)
	395	history_string = &
	396	'Created on '// date_string // &
	397	' at ' // time_string(1:2) // ':' // time_string(3:4) // &
	398	' (UTC' // zone_string // ')'
	399
[2696]	400	CALL check(nf90_put_att(ncid, NF90_GLOBAL, 'title', 'PALM input file for scenario ...'))
	401	CALL check(nf90_put_att(ncid, NF90_GLOBAL, 'institution', 'Deutscher Wetterdienst, Offenbach'))
	402	CALL check(nf90_put_att(ncid, NF90_GLOBAL, 'author', 'Eckhard Kadasch, eckhard.kadasch@dwd.de'))
[3182]	403	CALL check(nf90_put_att(ncid, NF90_GLOBAL, 'history', TRIM(history_string)))
[2696]	404	CALL check(nf90_put_att(ncid, NF90_GLOBAL, 'references', '--'))
	405	CALL check(nf90_put_att(ncid, NF90_GLOBAL, 'comment', '--'))
[3182]	406	CALL check(nf90_put_att(ncid, NF90_GLOBAL, 'origin_lat', TRIM(real_to_str(origin_lat*TO_DEGREES, '(F18.13)'))))
	407	CALL check(nf90_put_att(ncid, NF90_GLOBAL, 'origin_lon', TRIM(real_to_str(origin_lon*TO_DEGREES, '(F18.13)'))))
[2696]	408	CALL check(nf90_put_att(ncid, NF90_GLOBAL, 'inifor_version', TRIM(VERSION)))
	409	CALL check(nf90_put_att(ncid, NF90_GLOBAL, 'palm_version', '--'))
	410
	411	!
	412	!------------------------------------------------------------------------------
	413	!- Section 2: Define NetCDF dimensions and coordinates
	414	!------------------------------------------------------------------------------
	415	nt = SIZE(output_file % time)
	416	nx = palm_grid % nx
	417	ny = palm_grid % ny
	418	nz = palm_grid % nz
	419	z0 = palm_grid % z0
	420
	421	!
	422	!------------------------------------------------------------------------------
	423	!- Section 2a: Define dimensions for cell centers (scalars in soil and atmosph.)
	424	!------------------------------------------------------------------------------
	425	dimids = (/0, 0, 0/) ! reset dimids
	426	CALL check( nf90_def_dim(ncid, "x", nx+1, dimids(1)) )
	427	CALL check( nf90_def_dim(ncid, "y", ny+1, dimids(2)) )
[3182]	428	CALL check( nf90_def_dim(ncid, "z", nz, dimids(3)) )
[2696]	429	output_file % dimids_scl = dimids ! save dimids for later
	430
	431	dimvarids = (/0, 0, 0/) ! reset dimvarids
	432	CALL check(nf90_def_var(ncid, "x", NF90_FLOAT, dimids(1), dimvarids(1)))
	433	CALL check(nf90_put_att(ncid, dimvarids(1), "standard_name", "x coordinate of cell centers"))
	434	CALL check(nf90_put_att(ncid, dimvarids(1), "units", "m"))
	435
	436	CALL check(nf90_def_var(ncid, "y", NF90_FLOAT, dimids(2), dimvarids(2)))
	437	CALL check(nf90_put_att(ncid, dimvarids(2), "standard_name", "y coordinate of cell centers"))
	438	CALL check(nf90_put_att(ncid, dimvarids(2), "units", "m"))
	439
	440	CALL check(nf90_def_var(ncid, "z", NF90_FLOAT, dimids(3), dimvarids(3)))
	441	CALL check(nf90_put_att(ncid, dimvarids(3), "standard_name", "z coordinate of cell centers"))
	442	CALL check(nf90_put_att(ncid, dimvarids(3), "units", "m"))
	443	output_file % dimvarids_scl = dimvarids ! save dimvarids for later
	444
	445	! overwrite third dimid with the one of depth
[3182]	446	CALL check(nf90_def_dim(ncid, "zsoil", SIZE(palm_grid % depths), dimids(3)) )
[2696]	447	output_file % dimids_soil = dimids ! save dimids for later
	448
	449	! overwrite third dimvarid with the one of depth
[3182]	450	CALL check(nf90_def_var(ncid, "zsoil", NF90_FLOAT, output_file % dimids_soil(3), dimvarids(3)))
[2696]	451	CALL check(nf90_put_att(ncid, dimvarids(3), "standard_name", "depth_below_land"))
	452	CALL check(nf90_put_att(ncid, dimvarids(3), "positive", "down"))
	453	CALL check(nf90_put_att(ncid, dimvarids(3), "units", "m"))
	454	output_file % dimvarids_soil = dimvarids ! save dimvarids for later
	455	!
	456	!------------------------------------------------------------------------------
	457	!- Section 2b: Define dimensions for cell faces/velocities
	458	!------------------------------------------------------------------------------
	459	dimids = (/0, 0, 0/) ! reset dimids
	460	CALL check(nf90_def_dim(ncid, "xu", nx, dimids(1)) )
	461	CALL check(nf90_def_dim(ncid, "yv", ny, dimids(2)) )
[3182]	462	CALL check(nf90_def_dim(ncid, "zw", nz-1, dimids(3)) )
[2696]	463	output_file % dimids_vel = dimids ! save dimids for later
	464
	465	dimvarids = (/0, 0, 0/) ! reset dimvarids
	466	CALL check(nf90_def_var(ncid, "xu", NF90_FLOAT, dimids(1), dimvarids(1)))
	467	CALL check(nf90_put_att(ncid, dimvarids(1), "standard_name", "x coordinate of cell faces"))
	468	CALL check(nf90_put_att(ncid, dimvarids(1), "units", "m"))
	469
	470	CALL check(nf90_def_var(ncid, "yv", NF90_FLOAT, dimids(2), dimvarids(2)))
	471	CALL check(nf90_put_att(ncid, dimvarids(2), "standard_name", "y coordinate of cell faces"))
	472	CALL check(nf90_put_att(ncid, dimvarids(2), "units", "m"))
	473
	474	CALL check(nf90_def_var(ncid, "zw", NF90_FLOAT, dimids(3), dimvarids(3)))
	475	CALL check(nf90_put_att(ncid, dimvarids(3), "standard_name", "z coordinate of cell faces"))
	476	CALL check(nf90_put_att(ncid, dimvarids(3), "units", "m"))
	477	output_file % dimvarids_vel = dimvarids ! save dimvarids for later
	478
	479	!
	480	!------------------------------------------------------------------------------
	481	!- Section 2c: Define time dimension
	482	!------------------------------------------------------------------------------
	483	CALL check(nf90_def_dim(ncid, "time", nt, output_file % dimid_time) )
	484	CALL check(nf90_def_var(ncid, "time", NF90_FLOAT, &
	485	output_file % dimid_time, &
	486	output_file % dimvarid_time))
	487	CALL check(nf90_put_att(ncid, output_file % dimvarid_time, "standard_name", "time"))
	488	CALL check(nf90_put_att(ncid, output_file % dimvarid_time, "long_name", "time"))
[3182]	489	CALL check(nf90_put_att(ncid, output_file % dimvarid_time, "units", &
	490	"seconds since " // start_date_string // " UTC"))
[2696]	491
	492	CALL check(nf90_enddef(ncid))
	493
	494	!
	495	!------------------------------------------------------------------------------
	496	!- Section 3: Write grid coordinates
	497	!------------------------------------------------------------------------------
	498	CALL check(nf90_put_var(ncid, output_file % dimvarids_scl(1), palm_grid%x))
	499	CALL check(nf90_put_var(ncid, output_file % dimvarids_scl(2), palm_grid%y))
	500	CALL check(nf90_put_var(ncid, output_file % dimvarids_scl(3), palm_grid%z))
	501
	502	CALL check(nf90_put_var(ncid, output_file % dimvarids_vel(1), palm_grid%xu))
	503	CALL check(nf90_put_var(ncid, output_file % dimvarids_vel(2), palm_grid%yv))
	504	CALL check(nf90_put_var(ncid, output_file % dimvarids_vel(3), palm_grid%zw))
	505
	506	! TODO Read in soil depths from input file before this.
	507	CALL check(nf90_put_var(ncid, output_file % dimvarids_soil(3), palm_grid%depths))
	508
	509	! Write time vector
	510	CALL check(nf90_put_var(ncid, output_file % dimvarid_time, output_file % time))
	511
	512	! Close the file
	513	CALL check(nf90_close(ncid))
	514
	515	END SUBROUTINE setup_netcdf_dimensions
	516
	517
	518	SUBROUTINE setup_netcdf_variables(filename, output_variable_table)
	519
	520	CHARACTER (LEN=*), INTENT(IN) :: filename
	521	TYPE(nc_var), INTENT(INOUT), TARGET :: output_variable_table(:)
	522	TYPE(nc_var), POINTER :: var
	523	INTEGER :: i, ncid
	524
[3182]	525	message = "Defining variables in dynamic driver '" // TRIM(filename) // "'."
[2696]	526	CALL report('setup_netcdf_variables', message)
	527
	528	CALL check(nf90_open(TRIM(filename), NF90_WRITE, ncid))
	529	CALL check(nf90_redef(ncid))
	530
	531	DO i = 1, SIZE(output_variable_table)
	532
	533	var => output_variable_table(i)
	534
	535	IF ( var % to_be_processed ) THEN
[3182]	536	message = " variable #" // TRIM(str(i)) // " '" // TRIM(var%name) // "'."
[2696]	537	CALL report('setup_netcdf_variables', message)
	538
	539	CALL netcdf_define_variable(var, ncid)
	540	CALL netcdf_get_dimensions(var, ncid)
	541	END IF
	542
	543	END DO
	544
	545	CALL check(nf90_enddef(ncid))
	546	CALL check(nf90_close(ncid))
	547
[3182]	548	message = "Dynamic driver '" // TRIM(filename) // "' initialized successfully."
[2696]	549	CALL report('setup_netcdf_variables', message)
	550
	551	END SUBROUTINE setup_netcdf_variables
	552
	553
	554	!------------------------------------------------------------------------------!
	555	! Description:
	556	! ------------
	557	!> This routine reads and returns all input variables of the given IO group
	558	!> It accomodates the data by allocating a container variable that represents a
	559	!> list of arrays of the same length as the groups variable list. (This list
	560	!> will typically contain one or two items.) After the container, its members
	561	!> are allocated one by one with the appropriate, possibly different,
	562	!> dimensions.
	563	!>
	564	!> The 'group' is an INTENT(INOUT) variable so that 'get_netcdf_variable()' can
	565	!> record netCDF IDs in the 'in_var_list()' member variable.
	566	!------------------------------------------------------------------------------!
	567	SUBROUTINE read_input_variables(group, iter, buffer)
	568	TYPE(io_group), INTENT(INOUT), TARGET :: group
	569	INTEGER, INTENT(IN) :: iter
	570	TYPE(container), ALLOCATABLE, INTENT(INOUT) :: buffer(:)
	571	INTEGER :: hour, buf_id
	572	TYPE(nc_var), POINTER :: input_var
	573	CHARACTER(LEN=PATH), POINTER :: input_file
	574	INTEGER :: ivar, nbuffers
	575
	576	message = "Reading data for I/O group '" // TRIM(group % in_var_list(1) % name) // "'."
	577	CALL report('read_input_variables', message)
	578
	579	input_file => group % in_files(iter)
	580
	581	!
	582	!------------------------------------------------------------------------------
	583	!- Section 1: Load input buffers for accumulated variables
	584	!------------------------------------------------------------------------------
	585	IF (group % kind == 'running average' .OR. &
	586	group % kind == 'accumulated') THEN ! radiation budgets, precipitation
	587
	588	IF (SIZE(group % in_var_list) .GT. 1 ) THEN
	589	message = "I/O groups may not contain more than one " // &
	590	"accumulated variable. Group '" // TRIM(group % kind) //&
	591	"' contains " // &
	592	TRIM( str(SIZE(group % in_var_list)) ) // "."
	593	CALL abort('read_input_variables \| accumulation', message)
	594	END IF
	595
	596	! use two buffer arrays
	597	nbuffers = 2
	598	IF ( .NOT. ALLOCATED( buffer ) ) ALLOCATE( buffer(nbuffers) )
	599
	600	! chose correct buffer array
	601	hour = iter - 1! hour of the day
	602	buf_id = select_buffer(hour)
	603
	604	CALL run_control('time', 'read')
	605	IF ( ALLOCATED(buffer(buf_id) % array) ) DEALLOCATE(buffer(buf_id) % array)
	606	CALL run_control('time', 'alloc')
	607
	608	input_var => group % in_var_list(1)
[3182]	609	CALL get_netcdf_variable(input_file, input_var, buffer(buf_id) % array)
[2696]	610	CALL report('read_input_variables', "Read accumulated " // TRIM(group % in_var_list(1) % name))
	611
	612	IF ( input_var % is_upside_down ) CALL reverse(buffer(buf_id) % array)
	613	CALL run_control('time', 'comp')
	614
	615	!------------------------------------------------------------------------------
	616	!- Section 2: Load input buffers for normal I/O groups
	617	!------------------------------------------------------------------------------
	618	ELSE
	619
	620	nbuffers = SIZE( group % in_var_list )
	621	ALLOCATE( buffer(nbuffers) )
	622	CALL run_control('time', 'alloc')
	623
	624	DO ivar = 1, nbuffers
	625
	626	input_var => group % in_var_list(ivar)
	627
	628	! Check wheather P or PP is present in input file
	629	IF (input_var % name == 'P') THEN
	630	input_var % name = TRIM( get_pressure_var(input_file) )
	631	CALL run_control('time', 'read')
	632	END IF
	633
[3182]	634	CALL get_netcdf_variable(input_file, input_var, buffer(ivar) % array)
[2696]	635
	636	IF ( input_var % is_upside_down ) CALL reverse(buffer(ivar) % array)
	637	CALL run_control('time', 'comp')
	638
	639	END DO
	640	END IF
	641
	642	END SUBROUTINE read_input_variables
	643
	644
	645	INTEGER FUNCTION select_buffer(hour)
	646	INTEGER, INTENT(IN) :: hour
	647	INTEGER :: step
	648
	649	select_buffer = 0
	650	step = MODULO(hour, 3) + 1
	651
	652	SELECT CASE(step)
	653	CASE(1, 3)
	654	select_buffer = 1
	655	CASE(2)
	656	select_buffer = 2
	657	CASE DEFAULT
	658	message = "Invalid step '" // TRIM(str(step))
	659	CALL abort('select_buffer', message)
	660	END SELECT
	661	END FUNCTION select_buffer
	662
	663
	664	!------------------------------------------------------------------------------!
	665	! Description:
	666	! ------------
	667	!> Checks if the input_file contains the absolute pressure, 'P', or the pressure
	668	!> perturbation, 'PP', and returns the appropriate string.
	669	!------------------------------------------------------------------------------!
	670	CHARACTER(LEN=2) FUNCTION get_pressure_var(input_file) RESULT(var)
	671	CHARACTER(LEN=*) :: input_file
	672	INTEGER :: ncid, varid
	673
	674	CALL check(nf90_open( TRIM(input_file), NF90_NOWRITE, ncid ))
	675	IF ( nf90_inq_varid( ncid, 'P', varid ) .EQ. NF90_NOERR ) THEN
	676
	677	var = 'P'
	678
	679	ELSE IF ( nf90_inq_varid( ncid, 'PP', varid ) .EQ. NF90_NOERR ) THEN
	680
	681	var = 'PP'
	682	CALL report('get_pressure_var', 'Using PP instead of P')
	683
	684	ELSE
	685
	686	message = "Failed to read '" // TRIM(var) // &
	687	"' from file '" // TRIM(input_file) // "'."
	688	CALL abort('get_pressure_var', message)
	689
	690	END IF
	691
	692	CALL check(nf90_close(ncid))
	693
	694	END FUNCTION get_pressure_var
	695
	696
	697	FUNCTION get_netcdf_attribute(filename, attribute) RESULT(attribute_value)
	698
	699	CHARACTER(LEN=*), INTENT(IN) :: filename, attribute
	700	REAL(dp) :: attribute_value
	701
	702	INTEGER :: ncid
	703
	704	IF ( nf90_open( TRIM(filename), NF90_NOWRITE, ncid ) == NF90_NOERR ) THEN
	705
	706	CALL check(nf90_get_att(ncid, NF90_GLOBAL, TRIM(attribute), attribute_value))
[3182]	707	CALL check(nf90_close(ncid))
[2696]	708
	709	ELSE
	710
	711	message = "Failed to read '" // TRIM(attribute) // &
	712	"' from file '" // TRIM(filename) // "'."
	713	CALL abort('get_netcdf_attribute', message)
	714
	715	END IF
	716
	717	END FUNCTION get_netcdf_attribute
	718
	719
	720	SUBROUTINE update_output(var, array, iter, output_file)
	721	TYPE(nc_var), INTENT(IN) :: var
	722	REAL(dp), INTENT(IN) :: array(:,:,:)
	723	INTEGER, INTENT(IN) :: iter
	724	TYPE(nc_file), INTENT(IN) :: output_file
	725
	726	INTEGER :: ncid, ndim, start(4), count(4)
	727	LOGICAL :: var_is_time_dependent
	728
	729	var_is_time_dependent = ( &
	730	var % dimids( var % ndim ) == output_file % dimid_time &
	731	)
	732
	733	! Skip time dimension for output
[3182]	734	ndim = var % ndim
	735	IF ( var_is_time_dependent ) ndim = var % ndim - 1
[2696]	736
	737	start(:) = (/1,1,1,1/)
	738	start(ndim+1) = iter
	739	count(1:ndim) = var%dimlen(1:ndim)
	740
	741	CALL check(nf90_open(output_file % name, NF90_WRITE, ncid))
	742
	743	! Reduce dimension of output array according to variable kind
	744	SELECT CASE (TRIM(var % kind))
	745
	746	CASE ( 'init scalar profile', 'init u profile', 'init v profile', &
	747	'init w profile' )
	748
	749	CALL check(nf90_put_var( ncid, var%varid, array(1,1,:) ) )
	750
	751	CASE ( 'init soil', 'init scalar', 'init u', 'init v', 'init w' )
	752
	753	CALL check(nf90_put_var( ncid, var%varid, array(:,:,:) ) )
	754
	755	CASE( 'left scalar', 'right scalar', 'left w', 'right w' )
	756
	757	CALL check(nf90_put_var( ncid, var%varid, array(1,:,:), &
	758	start=start(1:ndim+1), &
	759	count=count(1:ndim) ) )
	760
	761
	762	IF (.NOT. SIZE(array, 2) .EQ. var % dimlen(1)) THEN
	763	PRINT *, "inifor: update_output: Dimension ", 1, " of variable ", &
	764	TRIM(var % name), " (", var % dimlen(1), &
	765	") does not match the dimension of the output array (", &
	766	SIZE(array, 2), ")."
	767	STOP
	768	END IF
	769
	770
	771	CASE( 'north scalar', 'south scalar', 'north w', 'south w' )
	772
	773	CALL check(nf90_put_var( ncid, var%varid, array(:,1,:), &
	774	start=start(1:ndim+1), &
	775	count=count(1:ndim) ) )
	776
	777
	778	CASE( 'surface forcing', 'top scalar', 'top w' )
	779
	780	CALL check(nf90_put_var( ncid, var%varid, array(:,:,1), &
	781	start=start(1:ndim+1), &
	782	count=count(1:ndim) ) )
	783
	784	CASE ( 'left u', 'right u', 'left v', 'right v' )
	785
	786	CALL check(nf90_put_var( ncid, var%varid, array(1,:,:), &
	787	start=start(1:ndim+1), &
	788	count=count(1:ndim) ) )
	789
	790	CASE ( 'north u', 'south u', 'north v', 'south v' )
	791
	792	CALL check(nf90_put_var( ncid, var%varid, array(:,1,:), &
	793	start=start(1:ndim+1), &
	794	count=count(1:ndim) ) )
	795
	796	CASE ( 'top u', 'top v' )
	797
	798	CALL check(nf90_put_var( ncid, var%varid, array(:,:,1), &
	799	start=start(1:ndim+1), &
	800	count=count(1:ndim) ) )
	801
	802	CASE ( 'time series' )
	803
	804	CALL check(nf90_put_var( ncid, var%varid, array(1,1,1), &
	805	start=start(1:ndim+1) ) )
	806
[3182]	807	CASE ( 'constant scalar profile' )
[2696]	808
	809	CALL check(nf90_put_var( ncid, var%varid, array(1,1,:), &
	810	start=start(1:ndim+1), &
	811	count=count(1:ndim) ) )
	812
[3182]	813	CASE ( 'large-scale scalar forcing', 'large-scale w forcing' )
	814
	815	message = "Doing nothing in terms of writing large-scale forings."
	816	CALL report('update_output', message)
	817
[2696]	818	CASE DEFAULT
	819
	820	message = "Variable kind '" // TRIM(var % kind) // &
	821	"' not recognized."
	822	CALL abort('update_output', message)
	823
	824	END SELECT
	825
	826	CALL check(nf90_close(ncid))
	827
	828	END SUBROUTINE update_output
	829
	830
	831	SUBROUTINE write_netcdf_variable_2d(var, iter, output_file, buffer)
	832	TYPE(nc_var), INTENT(IN) :: var
	833	INTEGER, INTENT(IN) :: iter
	834	TYPE(nc_file), INTENT(IN) :: output_file
	835	REAL(dp), INTENT(IN) :: buffer(:,:,:)
	836
	837	INTEGER :: ncid, ndim_out, start(4), count(4)
	838	LOGICAL :: last_dimension_is_time
	839
	840	ndim_out = var % ndim
	841
	842	last_dimension_is_time = var % dimids( var % ndim ) == output_file % dimid_time
	843	IF ( last_dimension_is_time ) THEN
	844	ndim_out = ndim_out - 1
	845	END IF
	846
	847	start(:) = (/1,1,1,iter/)
	848	count(1:ndim_out) = var%dimlen(1:ndim_out)
	849
	850	CALL check(nf90_open(output_file % name, NF90_WRITE, ncid))
	851
	852	IF (TRIM(var % kind) .EQ. 'left/right scalar') THEN
	853
	854	CALL check(nf90_put_var( ncid, var%varid, buffer(1,:,:) ) )
	855
	856	ELSE IF (TRIM(var % kind) .EQ. 'north/south scalar') THEN
	857
	858	CALL check(nf90_put_var( ncid, var%varid, buffer(:,1,:) ) )
	859
	860	ELSE IF (TRIM(var % kind) .EQ. 'top scalar') THEN
	861
	862	CALL check(nf90_put_var( ncid, var%varid, buffer(:,:,1) ) )
	863	ELSE
	864
	865	CALL check(nf90_put_var( ncid, var%varid, buffer ) )
	866
	867	END IF
	868	CALL check(nf90_close(ncid))
	869
	870	END SUBROUTINE write_netcdf_variable_2d
	871
	872
	873	SUBROUTINE check(status)
	874
	875	INTEGER, INTENT(IN) :: status
	876
	877	IF (status /= nf90_noerr) THEN
	878	message = "NetCDF API call failed with error: " // &
	879	TRIM( nf90_strerror(status) )
	880	CALL abort('io.check', message)
	881	END IF
	882
	883	END SUBROUTINE check
	884
	885	END MODULE io

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |