Home

Context Navigation

io.f90 @ 3402

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

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

Context Navigation

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

Download in other formats: