Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |