Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |