Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |