Home

Context Navigation

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

Last change on this file since 3507 was 3456, checked in by eckhard, 6 years ago
inifor: Removed surface forcing and internal arrays from netCDF output
Property svn:keywords set to `Id`
File size: 42.6 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |