Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |