Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |