Home

Context Navigation

io.f90 @ 3401

Last change on this file since 3401 was 3395, checked in by eckhard, 6 years ago
inifor: Added computation of geostrophic winds from COSMO input
Property svn:keywords set to `Id`
File size: 37.7 KB

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

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

Context Navigation

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

Download in other formats: