Home

Context Navigation

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

Last change on this file since 4798 was 4790, checked in by eckhard, 4 years ago
inifor: Validate netCDF dimensions of COSMO input files
Property svn:keywords set to `Id`
File size: 70.2 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-2020 Leibniz Universitaet Hannover
18	! Copyright 2017-2020 Deutscher Wetterdienst Offenbach
19	!------------------------------------------------------------------------------!
20	!
21	! Current revisions:
22	! -----------------
23	!
24	!
25	! Former revisions:
26	! -----------------
27	! $Id: inifor_io.f90 4790 2020-11-18 13:05:58Z schwenkel $
28	! Validate mesoscale input data by the number of dimensions and their lengths
29	!
30	!
31	! 4675 2020-09-11 10:00:26Z eckhard
32	! New command-line options for soil profile initialization and homogeneous
33	! boundary conditions
34	!
35	!
36	! 4659 2020-08-31 11:21:17Z eckhard
37	! New command-line option '--precipitation' enables output of precipitation
38	! surface forcing
39	! Only define netCDF variables in enabled IO groups
40	! Improved code formatting
41	!
42	!
43	! 4569 2020-06-19 12:46:17Z eckhard
44	! Removed unused variable
45	!
46	!
47	! 4568 2020-06-19 11:56:30Z eckhard
48	! Handle COSMO soil data with and without additional surface temperature
49	!
50	!
51	! 4553 2020-06-03 16:34:15Z eckhard
52	! Option --help now points the user to INIFOR's wiki page
53	! Remove deprecated command-line options -clon and -clat
54	!
55	!
56	! 4538 2020-05-18 13:45:35Z eckhard
57	! Make setting the vertical PALM origin mandatory
58	!
59	!
60	! 4523 2020-05-07 15:58:16Z eckhard
61	! respect integer working precision (iwp) specified in inifor_defs.f90
62	!
63	!
64	! 4481 2020-03-31 18:55:54Z maronga
65	! Pass hhl_file directly instead of entire INIFOR configuration
66	!
67	!
68	! 3997 2019-05-23 12:35:57Z eckhard
69	! Added boolean indicator for --elevation option invocation
70	! Stop INIFOR if no command-line options given
71	!
72	!
73	! 3866 2019-04-05 14:25:01Z eckhard
74	! Use PALM's working precision
75	! Improved coding style
76	!
77	!
78	! 3801 2019-03-15 17:14:25Z eckhard
79	! Added routine get_cosmo_grid() to read in COSMO rotated pole from COSMO domain
80	! Moved get_soil_layer_thickness() here from inifor_grid
81	!
82	! 3785 2019-03-06 10:41:14Z eckhard
83	! Temporariliy disabled height-based geostrophic wind averaging
84	! Improved variable naming
85	!
86	!
87	! 3764 2019-02-26 13:42:09Z eckhard
88	! Removed dependency on radiation input files
89	!
90	!
91	! 3716 2019-02-05 17:02:38Z eckhard
92	! Removed dependency on soilmoisture input files
93	!
94	!
95	! 3680 2019-01-18 14:54:12Z knoop
96	! Moved get_input_file_list() here from grid module, added check for presence of
97	! input files
98	!
99	!
100	! 3618 2018-12-10 13:25:22Z eckhard
101	! Prefixed all INIFOR modules with inifor_
102	!
103	!
104	! 3615 2018-12-10 07:21:03Z raasch
105	! bugfix: abort replaced by inifor_abort
106	!
107	! 3557 2018-11-22 16:01:22Z eckhard
108	! Updated documentation, removed unused subroutine write_netcdf_variable_2d()
109	!
110	!
111	! 3537 2018-11-20 10:53:14Z eckhard
112	! New routine get_netcdf_dim_vector()
113	!
114	!
115	! 3534 2018-11-19 15:35:16Z raasch
116	! bugfix: INTENT attribute changed
117	!
118	! 3456 2018-10-30 14:29:54Z eckhard
119	! NetCDf output of internal arrays only with --debug option
120	!
121	!
122	! 3447 2018-10-29 15:52:54Z eckhard
123	! Removed INCLUDE statement for get_netcdf_variable()
124	! Renamed source files for compatibilty with PALM build system
125	!
126	!
127	! 3395 2018-10-22 17:32:49Z eckhard
128	! Added command-line options for configuring the computation of geostrophic
129	! winds (--averagin-mode, --averaging-angle)
130	! Added command-line option --input-prefix for setting input file prefixes all
131	! at once
132	! Added --debug option for more verbose terminal output
133	! Added option-specific *_is_set LOGICALs to indicate invocation from the
134	! command-line
135	! Improved error messages in case of empty file-name strings
136	! Improved routine naming
137	!
138	! 3183 2018-07-27 14:25:55Z suehring
139	! Introduced new PALM grid stretching
140	! Updated variable names and metadata for PIDS v1.9 compatibility
141	! Improved handling of the start date string
142	! Better compatibility with older Intel compilers:
143	! - avoiding implicit array allocation with new get_netcdf_variable()
144	! subroutine instead of function
145	! Improved command line interface:
146	! - Added configuration validation
147	! - New options to configure input file prefixes
148	! - GNU-style short and long option names
149	! - Added version and copyright output
150	!
151	!
152	! 3182 2018-07-27 13:36:03Z suehring
153	! Initial revision
154	!
155	!
156	!
157	! Authors:
158	! --------
159	!> @author Eckhard Kadasch (Deutscher Wetterdienst, Offenbach)
160	!
161	! Description:
162	! ------------
163	!> The io module contains the functions needed to read and write netCDF data in
164	!> INIFOR.
165	!------------------------------------------------------------------------------!
166	#if defined ( __netcdf )
167	MODULE inifor_io
168
169	USE inifor_control
170	USE inifor_defs, &
171	ONLY: CFG_INIT_PROFILE, CFG_INIT_VOLUME, &
172	CFG_INIT_SOIL_PROFILE, CFG_INIT_SOIL_VOLUME, &
173	CFG_FORCING_HETERO, CFG_FORCING_HOMO, CFG_FORCING_NUDGING, &
174	DATE, SNAME, PATH, PI, TO_RADIANS, TO_DEGREES, VERSION, &
175	NC_DEPTH_DIM_IDX, NC_DEPTH_NAME, NC_HHL_NAME, NC_RLAT_NAME, &
176	NC_RLON_NAME, NC_ROTATED_POLE_NAME, NC_POLE_LATITUDE_NAME, &
177	NC_POLE_LONGITUDE_NAME, RHO_L, iwp, wp, &
178	PIDS_ORIGIN_LON, PIDS_ORIGIN_LAT, PIDS_ORIGIN_Z
179	USE inifor_types
180	USE inifor_util, &
181	ONLY: add_hours_to, nearly_equal, reverse, str, real_to_str
182	USE netcdf
183
184	IMPLICIT NONE
185
186	!------------------------------------------------------------------------------!
187	! Description:
188	! ------------
189	!> get_netcdf_variable() reads the netCDF data and metadate for the netCDF
190	!> variable 'in_var%name' from the file 'in_file'. The netCDF data array is
191	!> stored in the 'buffer' array and metadata added to the respective members of
192	!> the given 'in_var'.
193	!------------------------------------------------------------------------------!
194	INTERFACE get_netcdf_variable
195	MODULE PROCEDURE get_netcdf_variable_int
196	MODULE PROCEDURE get_netcdf_variable_real
197	END INTERFACE get_netcdf_variable
198
199	PRIVATE :: get_netcdf_variable_int, get_netcdf_variable_real
200
201	CONTAINS
202
203	!------------------------------------------------------------------------------!
204	! Description:
205	! ------------
206	!> get_netcdf_variable_int() implements the integer variant for the
207	!> get_netcdf_variable interface.
208	!------------------------------------------------------------------------------!
209	SUBROUTINE get_netcdf_variable_int(in_file, in_var, buffer)
210
211	CHARACTER(LEN=PATH), INTENT(IN) :: in_file
212	TYPE(nc_var), INTENT(INOUT) :: in_var
213	INTEGER(iwp), ALLOCATABLE, INTENT(INOUT) :: buffer(:,:,:)
214
215	INTEGER :: ncid
216	INTEGER, DIMENSION(3) :: start, count
217
218	IF ( nf90_open( TRIM(in_file), NF90_NOWRITE, ncid ) .EQ. NF90_NOERR .AND. &
219	nf90_inq_varid( ncid, in_var%name, in_var%varid ) .EQ. NF90_NOERR ) THEN
220
221	CALL get_input_dimensions(in_var, ncid)
222
223	CALL get_netcdf_start_and_count(in_var, start, count)
224	CALL log_runtime('time', 'read')
225
226	ALLOCATE( buffer( count(1), count(2), count(3) ) )
227	CALL log_runtime('time', 'alloc')
228
229	CALL check(nf90_get_var( ncid, in_var%varid, buffer, &
230	start = start, &
231	count = count ))
232
233	ELSE
234
235	message = "Failed to read '" // TRIM(in_var%name) // &
236	"' from file '" // TRIM(in_file) // "'."
237	CALL inifor_abort('get_netcdf_variable', message)
238
239	ENDIF
240
241	CALL check(nf90_close(ncid))
242	CALL log_runtime('time', 'read')
243
244	END SUBROUTINE get_netcdf_variable_int
245
246
247	!------------------------------------------------------------------------------!
248	! Description:
249	! ------------
250	!> get_netcdf_variable_real() implements the real variant for the
251	!> get_netcdf_variable interface.
252	!------------------------------------------------------------------------------!
253	SUBROUTINE get_netcdf_variable_real(in_file, in_var, buffer)
254
255	CHARACTER(LEN=PATH), INTENT(IN) :: in_file
256	TYPE(nc_var), INTENT(INOUT) :: in_var
257	REAL(wp), ALLOCATABLE, INTENT(INOUT) :: buffer(:,:,:)
258
259	INTEGER :: ncid
260	INTEGER, DIMENSION(3) :: start, count
261
262	IF ( nf90_open( TRIM(in_file), NF90_NOWRITE, ncid ) .EQ. NF90_NOERR .AND. &
263	nf90_inq_varid( ncid, in_var%name, in_var%varid ) .EQ. NF90_NOERR ) THEN
264
265	CALL get_input_dimensions(in_var, ncid)
266
267	CALL get_netcdf_start_and_count(in_var, start, count)
268	CALL log_runtime('time', 'read')
269
270	ALLOCATE( buffer( count(1), count(2), count(3) ) )
271	CALL log_runtime('time', 'alloc')
272
273	CALL check(nf90_get_var( ncid, in_var%varid, buffer, &
274	start = start, &
275	count = count ))
276
277	ELSE
278
279	message = "Failed to read '" // TRIM(in_var%name) // &
280	"' from file '" // TRIM(in_file) // "'."
281	CALL inifor_abort('get_netcdf_variable', message)
282
283	ENDIF
284
285	CALL check(nf90_close(ncid))
286	CALL log_runtime('time', 'read')
287
288	END SUBROUTINE get_netcdf_variable_real
289
290
291	!------------------------------------------------------------------------------!
292	! Description:
293	! ------------
294	!> get_netcdf_dim_vector() reads the coordinate array 'coordname' from the
295	!> netCDF file 'filename'.
296	!------------------------------------------------------------------------------!
297	SUBROUTINE get_netcdf_dim_vector(filename, coordname, coords)
298
299	CHARACTER(LEN=*), INTENT(IN) :: filename
300	CHARACTER(LEN=*), INTENT(IN) :: coordname
301	REAL(wp), ALLOCATABLE, INTENT(INOUT) :: coords(:)
302
303	INTEGER :: ncid, varid, dimlen
304	INTEGER :: dimids(NF90_MAX_VAR_DIMS)
305
306	IF ( nf90_open( TRIM(filename), NF90_NOWRITE, ncid ) .EQ. NF90_NOERR .AND. &
307	nf90_inq_varid( ncid, coordname, varid ) .EQ. NF90_NOERR ) THEN
308
309	CALL check(nf90_inquire_variable( ncid, varid, dimids = dimids ))
310	CALL check(nf90_inquire_dimension( ncid, dimids(1), len = dimlen ))
311
312	ALLOCATE(coords(dimlen))
313	CALL check(nf90_get_var( ncid, varid, coords))
314
315	ELSE
316
317	message = "Failed to read '" // TRIM(coordname) // &
318	"' from file '" // TRIM(filename) // "'."
319	CALL inifor_abort('get_netcdf_dim_vector', message)
320
321	ENDIF
322
323	END SUBROUTINE get_netcdf_dim_vector
324
325
326	!------------------------------------------------------------------------------!
327	! Description:
328	! ------------
329	!> get_input_dimensions() reads dimensions metadata of the netCDF variable given
330	!> by 'in_var%name' into 'in_var' data structure.
331	!------------------------------------------------------------------------------!
332	SUBROUTINE get_input_dimensions(in_var, ncid)
333
334	TYPE(nc_var), INTENT(INOUT) :: in_var
335	INTEGER, INTENT(IN) :: ncid
336
337	INTEGER :: i
338
339	CALL check(nf90_get_att( ncid, in_var%varid, "long_name", &
340	in_var%long_name))
341
342	CALL check(nf90_get_att( ncid, in_var%varid, "units", &
343	in_var%units ))
344
345	CALL check(nf90_inquire_variable( ncid, in_var%varid, &
346	ndims = in_var%ndim, &
347	dimids = in_var%dimids ))
348
349	DO i = 1, in_var%ndim
350	CALL check(nf90_inquire_dimension( ncid, in_var%dimids(i), &
351	name = in_var%dimname(i), &
352	len = in_var%dimlen(i) ))
353	ENDDO
354
355	END SUBROUTINE get_input_dimensions
356
357
358	!------------------------------------------------------------------------------!
359	! Description:
360	! ------------
361	!> get_netcdf_start_and_count() gets the start position and element counts for
362	!> the given netCDF file. This information is used in get_netcdf_variable_int()
363	!> and _real() for reading input variables..
364	!------------------------------------------------------------------------------!
365	SUBROUTINE get_netcdf_start_and_count(in_var, start, count)
366
367	TYPE(nc_var), INTENT(INOUT) :: in_var
368	INTEGER, DIMENSION(3), INTENT(OUT) :: start, count
369
370	INTEGER :: ndim
371
372	IF ( in_var%ndim .LT. 2 .OR. in_var%ndim .GT. 4 ) THEN
373
374	message = "Failed reading NetCDF variable " // &
375	TRIM(in_var%name) // " with " // &
376	TRIM(str(INT(in_var%ndim, kind=iwp))) // &
377	" dimensions because only two- and and three-dimensional" // &
378	" variables are supported."
379	CALL inifor_abort('get_netcdf_start_and_count', message)
380
381	ENDIF
382
383	start = (/ 1, 1, 1 /)
384	IF ( TRIM(in_var%name) .EQ. 'T_SO' .AND. &
385	in_var%has_redundant_first_level ) THEN
386
387	!
388	!-- Skip depth = 0.0 for T_SO and reduce number of depths from 9 to 8
389	in_var%dimlen(3) = in_var%dimlen(3) - 1
390
391	!
392	!-- Start reading from second level, e.g. depth = 0.005 instead of 0.0
393	start(3) = 2
394	ENDIF
395
396	IF (in_var%ndim .EQ. 2) THEN
397	in_var%dimlen(3) = 1
398	ENDIF
399
400	ndim = MIN(in_var%ndim, 3)
401	count = (/ 1, 1, 1 /)
402	count(1:ndim) = in_var%dimlen(1:ndim)
403
404	END SUBROUTINE get_netcdf_start_and_count
405
406
407	!------------------------------------------------------------------------------!
408	! Description:
409	! ------------
410	!> Routine for defining netCDF variables in the dynamic driver, INIFOR's netCDF
411	!> output.
412	!------------------------------------------------------------------------------!
413	SUBROUTINE netcdf_define_variable(var, ncid)
414
415	TYPE(nc_var), INTENT(INOUT) :: var
416	INTEGER, INTENT(IN) :: ncid
417
418	CALL check(nf90_def_var(ncid, var%name, NF90_FLOAT, var%dimids(1:var%ndim), var%varid))
419	CALL check(nf90_put_att(ncid, var%varid, "long_name", var%long_name))
420	CALL check(nf90_put_att(ncid, var%varid, "units", var%units))
421	IF ( var%lod .GE. 0 ) THEN
422	CALL check(nf90_put_att(ncid, var%varid, "lod", var%lod))
423	ENDIF
424	CALL check(nf90_put_att(ncid, var%varid, "source", var%source))
425	CALL check(nf90_put_att(ncid, var%varid, "_FillValue", NF90_FILL_REAL))
426
427	END SUBROUTINE netcdf_define_variable
428
429
430	!------------------------------------------------------------------------------!
431	! Description:
432	! ------------
433	!> netcdf_get_dimensions() reads in all dimensions and their lengths and stores
434	!> them in the given the 'var' data structure. This information is used later
435	!> for writing output variables in update_output().
436	!------------------------------------------------------------------------------!
437	SUBROUTINE netcdf_get_dimensions(var, ncid)
438
439	TYPE(nc_var), INTENT(INOUT) :: var
440	INTEGER, INTENT(IN) :: ncid
441	INTEGER :: i
442	CHARACTER(SNAME) :: null
443
444	DO i = 1, var%ndim
445	CALL check(nf90_inquire_dimension(ncid, var%dimids(i), &
446	name = null, &
447	len = var%dimlen(i) ) )
448	ENDDO
449
450	END SUBROUTINE netcdf_get_dimensions
451
452
453	!------------------------------------------------------------------------------!
454	! Description:
455	! ------------
456	!> This routine parses and interpretes the command-line options and stores the
457	!> resulting settings in the 'cfg' data structure.
458	!------------------------------------------------------------------------------!
459	SUBROUTINE parse_command_line_arguments( cfg )
460
461	TYPE(inifor_config), INTENT(INOUT) :: cfg
462
463	CHARACTER(LEN=PATH) :: option, arg
464	INTEGER :: arg_count, i
465
466	cfg%flow_prefix_is_set = .FALSE.
467	cfg%input_prefix_is_set = .FALSE.
468	cfg%p0_is_set = .FALSE.
469	cfg%radiation_prefix_is_set = .FALSE.
470	cfg%soil_prefix_is_set = .FALSE.
471	cfg%precipitation_prefix_is_set = .FALSE.
472	cfg%process_precipitation = .FALSE.
473	cfg%static_driver_is_set = .FALSE.
474	cfg%ug_defined_by_user = .FALSE.
475	cfg%vg_defined_by_user = .FALSE.
476	cfg%z0_is_set = .FALSE.
477
478	arg_count = COMMAND_ARGUMENT_COUNT()
479	IF (arg_count .GT. 0) THEN
480
481	i = 1
482	DO WHILE (i .LE. arg_count)
483
484	CALL GET_COMMAND_ARGUMENT( i, option )
485
486	SELECT CASE( TRIM(option) )
487
488	CASE( '--averaging-mode' )
489	CALL get_option_argument( i, arg )
490	cfg%averaging_mode = TRIM(arg)
491
492	CASE( '-date', '-d', '--date' )
493	CALL get_option_argument( i, arg )
494	cfg%start_date = TRIM(arg)
495
496	CASE( '--debug' )
497	cfg%debug = .TRUE.
498
499	CASE( '-z0', '-z', '--elevation' )
500	cfg%z0_is_set = .TRUE.
501	CALL get_option_argument( i, arg )
502	READ(arg, *) cfg%z0
503
504	CASE( '-p0', '-r', '--surface-pressure' )
505	cfg%p0_is_set = .TRUE.
506	CALL get_option_argument( i, arg )
507	READ(arg, *) cfg%p0
508
509	CASE( '-ug', '-u', '--geostrophic-u' )
510	cfg%ug_defined_by_user = .TRUE.
511	CALL get_option_argument( i, arg )
512	READ(arg, *) cfg%ug
513
514	CASE( '-vg', '-v', '--geostrophic-v' )
515	cfg%vg_defined_by_user = .TRUE.
516	CALL get_option_argument( i, arg )
517	READ(arg, *) cfg%vg
518
519	CASE( '-path', '-p', '--path' )
520	CALL get_option_argument( i, arg )
521	cfg%input_path = TRIM(arg)
522
523	CASE( '-hhl', '-l', '--hhl-file' )
524	CALL get_option_argument( i, arg )
525	cfg%hhl_file = TRIM(arg)
526
527	CASE( '--input-prefix')
528	CALL get_option_argument( i, arg )
529	cfg%input_prefix = TRIM(arg)
530	cfg%input_prefix_is_set = .TRUE.
531
532	CASE( '-a', '--averaging-angle' )
533	CALL get_option_argument( i, arg )
534	READ(arg, *) cfg%averaging_angle
535
536	CASE( '-static', '-t', '--static-driver' )
537	cfg%static_driver_is_set = .TRUE.
538	CALL get_option_argument( i, arg )
539	cfg%static_driver_file = TRIM(arg)
540
541	CASE( '-soil', '-s', '--soil-file')
542	CALL get_option_argument( i, arg )
543	cfg%soiltyp_file = TRIM(arg)
544
545	CASE( '--flow-prefix')
546	CALL get_option_argument( i, arg )
547	cfg%flow_prefix = TRIM(arg)
548	cfg%flow_prefix_is_set = .TRUE.
549
550	CASE( '--radiation-prefix')
551	CALL get_option_argument( i, arg )
552	cfg%radiation_prefix = TRIM(arg)
553	cfg%radiation_prefix_is_set = .TRUE.
554
555	CASE( '--soil-prefix')
556	CALL get_option_argument( i, arg )
557	cfg%soil_prefix = TRIM(arg)
558	cfg%soil_prefix_is_set = .TRUE.
559
560	CASE( '--precipitation-prefix')
561	CALL get_option_argument( i, arg )
562	cfg%precipitation_prefix = TRIM(arg)
563	cfg%precipitation_prefix_is_set = .TRUE.
564
565	CASE( '--precipitation')
566	cfg%process_precipitation = .TRUE.
567
568	CASE( '-o', '--output' )
569	CALL get_option_argument( i, arg )
570	cfg%output_file = TRIM(arg)
571
572	CASE( '-n', '--namelist' )
573	CALL get_option_argument( i, arg )
574	cfg%namelist_file = TRIM(arg)
575
576	CASE( '-i', '--init-mode' )
577	CALL get_option_argument( i, arg )
578	cfg%ic_mode = TRIM(arg)
579
580	CASE( '--soil-init-mode' )
581	CALL get_option_argument( i, arg )
582	cfg%isc_mode = TRIM(arg)
583
584	CASE( '-f', '--forcing-mode' )
585	CALL get_option_argument( i, arg )
586	cfg%bc_mode = TRIM(arg)
587
588	CASE( '--version' )
589	CALL print_version
590	STOP
591
592	CASE( '--help' )
593	CALL print_version
594	PRINT *, ""
595	PRINT *, &
596	"For documentation and a list of available command-line options " // NEW_LINE( " " ) // &
597	" please visit https://palm.muk.uni-hannover.de/trac/wiki/doc/app/iofiles/inifor."
598	STOP
599
600	CASE DEFAULT
601	message = "unknown option '" // TRIM(option) // "'."
602	CALL inifor_abort('parse_command_line_arguments', message)
603
604	END SELECT
605
606	i = i + 1
607
608	ENDDO
609
610	ELSE
611
612	CALL print_version
613	CALL report( 'parse_command_line_arguments', 'No arguments present, exiting.' )
614	STOP
615
616	ENDIF
617
618	END SUBROUTINE parse_command_line_arguments
619
620
621
622	SUBROUTINE get_datetime_file_list( start_date_string, start_hour, end_hour, &
623	step_hour, input_path, prefix, suffix, &
624	file_list )
625
626	CHARACTER (LEN=DATE), INTENT(IN) :: start_date_string
627	CHARACTER (LEN=*), INTENT(IN) :: prefix, suffix, input_path
628	INTEGER(iwp), INTENT(IN) :: start_hour, end_hour, step_hour
629	CHARACTER(LEN=*), ALLOCATABLE, INTENT(INOUT) :: file_list(:)
630
631	INTEGER(iwp) :: number_of_intervals, hour, i
632	CHARACTER(LEN=DATE) :: date_string
633
634	number_of_intervals = CEILING( REAL(end_hour - start_hour) / step_hour )
635	ALLOCATE( file_list(number_of_intervals + 1) )
636
637	DO i = 0, number_of_intervals
638
639	hour = start_hour + i * step_hour
640	date_string = add_hours_to(start_date_string, hour)
641
642	file_list(i+1) = TRIM(input_path) // TRIM(prefix) // &
643	TRIM(date_string) // TRIM(suffix) // '.nc'
644
645	ENDDO
646
647	END SUBROUTINE get_datetime_file_list
648
649	!------------------------------------------------------------------------------!
650	! Description:
651	! ------------
652	!> Establish a list of files based on the given start and end times and file
653	!> prefixes and suffixes.
654	!------------------------------------------------------------------------------!
655	SUBROUTINE get_input_file_list( start_date_string, start_hour, end_hour, &
656	step_hour, input_path, prefix, suffix, &
657	file_list )
658
659	CHARACTER (LEN=DATE), INTENT(IN) :: start_date_string
660	CHARACTER (LEN=*), INTENT(IN) :: prefix, suffix, input_path
661	INTEGER(iwp), INTENT(IN) :: start_hour, end_hour, step_hour
662	CHARACTER(LEN=*), ALLOCATABLE, INTENT(INOUT) :: file_list(:)
663
664	CALL get_datetime_file_list( start_date_string, start_hour, end_hour, &
665	step_hour, input_path, prefix, suffix, &
666	file_list )
667
668	END SUBROUTINE get_input_file_list
669
670
671	!------------------------------------------------------------------------------!
672	! Description:
673	! ------------
674	!> Abort INIFOR if the given file is not present.
675	!------------------------------------------------------------------------------!
676	LOGICAL FUNCTION file_is_present( filename, file_kind, message )
677
678	CHARACTER(LEN=*), INTENT(IN) :: filename, file_kind
679	CHARACTER(LEN=*), INTENT(OUT) :: message
680
681	file_is_present = file_present(filename)
682
683	IF (.NOT. file_is_present) THEN
684
685	IF (LEN( TRIM( filename ) ) == 0) THEN
686	message = "No name was given for the " // TRIM( file_kind ) // " file."
687	ELSE
688	message = "The " // TRIM( file_kind ) // " file '" // &
689	TRIM( filename ) // "' was not found."
690	ENDIF
691
692	ENDIF
693
694	END FUNCTION file_is_present
695
696
697	!------------------------------------------------------------------------------!
698	! Description:
699	! ------------
700	!> Abort INIFOR if the given file is not present.
701	!------------------------------------------------------------------------------!
702	SUBROUTINE verify_file( file_name, file_kind )
703
704	CHARACTER(LEN=*), INTENT(IN) :: file_name, file_kind
705
706	IF (.NOT. file_is_present( file_name, file_kind, message )) THEN
707
708	CALL inifor_abort( 'verify_file', message )
709
710	ENDIF
711
712	message = "Set up input file name '" // TRIM( file_name ) // "'"
713	CALL report( 'verify_file', message )
714
715	END SUBROUTINE verify_file
716	!------------------------------------------------------------------------------!
717	! Description:
718	! ------------
719	!> Get the argument of the i'th command line option, which is at the location
720	!> i+1 of the argument list.
721	!------------------------------------------------------------------------------!
722	SUBROUTINE get_option_argument( i, arg )
723	CHARACTER(LEN=PATH), INTENT(INOUT) :: arg
724	INTEGER, INTENT(INOUT) :: i
725
726	i = i + 1
727	CALL GET_COMMAND_ARGUMENT( i, arg )
728
729	END SUBROUTINE
730
731
732	!------------------------------------------------------------------------------!
733	! Description:
734	! ------------
735	!> Checks the INIFOR configuration 'cfg' for plausibility.
736	!------------------------------------------------------------------------------!
737	SUBROUTINE validate_config( cfg )
738	TYPE(inifor_config), INTENT(IN) :: cfg
739
740	CALL verify_file( cfg%hhl_file, 'HHL' )
741	CALL verify_file( cfg%namelist_file, 'NAMELIST' )
742	CALL verify_file( cfg%soiltyp_file, 'SOILTYP' )
743
744	!
745	!-- Only check optional static driver file name, if it has been given.
746	IF (TRIM( cfg%static_driver_file ) .NE. '') THEN
747	CALL verify_file( cfg%static_driver_file, 'static driver' )
748	ENDIF
749
750	SELECT CASE( TRIM( cfg%ic_mode ) )
751	CASE( CFG_INIT_PROFILE, CFG_INIT_VOLUME )
752	CASE DEFAULT
753	message = "Initialization mode '" // TRIM( cfg%ic_mode ) //&
754	"' is not supported. " //&
755	"Please select either '" // CFG_INIT_PROFILE //"' or '" // &
756	CFG_INIT_VOLUME //"', " //&
757	"or omit the -i/--init-mode option entirely, which corresponds "//&
758	"to the latter."
759	CALL inifor_abort( 'validate_config', message )
760	END SELECT
761
762	SELECT CASE( TRIM( cfg%isc_mode ) )
763	CASE( CFG_INIT_SOIL_PROFILE, CFG_INIT_SOIL_VOLUME )
764	CASE DEFAULT
765	message = "Soil initialization mode '" // TRIM( cfg%isc_mode ) //&
766	"' is not supported. " //&
767	"Please select either '" // CFG_INIT_SOIL_PROFILE //"' or '" // &
768	CFG_INIT_SOIL_VOLUME //"', " //&
769	"or omit the --soil-init-mode option entirely, which corresponds "//&
770	"to the latter."
771	CALL inifor_abort( 'validate_config', message )
772	END SELECT
773
774	SELECT CASE( TRIM(cfg%bc_mode) )
775	CASE( CFG_FORCING_HOMO, CFG_FORCING_HETERO, CFG_FORCING_NUDGING )
776	CASE DEFAULT
777	message = "Forcing mode '" // TRIM( cfg%bc_mode ) //&
778	"' is not supported. " //&
779	"Please select either '" // CFG_FORCING_NUDGING // &
780	"', '" // CFG_FORCING_HOMO // "', or '" // &
781	CFG_FORCING_HETERO // "' " //&
782	"or omit the -f/--forcing-mode option entirely, which corresponds "//&
783	"to the latter."
784	CALL inifor_abort( 'validate_config', message )
785	END SELECT
786
787	SELECT CASE( TRIM( cfg%averaging_mode ) )
788	CASE( 'level' )
789	CASE( 'height' )
790	message = "Averaging mode '" // TRIM( cfg%averaging_mode ) //&
791	"' is currently not supported. " //&
792	"Please use level-based averaging by selecting 'level', " //&
793	"or by omitting the --averaging-mode option entirely."
794	CALL inifor_abort( 'validate_config', message )
795	CASE DEFAULT
796	message = "Averaging mode '" // TRIM( cfg%averaging_mode ) //&
797	"' is not supported. " //&
798	! "Please select either 'height' or 'level', " //&
799	! "or omit the --averaging-mode option entirely, which corresponds "//&
800	! "to the latter."
801	"Please use level-based averaging by selecting 'level', " //&
802	"or by omitting the --averaging-mode option entirely."
803	CALL inifor_abort( 'validate_config', message )
804	END SELECT
805
806	IF ( cfg%ug_defined_by_user .NEQV. cfg%vg_defined_by_user ) THEN
807	message = "You specified only one component of the geostrophic " // &
808	"wind. Please specify either both or none."
809	CALL inifor_abort( 'validate_config', message )
810	ENDIF
811
812	IF ( .NOT. cfg%static_driver_is_set .AND. .NOT. cfg%z0_is_set ) THEN
813	message = &
814	"The vertical origin of the PALM grid has not been defined. " // NEW_LINE( " " ) // &
815	"Please specify the right value for your setup by either " // NEW_LINE( " " ) // &
816	" - using the command-line option --elevation <height above sea level>, or by" // NEW_LINE( " " ) // &
817	" - specifying a static driver file using --static <filename> in order to use " // NEW_LINE( " " ) // &
818	" use the value of origin_z (and origin_lon and origin_lat) specifed therein."
819	CALL inifor_abort( 'validate_config', message )
820	ENDIF
821
822	END SUBROUTINE validate_config
823
824
825	!------------------------------------------------------------------------------!
826	! Description:
827	! ------------
828	!> Checks wheather the COSMO grid matches the shape of the meteorological input
829	!> data by comparing the number of netCDF dimensions and their lengths in the
830	!> hhl.nc and the first of the *-flow files.
831	!------------------------------------------------------------------------------!
832	SUBROUTINE validate_dataset(flow_files, hhl_file)
833	CHARACTER(LEN=PATH), INTENT(IN) :: flow_files(:) !< paths to files containing atmospheric variables
834	CHARACTER(LEN=PATH), INTENT(IN) :: hhl_file !< path to file containing the HHL variable (height of half layers)
835
836	CHARACTER(SNAME), PARAMETER :: NC_W_NAME = 'W'
837	TYPE(nc_var) :: hhl_var, flow_var
838	INTEGER :: dim_idx, ncid, ndims_flow, ndims_hhl, varid
839	REAL(wp), ALLOCATABLE :: hhl_dim_vector(:), flow_dim_vector(:)
840	LOGICAL :: dims_have_same_length
841
842	IF ( nf90_open( TRIM( flow_files(1) ), NF90_NOWRITE, ncid ) .EQ. NF90_NOERR ) THEN
843
844	CALL check( nf90_inq_varid( ncid, NC_W_NAME, varid=varid ) )
845	CALL check( nf90_inquire_variable( ncid, varid, ndims=ndims_flow ) )
846	CALL netcdf_get_dimensions( flow_var, ncid )
847	CALL check( nf90_close( ncid ) )
848
849	ELSE
850
851	message = "Failed to read netCDF dimensions'" // &
852	"' from file '" // TRIM( flow_files(1) ) // "'."
853	CALL inifor_abort( 'validate_dataset', message )
854
855	ENDIF
856
857	IF ( nf90_open( TRIM( hhl_file ), NF90_NOWRITE, ncid ) .EQ. NF90_NOERR ) THEN
858
859	CALL check( nf90_inq_varid( ncid, NC_HHL_NAME, varid=varid ) )
860	CALL check( nf90_inquire_variable( ncid, varid, ndims=ndims_hhl ) )
861	CALL netcdf_get_dimensions( hhl_var, ncid )
862	CALL check( nf90_close( ncid ) )
863
864	ELSE
865
866	message = "Failed to read netCDF dimensions'" // &
867	"' from file '" // TRIM(hhl_file) // "'."
868	CALL inifor_abort( 'validate_dataset', message )
869
870	ENDIF
871
872	!
873	!-- Compare number dimensions of 'HHL' in hhl file and 'W' in first flow file
874	IF ( .NOT. ndims_flow .EQ. ndims_hhl ) THEN
875	message = "Mesoscale data inconsistent. Number of dimensions in the " //&
876	"hhl file does not match with the meteorologial fields " // &
877	"in the *-flow files (" // &
878	"HHL: ndims = " // TRIM( str( ndims_hhl ) ) // ", " // &
879	"W: ndims = " // TRIM( str( ndims_flow ) ) // ")."
880	CALL inifor_abort( 'validate_dataset', message )
881	ENDIF
882
883
884	!
885	!-- Compare lengths of each dimension, ignoring time (dim_idx = 1)
886	DO dim_idx = 2, ndims_hhl
887
888	CALL get_dimension_vector_of_variable( &
889	NC_HHL_NAME, &
890	dim_idx = dim_idx, &
891	filename = hhl_file, &
892	dim_vector = hhl_dim_vector &
893	)
894
895	CALL get_dimension_vector_of_variable( &
896	NC_W_NAME, &
897	dim_idx = dim_idx, &
898	filename = flow_files(1), &
899	dim_vector = flow_dim_vector &
900	)
901
902	dims_have_same_length = SIZE( flow_dim_vector ) .EQ. SIZE( hhl_dim_vector )
903	IF ( .NOT. dims_have_same_length ) THEN
904	message = &
905	"Mesoscale data inconsistent. Number of grid points " // &
906	"in dimension #" // TRIM( str( dim_idx ) ) // &
907	" do not match in the hhl and *-flow files (" // &
908	"HHL: len = " // TRIM( str( SIZE( hhl_dim_vector ) ) ) // ", " // &
909	"W: len = " // TRIM( str( SIZE( flow_dim_vector ) ) )// ")."
910	CALL inifor_abort( 'validate_dataset', message )
911	ENDIF
912
913	ENDDO
914
915	END SUBROUTINE validate_dataset
916
917	SUBROUTINE get_cosmo_grid( hhl_file, soil_file, rlon, rlat, hhl, hfl, depths, &
918	d_depth, d_depth_rho_inv, phi_n, lambda_n, &
919	phi_equat, &
920	lonmin_cosmo, lonmax_cosmo, &
921	latmin_cosmo, latmax_cosmo, &
922	nlon, nlat, nlev, ndepths )
923
924	CHARACTER(LEN=PATH), INTENT(IN) :: hhl_file !< path to file containing the HHL variable (height of half layers)
925	CHARACTER(LEN=PATH), INTENT(IN) :: soil_file !< path to one of the soil input files for reading soil layer depths
926	REAL(wp), DIMENSION(:), ALLOCATABLE, INTENT(OUT) :: rlon !< longitudes of COSMO-DE's rotated-pole grid
927	REAL(wp), DIMENSION(:), ALLOCATABLE, INTENT(OUT) :: rlat !< latitudes of COSMO-DE's rotated-pole grid
928	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, INTENT(OUT) :: hhl !< heights of half layers (cell faces) above sea level in COSMO-DE, read in from external file
929	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, INTENT(OUT) :: hfl !< heights of full layers (cell centres) above sea level in COSMO-DE, computed as arithmetic average of hhl
930	REAL(wp), DIMENSION(:), ALLOCATABLE, INTENT(OUT) :: depths !< COSMO-DE's TERRA-ML soil layer depths
931	REAL(wp), DIMENSION(:), ALLOCATABLE, INTENT(OUT) :: d_depth
932	REAL(wp), DIMENSION(:), ALLOCATABLE, INTENT(OUT) :: d_depth_rho_inv
933	REAL(wp), INTENT(OUT) :: phi_n
934	REAL(wp), INTENT(OUT) :: phi_equat
935	REAL(wp), INTENT(OUT) :: lambda_n
936	REAL(wp), INTENT(OUT) :: lonmin_cosmo !< Minimunm longitude of COSMO-DE's rotated-pole grid [COSMO rotated-pole rad]
937	REAL(wp), INTENT(OUT) :: lonmax_cosmo !< Maximum longitude of COSMO-DE's rotated-pole grid [COSMO rotated-pole rad]
938	REAL(wp), INTENT(OUT) :: latmin_cosmo !< Minimunm latitude of COSMO-DE's rotated-pole grid [COSMO rotated-pole rad]
939	REAL(wp), INTENT(OUT) :: latmax_cosmo !< Maximum latitude of COSMO-DE's rotated-pole grid [COSMO rotated-pole rad]
940	INTEGER(iwp), INTENT(OUT) :: nlon, nlat, nlev, ndepths
941
942	TYPE(nc_var) :: cosmo_var !< COSMO dummy variable, used for reading HHL, rlon, rlat
943	INTEGER(iwp) :: k
944
945	!
946	!-- Read in COSMO's heights of half layers (vertical cell faces)
947	cosmo_var%name = NC_HHL_NAME
948	CALL get_netcdf_variable( hhl_file, cosmo_var, hhl )
949	CALL get_netcdf_dim_vector( hhl_file, NC_RLON_NAME, rlon )
950	CALL get_netcdf_dim_vector( hhl_file, NC_RLAT_NAME, rlat )
951	CALL get_netcdf_dim_vector( soil_file, NC_DEPTH_NAME, depths)
952	CALL log_runtime( 'time', 'read' )
953
954	CALL reverse( hhl )
955	nlon = SIZE( hhl, 1 )
956	nlat = SIZE( hhl, 2 )
957	nlev = SIZE( hhl, 3 )
958	ndepths = SIZE( depths )
959
960	CALL log_runtime( 'time', 'comp' )
961
962	ALLOCATE( hfl( nlon, nlat, nlev-1 ) )
963	ALLOCATE( d_depth( ndepths ), d_depth_rho_inv( ndepths ) )
964	CALL log_runtime('time', 'alloc')
965
966	CALL get_soil_layer_thickness( depths, d_depth )
967	d_depth_rho_inv = 1.0_wp / ( d_depth * RHO_L )
968
969	!
970	!-- Compute COSMO's heights of full layers (cell centres)
971	DO k = 1, nlev-1
972	hfl(:,:,k) = 0.5_wp * ( hhl(:,:,k) + &
973	hhl(:,:,k+1) )
974	ENDDO
975	!
976	!-- COSMO rotated pole coordinates
977	phi_n = TO_RADIANS &
978	* get_netcdf_variable_attribute( hhl_file, &
979	NC_ROTATED_POLE_NAME, &
980	NC_POLE_LATITUDE_NAME )
981
982	lambda_n = TO_RADIANS &
983	* get_netcdf_variable_attribute( hhl_file, &
984	NC_ROTATED_POLE_NAME, &
985	NC_POLE_LONGITUDE_NAME )
986
987	phi_equat = 90.0_wp * TO_RADIANS - phi_n
988
989	lonmin_cosmo = MINVAL( rlon ) * TO_RADIANS
990	lonmax_cosmo = MAXVAL( rlon ) * TO_RADIANS
991	latmin_cosmo = MINVAL( rlat ) * TO_RADIANS
992	latmax_cosmo = MAXVAL( rlat ) * TO_RADIANS
993	CALL log_runtime('time', 'comp')
994
995	END SUBROUTINE get_cosmo_grid
996
997
998	!------------------------------------------------------------------------------!
999	! Description:
1000	! ------------
1001	!> Fills the thickness array of the COSMO soil layers. Since COSMO's (i.e.
1002	!> TERRA_ML's [1]) soil layer boundaries follow the rule
1003	!>
1004	!> depth(0) = 0.0, and
1005	!> depth(k) = 0.01 * 3**(k-1), k in [1,2,3,...,7]
1006	!>
1007	!> and full levels are defined as the midpoints between two layer boundaries,
1008	!> all except the first layer thicknesses equal the depth of the midpoint.
1009	!>
1010	!> [1] A Description of the Nonhydrostatic Regional COSMO Model Part II :
1011	!> Physical Parameterization*, Sect. 11 TERRA_ML.
1012	!> http://www.cosmo-model.org/content/model/documentation/core/cosmoPhysParamtr.pdf)
1013	!>
1014	!> Input parameters:
1015	!> -----------------
1016	!>
1017	!> depths: array of full soil layer depths (cell midpoints)
1018	!>
1019	!>
1020	!> Output parameters:
1021	!> ------------------
1022	!>
1023	!> d_depth: array of soil layer thicknesses
1024	!>
1025	!------------------------------------------------------------------------------!
1026	SUBROUTINE get_soil_layer_thickness( depths, d_depth )
1027
1028	REAL(wp), INTENT(IN) :: depths(:)
1029	REAL(wp), INTENT(OUT) :: d_depth(:)
1030
1031	d_depth(:) = depths(:)
1032	d_depth(1) = 2.0_wp * depths(1)
1033
1034	END SUBROUTINE get_soil_layer_thickness
1035	!------------------------------------------------------------------------------!
1036	! Description:
1037	! ------------
1038	!> Check whether the given file is present on the filesystem.
1039	!------------------------------------------------------------------------------!
1040	LOGICAL FUNCTION file_present( filename )
1041	CHARACTER(LEN=PATH), INTENT(IN) :: filename
1042
1043	INQUIRE( FILE=filename, EXIST=file_present )
1044
1045	END FUNCTION file_present
1046
1047
1048	!------------------------------------------------------------------------------!
1049	! Description:
1050	! ------------
1051	!> This routine initializes the dynamic driver file, i.e. INIFOR's netCDF output
1052	!> file.
1053	!>
1054	!> Besides writing metadata, such as global attributes, coordinates, variables,
1055	!> in the netCDF file, various netCDF IDs are saved for later, when INIFOR
1056	!> writes the actual data.
1057	!------------------------------------------------------------------------------!
1058	SUBROUTINE setup_netcdf_dimensions( output_file, palm_grid, &
1059	start_date_string, origin_lon, origin_lat )
1060
1061	TYPE(nc_file), INTENT(INOUT) :: output_file
1062	TYPE(grid_definition), INTENT(IN) :: palm_grid
1063	CHARACTER (LEN=DATE), INTENT(IN) :: start_date_string
1064	REAL(wp), INTENT(IN) :: origin_lon, origin_lat
1065
1066	CHARACTER (LEN=8) :: date_string
1067	CHARACTER (LEN=10) :: time_string
1068	CHARACTER (LEN=5) :: zone_string
1069	CHARACTER (LEN=SNAME) :: history_string
1070	INTEGER :: nx, ny, nz, nt
1071	INTEGER :: ncid, dimids(3), dimvarids(3)
1072	REAL(wp) :: z0
1073
1074	message = "Initializing PALM-4U dynamic driver file '" // &
1075	TRIM(output_file%name) // "' and setting up dimensions."
1076	CALL report('setup_netcdf_dimensions', message)
1077
1078	!
1079	!-- Create the netCDF file as in netCDF-4/HDF5 format if __netcdf4 preprocessor flag is given
1080	#if defined( __netcdf4 )
1081	CALL check(nf90_create(TRIM(output_file%name), OR(NF90_CLOBBER, NF90_HDF5), ncid))
1082	#else
1083	CALL check(nf90_create(TRIM(output_file%name), NF90_CLOBBER, ncid))
1084	#endif
1085
1086	!------------------------------------------------------------------------------
1087	!- Section 1: Define NetCDF dimensions and coordinates
1088	!------------------------------------------------------------------------------
1089	nt = SIZE(output_file%time)
1090	nx = palm_grid%nx
1091	ny = palm_grid%ny
1092	nz = palm_grid%nz
1093	z0 = palm_grid%z0
1094
1095
1096	!
1097	!------------------------------------------------------------------------------
1098	!- Section 2: Write global NetCDF attributes
1099	!------------------------------------------------------------------------------
1100	CALL date_and_time(DATE=date_string, TIME=time_string, ZONE=zone_string)
1101	history_string = &
1102	'Created on '// date_string // &
1103	' at ' // time_string(1:2) // ':' // time_string(3:4) // &
1104	' (UTC' // zone_string // ')'
1105
1106	CALL check(nf90_put_att(ncid, NF90_GLOBAL, 'title', 'PALM input file for scenario ...'))
1107	CALL check(nf90_put_att(ncid, NF90_GLOBAL, 'institution', 'Deutscher Wetterdienst, Offenbach'))
1108	CALL check(nf90_put_att(ncid, NF90_GLOBAL, 'author', 'Eckhard Kadasch, eckhard.kadasch@dwd.de'))
1109	CALL check(nf90_put_att(ncid, NF90_GLOBAL, 'history', TRIM(history_string)))
1110	CALL check(nf90_put_att(ncid, NF90_GLOBAL, 'references', '--'))
1111	CALL check(nf90_put_att(ncid, NF90_GLOBAL, 'comment', '--'))
1112	CALL check(nf90_put_att(ncid, NF90_GLOBAL, 'origin_lat', TRIM(real_to_str(origin_lat*TO_DEGREES, '(F18.13)'))))
1113	CALL check(nf90_put_att(ncid, NF90_GLOBAL, 'origin_lon', TRIM(real_to_str(origin_lon*TO_DEGREES, '(F18.13)'))))
1114	!
1115	!-- FIXME: This is the elevation relative to COSMO-DE/D2 sea level and does
1116	!-- FIXME: not necessarily comply with DHHN2016 (c.f. PALM Input Data
1117	!-- FIXME: Standard v1.9., origin_z)
1118	CALL check(nf90_put_att(ncid, NF90_GLOBAL, 'origin_z', TRIM(real_to_str(z0, '(F18.13)'))))
1119	CALL check(nf90_put_att(ncid, NF90_GLOBAL, 'inifor_version', TRIM(VERSION)))
1120	CALL check(nf90_put_att(ncid, NF90_GLOBAL, 'palm_version', '--'))
1121
1122	!
1123	!
1124	!------------------------------------------------------------------------------
1125	!- Section 2a: Define dimensions for cell centers (scalars in soil and atmosph.)
1126	!------------------------------------------------------------------------------
1127	!
1128	!-- reset dimids first
1129	dimids = (/0, 0, 0/)
1130	CALL check( nf90_def_dim(ncid, "x", nx+1, dimids(1)) )
1131	CALL check( nf90_def_dim(ncid, "y", ny+1, dimids(2)) )
1132	CALL check( nf90_def_dim(ncid, "z", nz, dimids(3)) )
1133	!
1134	!-- save dimids for later
1135	output_file%dimids_scl = dimids
1136
1137	!
1138	!-- reset dimvarids first
1139	dimvarids = (/0, 0, 0/)
1140	CALL check(nf90_def_var(ncid, "x", NF90_FLOAT, dimids(1), dimvarids(1)))
1141	CALL check(nf90_put_att(ncid, dimvarids(1), "standard_name", "x coordinate of cell centers"))
1142	CALL check(nf90_put_att(ncid, dimvarids(1), "units", "m"))
1143
1144	CALL check(nf90_def_var(ncid, "y", NF90_FLOAT, dimids(2), dimvarids(2)))
1145	CALL check(nf90_put_att(ncid, dimvarids(2), "standard_name", "y coordinate of cell centers"))
1146	CALL check(nf90_put_att(ncid, dimvarids(2), "units", "m"))
1147
1148	CALL check(nf90_def_var(ncid, "z", NF90_FLOAT, dimids(3), dimvarids(3)))
1149	CALL check(nf90_put_att(ncid, dimvarids(3), "standard_name", "z coordinate of cell centers"))
1150	CALL check(nf90_put_att(ncid, dimvarids(3), "units", "m"))
1151	!
1152	!-- save dimvarids for later
1153	output_file%dimvarids_scl = dimvarids
1154
1155	!
1156	!-- overwrite third dimid with the one of depth
1157	CALL check(nf90_def_dim(ncid, "zsoil", SIZE(palm_grid%depths), dimids(3)) )
1158	!
1159	!-- save dimids for later
1160	output_file%dimids_soil = dimids
1161
1162	!
1163	!-- overwrite third dimvarid with the one of depth
1164	CALL check(nf90_def_var(ncid, "zsoil", NF90_FLOAT, output_file%dimids_soil(3), dimvarids(3)))
1165	CALL check(nf90_put_att(ncid, dimvarids(3), "standard_name", "depth_below_land"))
1166	CALL check(nf90_put_att(ncid, dimvarids(3), "positive", "down"))
1167	CALL check(nf90_put_att(ncid, dimvarids(3), "units", "m"))
1168	!
1169	!-- save dimvarids for later
1170	output_file%dimvarids_soil = dimvarids
1171	!
1172	!------------------------------------------------------------------------------
1173	!- Section 2b: Define dimensions for cell faces/velocities
1174	!------------------------------------------------------------------------------
1175	!
1176	!-- reset dimids first
1177	dimids = (/0, 0, 0/)
1178	CALL check(nf90_def_dim(ncid, "xu", nx, dimids(1)) )
1179	CALL check(nf90_def_dim(ncid, "yv", ny, dimids(2)) )
1180	CALL check(nf90_def_dim(ncid, "zw", nz-1, dimids(3)) )
1181	!
1182	!-- save dimids for later
1183	output_file%dimids_vel = dimids
1184
1185	!
1186	!-- reset dimvarids first
1187	dimvarids = (/0, 0, 0/)
1188	CALL check(nf90_def_var(ncid, "xu", NF90_FLOAT, dimids(1), dimvarids(1)))
1189	CALL check(nf90_put_att(ncid, dimvarids(1), "standard_name", "x coordinate of cell faces"))
1190	CALL check(nf90_put_att(ncid, dimvarids(1), "units", "m"))
1191
1192	CALL check(nf90_def_var(ncid, "yv", NF90_FLOAT, dimids(2), dimvarids(2)))
1193	CALL check(nf90_put_att(ncid, dimvarids(2), "standard_name", "y coordinate of cell faces"))
1194	CALL check(nf90_put_att(ncid, dimvarids(2), "units", "m"))
1195
1196	CALL check(nf90_def_var(ncid, "zw", NF90_FLOAT, dimids(3), dimvarids(3)))
1197	CALL check(nf90_put_att(ncid, dimvarids(3), "standard_name", "z coordinate of cell faces"))
1198	CALL check(nf90_put_att(ncid, dimvarids(3), "units", "m"))
1199	!
1200	!-- save dimvarids for later
1201	output_file%dimvarids_vel = dimvarids
1202
1203	!
1204	!------------------------------------------------------------------------------
1205	!- Section 2c: Define time dimension
1206	!------------------------------------------------------------------------------
1207	CALL check(nf90_def_dim(ncid, "time", nt, output_file%dimid_time) )
1208	CALL check(nf90_def_var(ncid, "time", NF90_FLOAT, &
1209	output_file%dimid_time, &
1210	output_file%dimvarid_time))
1211	CALL check(nf90_put_att(ncid, output_file%dimvarid_time, "standard_name", "time"))
1212	CALL check(nf90_put_att(ncid, output_file%dimvarid_time, "long_name", "time"))
1213	CALL check(nf90_put_att(ncid, output_file%dimvarid_time, "units", &
1214	"seconds since " // start_date_string // " UTC"))
1215
1216	CALL check(nf90_enddef(ncid))
1217
1218	!
1219	!------------------------------------------------------------------------------
1220	!- Section 3: Write grid coordinates
1221	!------------------------------------------------------------------------------
1222	CALL check(nf90_put_var(ncid, output_file%dimvarids_scl(1), palm_grid%x))
1223	CALL check(nf90_put_var(ncid, output_file%dimvarids_scl(2), palm_grid%y))
1224	CALL check(nf90_put_var(ncid, output_file%dimvarids_scl(3), palm_grid%z))
1225
1226	CALL check(nf90_put_var(ncid, output_file%dimvarids_vel(1), palm_grid%xu))
1227	CALL check(nf90_put_var(ncid, output_file%dimvarids_vel(2), palm_grid%yv))
1228	CALL check(nf90_put_var(ncid, output_file%dimvarids_vel(3), palm_grid%zw))
1229
1230	!
1231	!-- TODO Read in soil depths from input file before this.
1232	CALL check(nf90_put_var(ncid, output_file%dimvarids_soil(3), palm_grid%depths))
1233
1234	!
1235	!-- Write time vector
1236	CALL check(nf90_put_var(ncid, output_file%dimvarid_time, output_file%time))
1237
1238	!
1239	!-- Close the file
1240	CALL check(nf90_close(ncid))
1241
1242	END SUBROUTINE setup_netcdf_dimensions
1243
1244
1245	!------------------------------------------------------------------------------!
1246	! Description:
1247	! ------------
1248	!> Defines the netCDF variables to be written to the dynamic driver file
1249	!------------------------------------------------------------------------------!
1250	SUBROUTINE setup_netcdf_variables(filename, io_group_list)
1251
1252	CHARACTER (LEN=*), INTENT(IN) :: filename
1253	TYPE(io_group), INTENT(IN), TARGET :: io_group_list(:)
1254
1255	TYPE(io_group), POINTER :: group
1256	TYPE(nc_var), POINTER :: var
1257	INTEGER(iwp) :: group_idx, var_idx, n_var
1258	INTEGER :: ncid
1259	LOGICAL :: to_be_written
1260
1261	message = "Defining variables in dynamic driver '" // TRIM(filename) // "'."
1262	CALL report('setup_netcdf_variables', message)
1263
1264	CALL check(nf90_open(TRIM(filename), NF90_WRITE, ncid))
1265	CALL check(nf90_redef(ncid))
1266
1267	n_var = 0
1268	DO group_idx = 1, SIZE( io_group_list )
1269
1270	group => io_group_list(group_idx)
1271	DO var_idx = 1, SIZE( group%out_vars )
1272
1273	to_be_written = .FALSE.
1274
1275	IF (ALLOCATED( group%out_vars )) THEN
1276	var => group%out_vars(var_idx)
1277	n_var = n_var + 1
1278
1279	to_be_written = ( &
1280	group%to_be_processed .AND. var%to_be_processed &
1281	.AND. .NOT. var%is_internal &
1282	)
1283	ENDIF
1284
1285	IF ( to_be_written ) THEN
1286	message = " variable #" // TRIM(str(n_var)) // " '" // TRIM(var%name) // "'."
1287	CALL report('setup_netcdf_variables', message)
1288
1289	CALL netcdf_define_variable(var, ncid)
1290	CALL netcdf_get_dimensions(var, ncid)
1291	ENDIF
1292
1293	ENDDO
1294
1295	ENDDO
1296
1297	CALL check(nf90_enddef(ncid))
1298	CALL check(nf90_close(ncid))
1299
1300	message = "Dynamic driver '" // TRIM(filename) // "' initialized successfully."
1301	CALL report('setup_netcdf_variables', message)
1302
1303	END SUBROUTINE setup_netcdf_variables
1304
1305
1306	!------------------------------------------------------------------------------!
1307	! Description:
1308	! ------------
1309	!> This routine reads and returns all input variables of the given IO group
1310	!> It accomodates the data by allocating a container variable that represents a
1311	!> list of arrays of the same length as the groups variable list. (This list
1312	!> will typically contain one or two items.) After the container, its members
1313	!> are allocated one by one with the appropriate, possibly different,
1314	!> dimensions.
1315	!>
1316	!> The 'group' is an INTENT(INOUT) variable so that 'get_netcdf_variable()' can
1317	!> record netCDF IDs in the 'in_var_list()' member variable.
1318	!------------------------------------------------------------------------------!
1319	SUBROUTINE read_input_variables(group, iter, buffer)
1320	TYPE(io_group), INTENT(INOUT), TARGET :: group
1321	INTEGER(iwp), INTENT(IN) :: iter
1322	TYPE(container), ALLOCATABLE, INTENT(INOUT) :: buffer(:)
1323	INTEGER(iwp) :: hour, buf_id
1324	TYPE(nc_var), POINTER :: input_var
1325	CHARACTER(LEN=PATH), POINTER :: input_file
1326	INTEGER(iwp) :: ivar, nbuffers
1327
1328	message = "Reading data for I/O group '" // TRIM(group%in_var_list(1)%name) // "'."
1329	CALL report('read_input_variables', message)
1330
1331	input_file => group%in_files(iter)
1332
1333	!
1334	!------------------------------------------------------------------------------
1335	!- Section 1: Load input buffers for accumulated variables
1336	!------------------------------------------------------------------------------
1337	!
1338	!-- radiation budgets, precipitation
1339	IF ( group%kind == 'running average' .OR. &
1340	group%kind == 'accumulated' ) THEN
1341
1342	IF ( SIZE( group%in_var_list ) .GT. 1 ) THEN
1343	message = "I/O groups may not contain more than one " // &
1344	"accumulated variable. Group '" // TRIM(group%kind) //&
1345	"' contains " // &
1346	TRIM( str( SIZE( group%in_var_list, kind=iwp ) ) ) // "."
1347	CALL inifor_abort( 'read_input_variables \| accumulation', message )
1348	ENDIF
1349
1350	!
1351	!-- use two buffer arrays
1352	nbuffers = 2
1353	IF ( .NOT. ALLOCATED( buffer ) ) ALLOCATE( buffer(nbuffers) )
1354
1355	!
1356	!-- hour of the day
1357	hour = iter - 1
1358	!
1359	!-- chose correct buffer array
1360	buf_id = select_buffer(hour)
1361
1362	CALL log_runtime('time', 'read')
1363	IF ( ALLOCATED(buffer(buf_id)%array) ) DEALLOCATE(buffer(buf_id)%array)
1364	CALL log_runtime('time', 'alloc')
1365
1366	input_var => group%in_var_list(1)
1367	CALL get_netcdf_variable(input_file, input_var, buffer(buf_id)%array)
1368	CALL report('read_input_variables', "Read accumulated " // TRIM(group%in_var_list(1)%name))
1369
1370	IF ( input_var%is_upside_down ) CALL reverse(buffer(buf_id)%array)
1371	CALL log_runtime('time', 'comp')
1372
1373	!------------------------------------------------------------------------------
1374	!- Section 2: Load input buffers for normal I/O groups
1375	!------------------------------------------------------------------------------
1376	ELSE
1377
1378	!
1379	!-- Allocate one input buffer per output quantity. If more quantities
1380	!-- have to be computed than input variables exist in this group,
1381	!-- allocate more buffers. For instance, in the thermodynamics group,
1382	!-- there are three input variabels (PP, T, Qv) and four quantities
1383	!-- necessart (P, Theta, Rho, qv) for the corresponding output fields
1384	!-- (p0, Theta, qv, ug, and vg)
1385	ALLOCATE( buffer(group%n_output_quantities) )
1386	CALL log_runtime('time', 'alloc')
1387
1388	!
1389	!-- Read in all input variables, leave extra buffers-if any-untouched.
1390	DO ivar = 1, group%n_inputs
1391
1392	input_var => group%in_var_list(ivar)
1393
1394	IF ( input_var%to_be_processed ) THEN
1395	! Check wheather P or PP is present in input file
1396	IF (input_var%name == 'P') THEN
1397	input_var%name = TRIM( get_pressure_varname(input_file) )
1398	CALL log_runtime('time', 'read')
1399	ENDIF
1400
1401	CALL get_netcdf_variable(input_file, input_var, buffer(ivar)%array)
1402
1403	IF ( input_var%is_upside_down ) CALL reverse(buffer(ivar)%array)
1404	CALL log_runtime('time', 'comp')
1405	ENDIF
1406
1407	ENDDO
1408	ENDIF
1409
1410	END SUBROUTINE read_input_variables
1411
1412
1413	!------------------------------------------------------------------------------!
1414	! Description:
1415	! ------------
1416	!> Select the appropriate buffer ID for accumulated COSMO input variables
1417	!> depending on the current hour.
1418	!------------------------------------------------------------------------------!
1419	INTEGER(iwp) FUNCTION select_buffer(hour)
1420	INTEGER(iwp), INTENT(IN) :: hour
1421	INTEGER(iwp) :: step
1422
1423	select_buffer = 0_iwp
1424	step = MODULO(hour, 3_iwp) + 1_iwp
1425
1426	SELECT CASE(step)
1427	CASE(1, 3)
1428	select_buffer = 1_iwp
1429	CASE(2)
1430	select_buffer = 2_iwp
1431	CASE DEFAULT
1432	message = "Invalid step '" // TRIM(str(step))
1433	CALL inifor_abort('select_buffer', message)
1434	END SELECT
1435	END FUNCTION select_buffer
1436
1437
1438	!------------------------------------------------------------------------------!
1439	! Description:
1440	! ------------
1441	!> Checks if the input_file contains the absolute pressure, 'P', or the pressure
1442	!> perturbation, 'PP', and returns the appropriate string.
1443	!------------------------------------------------------------------------------!
1444	CHARACTER(LEN=2) FUNCTION get_pressure_varname(input_file) RESULT(var)
1445	CHARACTER(LEN=*) :: input_file
1446	INTEGER :: ncid, varid
1447
1448	CALL check(nf90_open( TRIM(input_file), NF90_NOWRITE, ncid ))
1449	IF ( nf90_inq_varid( ncid, 'P', varid ) .EQ. NF90_NOERR ) THEN
1450
1451	var = 'P'
1452
1453	ELSE IF ( nf90_inq_varid( ncid, 'PP', varid ) .EQ. NF90_NOERR ) THEN
1454
1455	var = 'PP'
1456	CALL report('get_pressure_var', 'Using PP instead of P')
1457
1458	ELSE
1459
1460	message = "Failed to read '" // TRIM(var) // &
1461	"' from file '" // TRIM(input_file) // "'."
1462	CALL inifor_abort('get_pressure_var', message)
1463
1464	ENDIF
1465
1466	CALL check(nf90_close(ncid))
1467
1468	END FUNCTION get_pressure_varname
1469
1470
1471	!------------------------------------------------------------------------------!
1472	! Description:
1473	! ------------
1474	!> Read the given global attribute form the given netCDF file.
1475	!------------------------------------------------------------------------------!
1476	FUNCTION get_netcdf_attribute(filename, attribute) RESULT(attribute_value)
1477
1478	CHARACTER(LEN=*), INTENT(IN) :: filename, attribute
1479	REAL(wp) :: attribute_value
1480
1481	INTEGER :: ncid
1482
1483	IF ( nf90_open( TRIM(filename), NF90_NOWRITE, ncid ) == NF90_NOERR ) THEN
1484
1485	CALL check(nf90_get_att(ncid, NF90_GLOBAL, TRIM(attribute), attribute_value))
1486	CALL check(nf90_close(ncid))
1487
1488	ELSE
1489
1490	message = "Failed to read '" // TRIM(attribute) // &
1491	"' from file '" // TRIM(filename) // "'."
1492	CALL inifor_abort('get_netcdf_attribute', message)
1493
1494	ENDIF
1495
1496	END FUNCTION get_netcdf_attribute
1497
1498
1499	!------------------------------------------------------------------------------!
1500	! Description:
1501	! ------------
1502	!> Read the attribute of the given variable form the given netCDF file.
1503	!------------------------------------------------------------------------------!
1504	FUNCTION get_netcdf_variable_attribute(filename, varname, attribute) &
1505	RESULT(attribute_value)
1506
1507	CHARACTER(LEN=*), INTENT(IN) :: filename, varname, attribute
1508	REAL(wp) :: attribute_value
1509
1510	INTEGER :: ncid, varid
1511
1512	IF ( nf90_open( TRIM(filename), NF90_NOWRITE, ncid ) == NF90_NOERR ) THEN
1513
1514	CALL check( nf90_inq_varid( ncid, TRIM( varname ), varid ) )
1515	CALL check( nf90_get_att( ncid, varid, TRIM( attribute ), &
1516	attribute_value ) )
1517	CALL check( nf90_close( ncid ) )
1518
1519	ELSE
1520
1521	message = "Failed to read '" // TRIM( varname ) // ":" // &
1522	TRIM( attribute ) // "' from file '" // TRIM(filename) // "'."
1523	CALL inifor_abort('get_netcdf_variable_attribute', message)
1524
1525	ENDIF
1526
1527	END FUNCTION get_netcdf_variable_attribute
1528
1529	!------------------------------------------------------------------------------!
1530	! Description:
1531	! ------------
1532	!> Updates the dynamic driver with the interpolated field of the current
1533	!> variable at the current time step.
1534	!------------------------------------------------------------------------------!
1535	SUBROUTINE update_output(var, array, iter, output_file, cfg)
1536	TYPE(nc_var), INTENT(IN) :: var
1537	REAL(wp), INTENT(IN) :: array(:,:,:)
1538	INTEGER(iwp), INTENT(IN) :: iter
1539	TYPE(nc_file), INTENT(IN) :: output_file
1540	TYPE(inifor_config) :: cfg
1541
1542	INTEGER :: ncid, ndim, start(4), count(4)
1543	LOGICAL :: var_is_time_dependent
1544
1545	var_is_time_dependent = ( &
1546	var%dimids( var%ndim ) == output_file%dimid_time &
1547	)
1548
1549	!
1550	!-- Skip time dimension for output
1551	ndim = var%ndim
1552	IF ( var_is_time_dependent ) ndim = var%ndim - 1
1553
1554	start(:) = (/1,1,1,1/)
1555	start(ndim+1) = iter
1556	count(1:ndim) = var%dimlen(1:ndim)
1557
1558	CALL check(nf90_open(output_file%name, NF90_WRITE, ncid))
1559
1560	!
1561	!-- Reduce dimension of output array according to variable kind
1562	SELECT CASE (TRIM(var%kind))
1563
1564	CASE( 'init scalar profile', 'init u profile', 'init v profile', &
1565	'init w profile', 'init soil profile' )
1566
1567	CALL check(nf90_put_var( ncid, var%varid, array(1,1,:) ) )
1568
1569	CASE( 'init soil', 'init scalar', 'init u', 'init v', 'init w' )
1570
1571	CALL check(nf90_put_var( ncid, var%varid, array(:,:,:) ) )
1572
1573	CASE( 'left scalar', 'right scalar', 'left w', 'right w' )
1574
1575	CALL check(nf90_put_var( ncid, var%varid, array(1,:,:), &
1576	start=start(1:ndim+1), &
1577	count=count(1:ndim) ) )
1578
1579
1580	IF (.NOT. SIZE(array, 2) .EQ. var%dimlen(1)) THEN
1581	PRINT *, "inifor: update_output: Dimension ", 1, " of variable ", &
1582	TRIM(var%name), " (", var%dimlen(1), &
1583	") does not match the dimension of the output array (", &
1584	SIZE(array, 2), ")."
1585	STOP
1586	ENDIF
1587
1588
1589	CASE( 'north scalar', 'south scalar', 'north w', 'south w' )
1590
1591	CALL check(nf90_put_var( ncid, var%varid, array(:,1,:), &
1592	start=start(1:ndim+1), &
1593	count=count(1:ndim) ) )
1594
1595
1596	CASE( 'surface forcing', 'top scalar', 'top w' )
1597
1598	CALL check(nf90_put_var( ncid, var%varid, array(:,:,1), &
1599	start=start(1:ndim+1), &
1600	count=count(1:ndim) ) )
1601
1602	CASE( 'left u', 'right u', 'left v', 'right v' )
1603
1604	CALL check(nf90_put_var( ncid, var%varid, array(1,:,:), &
1605	start=start(1:ndim+1), &
1606	count=count(1:ndim) ) )
1607
1608	CASE( 'north u', 'south u', 'north v', 'south v' )
1609
1610	CALL check(nf90_put_var( ncid, var%varid, array(:,1,:), &
1611	start=start(1:ndim+1), &
1612	count=count(1:ndim) ) )
1613
1614	CASE( 'top u', 'top v' )
1615
1616	CALL check(nf90_put_var( ncid, var%varid, array(:,:,1), &
1617	start=start(1:ndim+1), &
1618	count=count(1:ndim) ) )
1619
1620	CASE( 'time series' )
1621
1622	CALL check(nf90_put_var( ncid, var%varid, array(1,1,1), &
1623	start=start(1:ndim+1) ) )
1624
1625	CASE( 'constant scalar profile', 'geostrophic', &
1626	'left scalar profile', 'right scalar profile', &
1627	'north scalar profile', 'south scalar profile', &
1628	'left u profile', 'right u profile', &
1629	'north u profile', 'south u profile', &
1630	'left v profile', 'right v profile', &
1631	'north v profile', 'south v profile', &
1632	'top scalar profile', 'top u profile', 'top v profile', &
1633	'top w profile', &
1634	'left w profile', 'right w profile', &
1635	'north w profile', 'south w profile' )
1636
1637	CALL check(nf90_put_var( ncid, var%varid, array(1,1,:), &
1638	start=start(1:ndim+1), &
1639	count=count(1:ndim) ) )
1640
1641	CASE ( 'internal profile' )
1642
1643	IF ( cfg%debug ) THEN
1644	CALL check(nf90_put_var( ncid, var%varid, array(1,1,:), &
1645	start=start(1:ndim+1), &
1646	count=count(1:ndim) ) )
1647	ENDIF
1648
1649	CASE ( 'large-scale scalar forcing', 'large-scale w forcing' )
1650
1651	message = "Doing nothing in terms of writing large-scale forings."
1652	CALL report('update_output', message)
1653
1654	CASE DEFAULT
1655
1656	message = "Variable kind '" // TRIM(var%kind) // &
1657	"' not recognized."
1658	CALL inifor_abort('update_output', message)
1659
1660	END SELECT
1661
1662	CALL check(nf90_close(ncid))
1663
1664	END SUBROUTINE update_output
1665
1666
1667	!------------------------------------------------------------------------------!
1668	! Description:
1669	! ------------
1670	!> Checks the status of a netCDF API call and aborts if an error occured
1671	!------------------------------------------------------------------------------!
1672	SUBROUTINE check(status)
1673
1674	INTEGER, INTENT(IN) :: status
1675
1676	IF (status /= nf90_noerr) THEN
1677	message = "NetCDF API call failed with error: " // &
1678	TRIM( nf90_strerror(status) )
1679	CALL inifor_abort('io.check', message)
1680	ENDIF
1681
1682	END SUBROUTINE check
1683
1684
1685	!------------------------------------------------------------------------------!
1686	! Description:
1687	! ------------
1688	!> Setup the origin of the PALM coordinate system based on what is given in the
1689	!> INIFOR namelist file and via an optional static driver.
1690	!------------------------------------------------------------------------------!
1691	SUBROUTINE set_palm_origin( cfg, namelist_longitude, namelist_latitude, &
1692	origin_lon, origin_lat, z0 )
1693
1694	TYPE(inifor_config), INTENT(IN) :: cfg
1695	REAL(wp), INTENT(IN) :: namelist_longitude, namelist_latitude
1696	REAL(wp), INTENT(OUT) :: origin_lon, origin_lat, z0
1697
1698	message = 'Reading PALM-4U origin from'
1699	IF ( TRIM( cfg%static_driver_file ) .NE. '' ) THEN
1700
1701	origin_lon = get_netcdf_attribute( cfg%static_driver_file, TRIM( PIDS_ORIGIN_LON ) )
1702	origin_lat = get_netcdf_attribute( cfg%static_driver_file, TRIM( PIDS_ORIGIN_LAT ) )
1703	z0 = get_netcdf_attribute( cfg%static_driver_file, TRIM( PIDS_ORIGIN_Z ) )
1704
1705	message = TRIM(message) // " static driver file '" &
1706	// TRIM( cfg%static_driver_file ) // "'"
1707
1708
1709	ELSE
1710
1711	origin_lon = namelist_longitude
1712	origin_lat = namelist_latitude
1713
1714	message = TRIM( message ) // " namlist file '" &
1715	// TRIM( cfg%namelist_file ) // "'"
1716
1717	ENDIF
1718	origin_lon = origin_lon * TO_RADIANS
1719	origin_lat = origin_lat * TO_RADIANS
1720
1721	CALL report('set_palm_origin', message)
1722
1723	IF ( cfg%z0_is_set ) THEN
1724	z0 = cfg%z0
1725	IF ( TRIM( cfg%static_driver_file ) .NE. '' ) THEN
1726	message = 'You specified both --static-driver/-t and --elevation/-z0. ' // &
1727	'Using the command line value (' // TRIM( real_to_str_f( cfg%z0 ) ) // &
1728	') instead of static driver value (' // TRIM( real_to_str_f( z0 ) ) // ').'
1729	CALL warn( 'set_palm_origin', message )
1730	ENDIF
1731	ENDIF
1732
1733	END SUBROUTINE set_palm_origin
1734
1735
1736	!------------------------------------------------------------------------------!
1737	! Description:
1738	! ------------
1739	! This function is meant to check weather a COSMO soil variable has an
1740	! additional and redunant surface value at depth = 0.0. For instance operational
1741	! DWD COSMO output contains the surface temperature in T_SO as a copy of the
1742	! values in the first soil layer.
1743	!------------------------------------------------------------------------------!
1744	LOGICAL FUNCTION has_surface_value( soil_var, filename )
1745
1746	TYPE(nc_var), INTENT(IN) :: soil_var
1747	CHARACTER(LEN=*), INTENT(IN) :: filename
1748
1749	REAL(wp), ALLOCATABLE :: depths(:)
1750
1751	CALL get_dimension_vector_of_variable( &
1752	soil_var%name, &
1753	dim_idx = NC_DEPTH_DIM_IDX, &
1754	filename = filename, &
1755	dim_vector = depths &
1756	)
1757
1758	has_surface_value = nearly_equal( depths(1), 0.0_wp, 10 * EPSILON(1.0_wp) )
1759
1760	END FUNCTION has_surface_value
1761
1762
1763	!------------------------------------------------------------------------------!
1764	! Description:
1765	! ------------
1766	! This routine reads the dim_idx-th dimension vector of the variable varname
1767	! from netCDF file filename. It is used for finding the depth coordinate vector
1768	! of COSMO soil variables without knowing its name.
1769	!------------------------------------------------------------------------------!
1770	SUBROUTINE get_dimension_vector_of_variable( varname, dim_idx, filename, dim_vector )
1771	CHARACTER(LEN=*), INTENT(IN) :: varname, filename
1772	INTEGER, INTENT(IN) :: dim_idx
1773
1774	REAL(wp), INTENT(OUT), ALLOCATABLE :: dim_vector(:)
1775
1776	INTEGER :: dimids(NF90_MAX_VAR_DIMS)
1777	INTEGER :: varid
1778	CHARACTER(LEN=NF90_MAX_NAME) :: dimname
1779
1780	INTEGER :: ncid
1781
1782	IF ( nf90_open( TRIM( filename ), NF90_NOWRITE, ncid ) .EQ. NF90_NOERR ) THEN
1783
1784	!
1785	!-- get id of variable varname
1786	CALL check( nf90_inq_varid( ncid, TRIM( varname ), varid ) )
1787
1788	!
1789	!-- get dimension ids of variable with varid
1790	CALL check( nf90_inquire_variable( ncid, varid, dimids = dimids ) )
1791
1792	!
1793	!-- get name of dim_idx-th dimension variable
1794	CALL check( nf90_inquire_dimension( ncid, dimids(dim_idx), name = dimname ) )
1795	CALL check( nf90_close( ncid ) )
1796
1797	ELSE
1798
1799	message = "Failed to open file '" // TRIM(filename) // "'."
1800	CALL inifor_abort('get_netcdf_variable', message)
1801
1802	ENDIF
1803
1804	! get dimension vector with dimname
1805	CALL get_netcdf_dim_vector( filename, dimname, dim_vector )
1806
1807	END SUBROUTINE get_dimension_vector_of_variable
1808
1809
1810	LOGICAL FUNCTION netcdf_variable_present_in_file( varname, filename )
1811	CHARACTER(LEN=SNAME), INTENT(IN) :: varname
1812	CHARACTER(LEN=PATH), INTENT(IN) :: filename
1813
1814	INTEGER :: ncid, varid
1815
1816	netcdf_variable_present_in_file = ( &
1817	nf90_open( TRIM( filename ), NF90_NOWRITE, ncid ) .EQ. NF90_NOERR .AND. &
1818	nf90_inq_varid( ncid, varname, varid ) .EQ. NF90_NOERR .AND. &
1819	nf90_close( ncid ) .EQ. NF90_NOERR &
1820	)
1821
1822	END FUNCTION netcdf_variable_present_in_file
1823
1824
1825	END MODULE inifor_io
1826	#endif

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |