Home

Context Navigation

source: palm/trunk/SOURCE/virtual_flight_mod.f90 @ 4540

Last change on this file since 4540 was 4535, checked in by raasch, 4 years ago
bugfix for restart data format query
Property svn:keywords set to `Id`
File size: 47.0 KB

Line
1	!> @file virtual_flights_mod.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 terms of the GNU General
6	! Public License as published by the Free Software Foundation, either version 3 of the License, or
7	! (at your option) any later version.
8	!
9	! PALM is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the
10	! implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General
11	! Public License for more details.
12	!
13	! You should have received a copy of the GNU General Public License along with PALM. If not, see
14	! <http://www.gnu.org/licenses/>.
15	!
16	! Copyright 1997-2020 Leibniz Universitaet Hannover
17	!--------------------------------------------------------------------------------------------------!
18	!
19	!
20	! Current revisions:
21	! -----------------
22	!
23	!
24	! Former revisions:
25	! -----------------
26	! $Id: virtual_flight_mod.f90 4535 2020-05-15 12:07:23Z raasch $
27	! bugfix for restart data format query
28	!
29	! 4522 2020-05-06 14:17:05Z suehring
30	! Modularize user_init_flight in order to provide an explicit interface.
31	!
32	! 4497 2020-04-15 10:20:51Z raasch
33	! file re-formatted to follow the PALM coding standard
34	!
35	!
36	! 4495 2020-04-13 20:11:20Z raasch
37	! restart data handling with MPI-IO added
38	!
39	! 4360 2020-01-07 11:25:50Z suehring
40	! Corrected "Former revisions" section
41	!
42	! 4004 2019-05-24 11:32:38Z suehring
43	! Allow variable start- and end locations also in return mode
44	!
45	! 3885 2019-04-11 11:29:34Z kanani
46	! Changes related to global restructuring of location messages and introduction of additional
47	! debug messages
48	!
49	! 3655 2019-01-07 16:51:22Z knoop
50	! variables documented
51	!
52	! 1957 2016-07-07 10:43:48Z suehring
53	! Initial revision
54	!
55	! Description:
56	! ------------
57	!> Module for virtual flight measurements.
58	!> @todo Err msg PA0438: flight can be inside topography -> extra check?
59	!--------------------------------------------------------------------------------------------------!
60	MODULE flight_mod
61
62	USE control_parameters, &
63	ONLY: debug_output, &
64	fl_max, num_leg, &
65	num_var_fl, &
66	num_var_fl_user, &
67	restart_data_format_output, &
68	virtual_flight
69
70	USE kinds
71
72	USE restart_data_mpi_io_mod, &
73	ONLY: rd_mpi_io_check_array, rrd_mpi_io_global_array, wrd_mpi_io_global_array
74
75	USE user_init_flight_mod, &
76	ONLY: user_init_flight
77
78
79	CHARACTER(LEN=6), DIMENSION(fl_max) :: leg_mode = 'cyclic' !< flight mode through the model domain, either 'cyclic' or
80	!<'return'
81
82	INTEGER(iwp) :: l !< index for flight leg
83	INTEGER(iwp) :: var_index !< index for measured variable
84
85	LOGICAL, DIMENSION(:), ALLOCATABLE :: cyclic_leg !< flag to identify fly mode
86
87	REAL(wp) :: flight_end = 9999999.9_wp !< end time of virtual flight
88	REAL(wp) :: flight_begin = 0.0_wp !< end time of virtual flight
89
90	REAL(wp), DIMENSION(fl_max) :: flight_angle = 45.0_wp !< angle determining the horizontal flight direction
91	REAL(wp), DIMENSION(fl_max) :: flight_level = 100.0_wp !< flight level
92	REAL(wp), DIMENSION(fl_max) :: max_elev_change = 0.0_wp !< maximum elevation change for the respective flight leg
93	REAL(wp), DIMENSION(fl_max) :: rate_of_climb = 0.0_wp !< rate of climb or descent
94	REAL(wp), DIMENSION(fl_max) :: speed_agl = 25.0_wp !< absolute horizontal flight speed above ground level (agl)
95	REAL(wp), DIMENSION(fl_max) :: x_start = 999999999.0_wp !< start x position
96	REAL(wp), DIMENSION(fl_max) :: x_end = 999999999.0_wp !< end x position
97	REAL(wp), DIMENSION(fl_max) :: y_start = 999999999.0_wp !< start y position
98	REAL(wp), DIMENSION(fl_max) :: y_end = 999999999.0_wp !< end y position
99
100	REAL(wp), DIMENSION(:), ALLOCATABLE :: u_agl !< u-component of flight speed
101	REAL(wp), DIMENSION(:), ALLOCATABLE :: v_agl !< v-component of flight speed
102	REAL(wp), DIMENSION(:), ALLOCATABLE :: w_agl !< w-component of flight speed
103	REAL(wp), DIMENSION(:), ALLOCATABLE :: x_pos !< aircraft x-position
104	REAL(wp), DIMENSION(:), ALLOCATABLE :: y_pos !< aircraft y-position
105	REAL(wp), DIMENSION(:), ALLOCATABLE :: z_pos !< aircraft z-position
106
107	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: sensor_l !< measured data on local PE
108	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: sensor !< measured data
109
110	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: var_u !< dummy array for possibly user-defined quantities
111
112	SAVE
113
114	PRIVATE
115
116	INTERFACE flight_header
117	MODULE PROCEDURE flight_header
118	END INTERFACE flight_header
119
120	INTERFACE flight_init
121	MODULE PROCEDURE flight_init
122	END INTERFACE flight_init
123
124	INTERFACE flight_init_output
125	MODULE PROCEDURE flight_init_output
126	END INTERFACE flight_init_output
127
128	INTERFACE flight_check_parameters
129	MODULE PROCEDURE flight_check_parameters
130	END INTERFACE flight_check_parameters
131
132	INTERFACE flight_parin
133	MODULE PROCEDURE flight_parin
134	END INTERFACE flight_parin
135
136	INTERFACE interpolate_xyz
137	MODULE PROCEDURE interpolate_xyz
138	END INTERFACE interpolate_xyz
139
140	INTERFACE flight_measurement
141	MODULE PROCEDURE flight_measurement
142	END INTERFACE flight_measurement
143
144	INTERFACE flight_rrd_global
145	MODULE PROCEDURE flight_rrd_global_ftn
146	MODULE PROCEDURE flight_rrd_global_mpi
147	END INTERFACE flight_rrd_global
148
149	INTERFACE flight_wrd_global
150	MODULE PROCEDURE flight_wrd_global
151	END INTERFACE flight_wrd_global
152
153	!
154	!-- Private interfaces
155	PRIVATE flight_check_parameters, &
156	flight_init_output, &
157	interpolate_xyz
158	!
159	!-- Public interfaces
160	PUBLIC flight_init, &
161	flight_header, &
162	flight_parin, &
163	flight_measurement, &
164	flight_wrd_global, &
165	flight_rrd_global
166	!
167	!-- Public variables
168	PUBLIC fl_max, &
169	sensor, &
170	x_pos, &
171	y_pos, &
172	z_pos
173
174	CONTAINS
175
176	!--------------------------------------------------------------------------------------------------!
177	! Description:
178	! ------------
179	!> Header output for flight module.
180	!--------------------------------------------------------------------------------------------------!
181	SUBROUTINE flight_header ( io )
182
183
184	IMPLICIT NONE
185
186	INTEGER(iwp), INTENT(IN) :: io !< Unit of the output file
187
188	WRITE ( io, 1 )
189	WRITE ( io, 2 )
190	WRITE ( io, 3 ) num_leg
191	WRITE ( io, 4 ) flight_begin
192	WRITE ( io, 5 ) flight_end
193
194	DO l=1, num_leg
195	WRITE ( io, 6 ) l
196	WRITE ( io, 7 ) speed_agl(l)
197	WRITE ( io, 8 ) flight_level(l)
198	WRITE ( io, 9 ) max_elev_change(l)
199	WRITE ( io, 10 ) rate_of_climb(l)
200	WRITE ( io, 11 ) leg_mode(l)
201	ENDDO
202
203
204	1 FORMAT (' Virtual flights: ----------------' )
205	2 FORMAT (' Output every timestep' )
206	3 FORMAT (' Number of flight legs:', I3 )
207	4 FORMAT (' Begin of measurements:', F10.1 , ' s' )
208	5 FORMAT (' End of measurements:', F10.1 , ' s' )
209	6 FORMAT (' Leg', I3/, ' ------' )
210	7 FORMAT (' Flight speed : ', F5.1, ' m/s' )
211	8 FORMAT (' Flight level : ', F5.1, ' m' )
212	9 FORMAT (' Maximum elevation change: ', F5.1, ' m/s' )
213	10 FORMAT (' Rate of climb / descent : ', F5.1, ' m/s' )
214	11 FORMAT (' Leg mode : ', A/ )
215
216	END SUBROUTINE flight_header
217
218	!--------------------------------------------------------------------------------------------------!
219	! Description:
220	! ------------
221	!> Reads the namelist flight_par.
222	!--------------------------------------------------------------------------------------------------!
223	SUBROUTINE flight_parin
224
225	USE control_parameters, &
226	ONLY: message_string
227
228	IMPLICIT NONE
229
230	CHARACTER(LEN=80) :: line !< dummy string that contains the current line of the parameter file
231
232	NAMELIST /flight_par/ flight_angle, &
233	flight_begin, &
234	flight_end, &
235	flight_level, &
236	leg_mode, &
237	max_elev_change, &
238	rate_of_climb, &
239	speed_agl, &
240	x_end, &
241	x_start, &
242	y_end, &
243	y_start
244
245
246	NAMELIST /virtual_flight_parameters/ flight_angle, &
247	flight_begin, &
248	flight_end, &
249	flight_level, &
250	leg_mode, &
251	max_elev_change, &
252	rate_of_climb, &
253	speed_agl, &
254	x_end, &
255	x_start, &
256	y_end, &
257	y_start
258	!
259	!-- Try to find the namelist flight_par
260	REWIND ( 11 )
261	line = ' '
262	DO WHILE ( INDEX( line, '&virtual_flight_parameters' ) == 0 )
263	READ ( 11, '(A)', END = 12 ) line
264	ENDDO
265	BACKSPACE ( 11 )
266
267	!
268	!-- Read namelist
269	READ ( 11, virtual_flight_parameters, ERR = 10 )
270	!
271	!-- Set switch so that virtual flights shall be carried out
272	virtual_flight = .TRUE.
273
274	GOTO 14
275
276	10 BACKSPACE( 11 )
277	READ( 11 , '(A)') line
278	CALL parin_fail_message( 'virtual_flight_parameters', line )
279	!
280	!-- Try to find the old namelist
281	12 REWIND ( 11 )
282	line = ' '
283	DO WHILE ( INDEX( line, '&flight_par' ) == 0 )
284	READ ( 11, '(A)', END = 14 ) line
285	ENDDO
286	BACKSPACE ( 11 )
287
288	!
289	!-- Read namelist
290	READ ( 11, flight_par, ERR = 13, END = 14 )
291
292	message_string = 'namelist flight_par is deprecated and will be ' // &
293	'removed in near future.& Please use namelist ' // &
294	'virtual_flight_parameters instead'
295	CALL message( 'flight_parin', 'PA0487', 0, 1, 0, 6, 0 )
296	!
297	!-- Set switch so that virtual flights shall be carried out
298	virtual_flight = .TRUE.
299
300	GOTO 14
301
302	13 BACKSPACE( 11 )
303	READ( 11 , '(A)') line
304	CALL parin_fail_message( 'flight_par', line )
305
306	14 CONTINUE
307
308	END SUBROUTINE flight_parin
309
310	!--------------------------------------------------------------------------------------------------!
311	! Description:
312	! ------------
313	!> Inititalization of required arrays, number of legs and flags. Moreover, initialize flight speed
314	!> and -direction, as well as start positions.
315	!--------------------------------------------------------------------------------------------------!
316	SUBROUTINE flight_init
317
318	USE basic_constants_and_equations_mod, &
319	ONLY: pi
320
321	USE control_parameters, &
322	ONLY: initializing_actions
323
324	USE indices, &
325	ONLY: nxlg, &
326	nxrg, &
327	nysg, &
328	nyng, &
329	nzb, &
330	nzt
331
332	IMPLICIT NONE
333
334	REAL(wp) :: distance !< distance between start and end position of a flight leg
335
336
337	IF ( debug_output ) CALL debug_message( 'flight_init', 'start' )
338	!
339	!-- Determine the number of flight legs
340	l = 1
341	DO WHILE ( x_start(l) /= 999999999.0_wp .AND. l <= SIZE(x_start) )
342	l = l + 1
343	ENDDO
344	num_leg = l-1
345	!
346	!-- Check for proper parameter settings
347	CALL flight_check_parameters
348	!
349	!-- Allocate and initialize logical array for flight pattern
350	ALLOCATE( cyclic_leg(1:num_leg) )
351	!
352	!-- Initialize flags for cyclic/return legs
353	DO l = 1, num_leg
354	cyclic_leg(l) = MERGE( .TRUE., .FALSE., TRIM( leg_mode(l) ) == 'cyclic' )
355	ENDDO
356	!
357	!-- Allocate and initialize arraxs for flight position and speed. In case of restart runs these data
358	!-- are read by the routine read_flight_restart_data instead.
359	IF ( TRIM( initializing_actions ) /= 'read_restart_data' ) THEN
360
361	ALLOCATE( x_pos(1:num_leg), y_pos(1:num_leg ), z_pos(1:num_leg) )
362	!
363	!-- Inititalize x-, y-, and z-positions with initial start position
364	x_pos(1:num_leg) = x_start(1:num_leg)
365	y_pos(1:num_leg) = y_start(1:num_leg)
366	z_pos(1:num_leg) = flight_level(1:num_leg)
367	!
368	!-- Allocate arrays for flight-speed components
369	ALLOCATE( u_agl(1:num_leg), &
370	v_agl(1:num_leg), &
371	w_agl(1:num_leg) )
372	!
373	!-- Inititalize u-, v- and w-component.
374	DO l = 1, num_leg
375	!
376	!-- In case of return-legs, the flight direction, i.e. the horizontal flight-speed components,
377	!-- are derived from the given start/end positions.
378	IF ( .NOT. cyclic_leg(l) ) THEN
379	distance = SQRT( ( x_end(l) - x_start(l) )2 + ( y_end(l) - y_start(l) )2 )
380	u_agl(l) = speed_agl(l) * ( x_end(l) - x_start(l) ) / distance
381	v_agl(l) = speed_agl(l) * ( y_end(l) - y_start(l) ) / distance
382	w_agl(l) = rate_of_climb(l)
383	!
384	!-- In case of cyclic-legs, flight direction is directly derived from the given flight angle.
385	ELSE
386	u_agl(l) = speed_agl(l) * COS( flight_angle(l) * pi / 180.0_wp )
387	v_agl(l) = speed_agl(l) * SIN( flight_angle(l) * pi / 180.0_wp )
388	w_agl(l) = rate_of_climb(l)
389	ENDIF
390
391	ENDDO
392
393	ENDIF
394	!
395	!-- Initialized data output
396	CALL flight_init_output
397	!
398	!-- Allocate array required for user-defined quantities if necessary.
399	IF ( num_var_fl_user > 0 ) ALLOCATE( var_u(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
400	!
401	!-- Allocate and initialize arrays containing the measured data
402	ALLOCATE( sensor_l(1:num_var_fl,1:num_leg) )
403	ALLOCATE( sensor(1:num_var_fl,1:num_leg) )
404	sensor_l = 0.0_wp
405	sensor = 0.0_wp
406
407	IF ( debug_output ) CALL debug_message( 'flight_init', 'end' )
408
409	END SUBROUTINE flight_init
410
411
412	!--------------------------------------------------------------------------------------------------!
413	! Description:
414	! ------------
415	!> Initialization of output-variable names and units.
416	!--------------------------------------------------------------------------------------------------!
417	SUBROUTINE flight_init_output
418
419	USE control_parameters, &
420	ONLY: cloud_droplets, &
421	humidity, &
422	neutral, &
423	passive_scalar
424
425	USE bulk_cloud_model_mod, &
426	ONLY: bulk_cloud_model
427
428	USE netcdf_interface
429
430	IMPLICIT NONE
431
432	CHARACTER(LEN=10) :: label_leg !< dummy argument to convert integer to string
433
434	INTEGER(iwp) :: i !< loop variable
435	INTEGER(iwp) :: id_pt !< identifyer for labeling
436	INTEGER(iwp) :: id_q !< identifyer for labeling
437	INTEGER(iwp) :: id_ql !< identifyer for labeling
438	INTEGER(iwp) :: id_s !< identifyer for labeling
439	INTEGER(iwp) :: id_u = 1 !< identifyer for labeling
440	INTEGER(iwp) :: id_v = 2 !< identifyer for labeling
441	INTEGER(iwp) :: id_w = 3 !< identifyer for labeling
442	INTEGER(iwp) :: k !< index variable
443
444	LOGICAL :: init = .TRUE. !< flag to distiquish calls of user_init_flight
445	!
446	!-- Define output quanities, at least three variables are measured (u,v,w)
447	num_var_fl = 3
448	IF ( .NOT. neutral ) THEN
449	num_var_fl = num_var_fl + 1
450	id_pt = num_var_fl
451	ENDIF
452	IF ( humidity ) THEN
453	num_var_fl = num_var_fl + 1
454	id_q = num_var_fl
455	ENDIF
456	IF ( bulk_cloud_model .OR. cloud_droplets ) THEN
457	num_var_fl = num_var_fl + 1
458	id_ql = num_var_fl
459	ENDIF
460	IF ( passive_scalar ) THEN
461	num_var_fl = num_var_fl + 1
462	id_s = num_var_fl
463	ENDIF
464	!
465	!-- Write output strings for dimensions x, y, z
466	DO l=1, num_leg
467
468	IF ( l < 10 ) WRITE( label_leg, '(I1)' ) l
469	IF ( l >= 10 .AND. l < 100 ) WRITE( label_leg, '(I2)' ) l
470	IF ( l >= 100 .AND. l < 1000 ) WRITE( label_leg, '(I3)' ) l
471
472	dofl_dim_label_x(l) = 'x_' // TRIM( label_leg )
473	dofl_dim_label_y(l) = 'y_' // TRIM( label_leg )
474	dofl_dim_label_z(l) = 'z_' // TRIM( label_leg )
475
476	ENDDO
477
478	!
479	!-- Call user routine to initialize further variables
480	CALL user_init_flight( init )
481	!
482	!-- Write output labels and units for the quanities
483	k = 1
484	DO l=1, num_leg
485
486	IF ( l < 10 ) WRITE( label_leg, '(I1)' ) l
487	IF ( l >= 10 .AND. l < 100 ) WRITE( label_leg, '(I2)' ) l
488	IF ( l >= 100 .AND. l < 1000 ) WRITE( label_leg, '(I3)' ) l
489
490	label_leg = 'leg_' // TRIM(label_leg)
491	DO i=1, num_var_fl
492
493	IF ( i == id_u ) THEN
494	dofl_label(k) = TRIM( label_leg ) // '_u'
495	dofl_unit(k) = 'm/s'
496	k = k + 1
497	ELSEIF ( i == id_v ) THEN
498	dofl_label(k) = TRIM( label_leg ) // '_v'
499	dofl_unit(k) = 'm/s'
500	k = k + 1
501	ELSEIF ( i == id_w ) THEN
502	dofl_label(k) = TRIM( label_leg ) // '_w'
503	dofl_unit(k) = 'm/s'
504	k = k + 1
505	ELSEIF ( i == id_pt ) THEN
506	dofl_label(k) = TRIM( label_leg ) // '_theta'
507	dofl_unit(k) = 'K'
508	k = k + 1
509	ELSEIF ( i == id_q ) THEN
510	dofl_label(k) = TRIM( label_leg ) // '_q'
511	dofl_unit(k) = 'kg/kg'
512	k = k + 1
513	ELSEIF ( i == id_ql ) THEN
514	dofl_label(k) = TRIM( label_leg ) // '_ql'
515	dofl_unit(k) = 'kg/kg'
516	k = k + 1
517	ELSEIF ( i == id_s ) THEN
518	dofl_label(k) = TRIM( label_leg ) // '_s'
519	dofl_unit(k) = 'kg/kg'
520	k = k + 1
521	ENDIF
522	ENDDO
523
524	DO i=1, num_var_fl_user
525	CALL user_init_flight( init, k, i, label_leg )
526	ENDDO
527
528	ENDDO
529	!
530	!-- Finally, set the total number of flight-output quantities.
531	num_var_fl = num_var_fl + num_var_fl_user
532
533	END SUBROUTINE flight_init_output
534
535	!--------------------------------------------------------------------------------------------------!
536	! Description:
537	! ------------
538	!> This routine calculates the current flight positions and calls the respective interpolation
539	!> routine to measure the data at the current flight position.
540	!--------------------------------------------------------------------------------------------------!
541	SUBROUTINE flight_measurement
542
543	USE arrays_3d, &
544	ONLY: ddzu, &
545	ddzw, &
546	pt, &
547	q, &
548	ql, &
549	s, &
550	u, &
551	v, &
552	w, &
553	zu, &
554	zw
555
556	USE control_parameters, &
557	ONLY: cloud_droplets, &
558	dt_3d, &
559	humidity, &
560	neutral, &
561	passive_scalar, &
562	simulated_time
563
564	USE cpulog, &
565	ONLY: cpu_log, &
566	log_point
567
568	USE grid_variables, &
569	ONLY: ddx, &
570	ddy, &
571	dx, &
572	dy
573
574	USE indices, &
575	ONLY: nx, &
576	nxl, &
577	nxr, &
578	ny, &
579	nys, &
580	nyn
581
582	USE bulk_cloud_model_mod, &
583	ONLY: bulk_cloud_model
584
585	USE pegrid
586
587	IMPLICIT NONE
588
589	INTEGER(iwp) :: i !< index of current grid box along x
590	INTEGER(iwp) :: j !< index of current grid box along y
591	INTEGER(iwp) :: n !< loop variable for number of user-defined output quantities
592
593	LOGICAL :: on_pe !< flag to check if current flight position is on current PE
594
595	REAL(wp) :: x !< distance between left edge of current grid box and flight position
596	REAL(wp) :: y !< distance between south edge of current grid box and flight position
597
598	CALL cpu_log( log_point(65), 'flight_measurement', 'start' )
599	!
600	!-- Perform flight measurement if start time is reached.
601	IF ( simulated_time >= flight_begin .AND. simulated_time <= flight_end ) THEN
602
603	sensor_l = 0.0_wp
604	sensor = 0.0_wp
605	!
606	!-- Loop over all flight legs
607	DO l = 1, num_leg
608	!
609	!-- Update location for each leg
610	x_pos(l) = x_pos(l) + u_agl(l) * dt_3d
611	y_pos(l) = y_pos(l) + v_agl(l) * dt_3d
612	z_pos(l) = z_pos(l) + w_agl(l) * dt_3d
613	!
614	!-- Check if location must be modified for return legs.
615	!-- Carry out horizontal reflection if required.
616	IF ( .NOT. cyclic_leg(l) ) THEN
617
618	IF ( x_start(l) <= x_end(l) ) THEN
619	!
620	!-- Outward flight, i.e. from start to end
621	IF ( u_agl(l) >= 0.0_wp .AND. x_pos(l) > x_end(l) ) THEN
622	x_pos(l) = 2.0_wp * x_end(l) - x_pos(l)
623	u_agl(l) = - u_agl(l)
624	!
625	!-- Return flight, i.e. from end to start
626	ELSEIF ( u_agl(l) < 0.0_wp .AND. x_pos(l) < x_start(l) ) THEN
627	x_pos(l) = 2.0_wp * x_start(l) - x_pos(l)
628	u_agl(l) = - u_agl(l)
629	ENDIF
630	ELSE
631	!
632	!-- Outward flight, i.e. from start to end
633	IF ( u_agl(l) < 0.0_wp .AND. x_pos(l) < x_end(l) ) THEN
634	x_pos(l) = 2.0_wp * x_end(l) - x_pos(l)
635	u_agl(l) = - u_agl(l)
636	!
637	!-- Return flight, i.e. from end to start
638	ELSEIF ( u_agl(l) >= 0.0_wp .AND. x_pos(l) > x_start(l) ) THEN
639	x_pos(l) = 2.0_wp * x_start(l) - x_pos(l)
640	u_agl(l) = - u_agl(l)
641	ENDIF
642	ENDIF
643
644	IF ( y_start(l) <= y_end(l) ) THEN
645	!
646	!-- Outward flight, i.e. from start to end
647	IF ( v_agl(l) >= 0.0_wp .AND. y_pos(l) > y_end(l) ) THEN
648	y_pos(l) = 2.0_wp * y_end(l) - y_pos(l)
649	v_agl(l) = - v_agl(l)
650	!
651	!-- Return flight, i.e. from end to start
652	ELSEIF ( v_agl(l) < 0.0_wp .AND. y_pos(l) < y_start(l) ) THEN
653	y_pos(l) = 2.0_wp * y_start(l) - y_pos(l)
654	v_agl(l) = - v_agl(l)
655	ENDIF
656	ELSE
657	!
658	!-- Outward flight, i.e. from start to end
659	IF ( v_agl(l) < 0.0_wp .AND. y_pos(l) < y_end(l) ) THEN
660	y_pos(l) = 2.0_wp * y_end(l) - y_pos(l)
661	v_agl(l) = - v_agl(l)
662	!
663	!-- Return flight, i.e. from end to start
664	ELSEIF ( v_agl(l) >= 0.0_wp .AND. y_pos(l) > y_start(l) ) THEN
665	y_pos(l) = 2.0_wp * y_start(l) - y_pos(l)
666	v_agl(l) = - v_agl(l)
667	ENDIF
668	ENDIF
669	!
670	!-- Check if flight position is outside the model domain and apply cyclic conditions if required
671	ELSEIF ( cyclic_leg(l) ) THEN
672	!
673	!-- Check if aircraft leaves the model domain at the right boundary
674	IF ( ( flight_angle(l) >= 0.0_wp .AND. &
675	flight_angle(l) <= 90.0_wp ) .OR. &
676	( flight_angle(l) >= 270.0_wp .AND. &
677	flight_angle(l) <= 360.0_wp ) ) THEN
678	IF ( x_pos(l) >= ( nx + 0.5_wp ) * dx ) &
679	x_pos(l) = x_pos(l) - ( nx + 1 ) * dx
680	!
681	!-- Check if aircraft leaves the model domain at the left boundary
682	ELSEIF ( flight_angle(l) > 90.0_wp .AND. flight_angle(l) < 270.0_wp ) THEN
683	IF ( x_pos(l) < -0.5_wp * dx ) &
684	x_pos(l) = ( nx + 1 ) * dx + x_pos(l)
685	ENDIF
686	!
687	!-- Check if aircraft leaves the model domain at the north boundary
688	IF ( flight_angle(l) >= 0.0_wp .AND. flight_angle(l) <= 180.0_wp ) THEN
689	IF ( y_pos(l) >= ( ny + 0.5_wp ) * dy ) &
690	y_pos(l) = y_pos(l) - ( ny + 1 ) * dy
691	!
692	!-- Check if aircraft leaves the model domain at the south boundary
693	ELSEIF ( flight_angle(l) > 180.0_wp .AND. flight_angle(l) < 360.0_wp ) THEN
694	IF ( y_pos(l) < -0.5_wp * dy ) &
695	y_pos(l) = ( ny + 1 ) * dy + y_pos(l)
696	ENDIF
697
698	ENDIF
699	!
700	!-- Check if maximum elevation change is already reached. If required reflect vertically.
701	IF ( rate_of_climb(l) /= 0.0_wp ) THEN
702	!
703	!-- First check if aircraft is too high
704	IF ( w_agl(l) > 0.0_wp .AND. z_pos(l) - flight_level(l) > max_elev_change(l) ) THEN
705	z_pos(l) = 2.0_wp * ( flight_level(l) + max_elev_change(l) ) - z_pos(l)
706	w_agl(l) = - w_agl(l)
707	!
708	!-- Check if aircraft is too low
709	ELSEIF ( w_agl(l) < 0.0_wp .AND. z_pos(l) < flight_level(l) ) THEN
710	z_pos(l) = 2.0_wp * flight_level(l) - z_pos(l)
711	w_agl(l) = - w_agl(l)
712	ENDIF
713
714	ENDIF
715	!
716	!-- Determine grid indices for flight position along x- and y-direction. Please note, there is
717	!-- a special treatment for the index along z-direction, which is due to vertical grid stretching.
718	i = ( x_pos(l) + 0.5_wp * dx ) * ddx
719	j = ( y_pos(l) + 0.5_wp * dy ) * ddy
720	!
721	!-- Check if indices are on current PE
722	on_pe = ( i >= nxl .AND. i <= nxr .AND. j >= nys .AND. j <= nyn )
723
724	IF ( on_pe ) THEN
725
726	var_index = 1
727	!
728	!-- Recalculate indices, as u is shifted by -0.5*dx.
729	i = x_pos(l) * ddx
730	j = ( y_pos(l) + 0.5_wp * dy ) * ddy
731	!
732	!-- Calculate distance from left and south grid-box edges.
733	x = x_pos(l) - ( 0.5_wp - i ) * dx
734	y = y_pos(l) - j * dy
735	!
736	!-- Interpolate u-component onto current flight position.
737	CALL interpolate_xyz( u, zu, ddzu, 1.0_wp, x, y, var_index, j, i )
738	var_index = var_index + 1
739	!
740	!-- Recalculate indices, as v is shifted by -0.5*dy.
741	i = ( x_pos(l) + 0.5_wp * dx ) * ddx
742	j = y_pos(l) * ddy
743
744	x = x_pos(l) - i * dx
745	y = y_pos(l) - ( 0.5_wp - j ) * dy
746	CALL interpolate_xyz( v, zu, ddzu, 1.0_wp, x, y, var_index, j, i )
747	var_index = var_index + 1
748	!
749	!-- Interpolate w and scalar quantities. Recalculate indices.
750	i = ( x_pos(l) + 0.5_wp * dx ) * ddx
751	j = ( y_pos(l) + 0.5_wp * dy ) * ddy
752	x = x_pos(l) - i * dx
753	y = y_pos(l) - j * dy
754	!
755	!-- Interpolate w-velocity component.
756	CALL interpolate_xyz( w, zw, ddzw, 0.0_wp, x, y, var_index, j, i )
757	var_index = var_index + 1
758	!
759	!-- Potential temerature
760	IF ( .NOT. neutral ) THEN
761	CALL interpolate_xyz( pt, zu, ddzu, 1.0_wp, x, y, var_index, j, i )
762	var_index = var_index + 1
763	ENDIF
764	!
765	!-- Humidity
766	IF ( humidity ) THEN
767	CALL interpolate_xyz( q, zu, ddzu, 1.0_wp, x, y, var_index, j, i )
768	var_index = var_index + 1
769	ENDIF
770	!
771	!-- Liquid water content
772	IF ( bulk_cloud_model .OR. cloud_droplets ) THEN
773	CALL interpolate_xyz( ql, zu, ddzu, 1.0_wp, x, y, var_index, j, i )
774	var_index = var_index + 1
775	ENDIF
776	!
777	!-- Passive scalar
778	IF ( passive_scalar ) THEN
779	CALL interpolate_xyz( s, zu, ddzu, 1.0_wp, x, y, var_index, j, i )
780	var_index = var_index + 1
781	ENDIF
782	!
783	!-- Treat user-defined variables if required
784	DO n = 1, num_var_fl_user
785	CALL user_flight( var_u, n )
786	CALL interpolate_xyz( var_u, zu, ddzu, 1.0_wp, x, y, var_index, j, i )
787	var_index = var_index + 1
788	ENDDO
789	ENDIF
790
791	ENDDO
792	!
793	!-- Write local data on global array.
794	#if defined( __parallel )
795	CALL MPI_ALLREDUCE( sensor_l(1,1), sensor(1,1), num_var_fl * num_leg, MPI_REAL, MPI_SUM, &
796	comm2d, ierr )
797	#else
798	sensor = sensor_l
799	#endif
800	ENDIF
801
802	CALL cpu_log( log_point(65), 'flight_measurement', 'stop' )
803
804	END SUBROUTINE flight_measurement
805
806	!--------------------------------------------------------------------------------------------------!
807	! Description:
808	! ------------
809	!> This subroutine bi-linearly interpolates the respective data onto the current flight position.
810	!--------------------------------------------------------------------------------------------------!
811	SUBROUTINE interpolate_xyz( var, z_uw, ddz_uw, fac, x, y, var_ind, j, i )
812
813	USE control_parameters, &
814	ONLY: dz, &
815	dz_stretch_level_start
816
817	USE grid_variables, &
818	ONLY: dx, &
819	dy
820
821	USE indices, &
822	ONLY: nzb, &
823	nzt, &
824	nxlg, &
825	nxrg, &
826	nysg, &
827	nyng
828
829	IMPLICIT NONE
830
831	INTEGER(iwp) :: i !< index along x
832	INTEGER(iwp) :: j !< index along y
833	INTEGER(iwp) :: k !< index along z
834	INTEGER(iwp) :: k1 !< dummy variable
835	INTEGER(iwp) :: var_ind !< index variable for output quantity
836
837	REAL(wp) :: aa !< dummy argument for horizontal interpolation
838	REAL(wp) :: bb !< dummy argument for horizontal interpolation
839	REAL(wp) :: cc !< dummy argument for horizontal interpolation
840	REAL(wp) :: dd !< dummy argument for horizontal interpolation
841	REAL(wp) :: gg !< dummy argument for horizontal interpolation
842	REAL(wp) :: fac !< flag to indentify if quantity is on zu or zw level
843	REAL(wp) :: var_int !< horizontal interpolated variable at current position
844	REAL(wp) :: var_int_l !< horizontal interpolated variable at k-level
845	REAL(wp) :: var_int_u !< horizontal interpolated variable at (k+1)-level
846	REAL(wp) :: x !< distance between left edge of current grid box and flight position
847	REAL(wp) :: y !< distance between south edge of current grid box and flight position
848
849	REAL(wp), DIMENSION(1:nzt+1) :: ddz_uw !< inverse vertical grid spacing
850	REAL(wp), DIMENSION(nzb:nzt+1) :: z_uw !< height level
851
852	REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) :: var !< treated quantity
853	!
854	!-- Calculate interpolation coefficients
855	aa = x2 + y2
856	bb = ( dx - x )2 + y2
857	cc = x2 + ( dy - y )2
858	dd = ( dx - x )2 + ( dy - y )2
859	gg = aa + bb + cc + dd
860	!
861	!-- Obtain vertical index. Special treatment for grid index along z-direction if flight position is
862	!-- above the vertical grid-stretching level. fac=1 if variable is on scalar grid level, fac=0 for
863	!-- w-component.
864	IF ( z_pos(l) < dz_stretch_level_start(1) ) THEN
865	k = ( z_pos(l) + fac * 0.5_wp * dz(1) ) / dz(1)
866	ELSE
867	!
868	!-- Search for k-index
869	DO k1 = nzb, nzt
870	IF ( z_pos(l) >= z_uw(k1) .AND. z_pos(l) < z_uw(k1+1) ) THEN
871	k = k1
872	EXIT
873	ENDIF
874	ENDDO
875	ENDIF
876	!
877	!-- (x,y)-interpolation at k-level
878	var_int_l = ( ( gg - aa ) * var(k,j,i) + &
879	( gg - bb ) * var(k,j,i+1) + &
880	( gg - cc ) * var(k,j+1,i) + &
881	( gg - dd ) * var(k,j+1,i+1) &
882	) / ( 3.0_wp * gg )
883	!
884	!-- (x,y)-interpolation on (k+1)-level
885	var_int_u = ( ( gg - aa ) * var(k+1,j,i) + &
886	( gg - bb ) * var(k+1,j,i+1) + &
887	( gg - cc ) * var(k+1,j+1,i) + &
888	( gg - dd ) * var(k+1,j+1,i+1) &
889	) / ( 3.0_wp * gg )
890	!
891	!-- z-interpolation onto current flight postion
892	var_int = var_int_l + ( z_pos(l) - z_uw(k) ) * ddz_uw(k+1) * (var_int_u - var_int_l )
893	!
894	!-- Store on local data array
895	sensor_l(var_ind,l) = var_int
896
897	END SUBROUTINE interpolate_xyz
898
899	!--------------------------------------------------------------------------------------------------!
900	! Description:
901	! ------------
902	!> Perform parameter checks.
903	!--------------------------------------------------------------------------------------------------!
904	SUBROUTINE flight_check_parameters
905
906	USE arrays_3d, &
907	ONLY: zu
908
909	USE control_parameters, &
910	ONLY: bc_lr_cyc, &
911	bc_ns_cyc, &
912	message_string
913
914	USE grid_variables, &
915	ONLY: dx, &
916	dy
917
918	USE indices, &
919	ONLY: nx, &
920	ny, &
921	nz
922
923	USE netcdf_interface, &
924	ONLY: netcdf_data_format
925
926	IMPLICIT NONE
927
928	!
929	!-- Check if start positions are properly set.
930	DO l = 1, num_leg
931	IF ( x_start(l) < 0.0_wp .OR. x_start(l) > ( nx + 1 ) * dx ) THEN
932	message_string = 'Start x position is outside the model domain'
933	CALL message( 'flight_check_parameters', 'PA0431', 1, 2, 0, 6, 0 )
934	ENDIF
935	IF ( y_start(l) < 0.0_wp .OR. y_start(l) > ( ny + 1 ) * dy ) THEN
936	message_string = 'Start y position is outside the model domain'
937	CALL message( 'flight_check_parameters', 'PA0432', 1, 2, 0, 6, 0 )
938	ENDIF
939
940	ENDDO
941	!
942	!-- Check for leg mode
943	DO l = 1, num_leg
944	!
945	!-- Check if leg mode matches the overall lateral model boundary conditions.
946	IF ( TRIM( leg_mode(l) ) == 'cyclic' ) THEN
947	IF ( .NOT. bc_lr_cyc .OR. .NOT. bc_ns_cyc ) THEN
948	message_string = 'Cyclic flight leg does not match lateral boundary condition'
949	CALL message( 'flight_check_parameters', 'PA0433', 1, 2, 0, 6, 0 )
950	ENDIF
951	!
952	!-- Check if end-positions are inside the model domain in case of return-legs.
953	ELSEIF ( TRIM( leg_mode(l) ) == 'return' ) THEN
954	IF ( x_end(l) > ( nx + 1 ) * dx .OR. y_end(l) > ( ny + 1 ) * dx ) THEN
955	message_string = 'Flight leg or parts of it are outside the model domain'
956	CALL message( 'flight_check_parameters', 'PA0434', 1, 2, 0, 6, 0 )
957	ENDIF
958	ELSE
959	message_string = 'Unknown flight mode'
960	CALL message( 'flight_check_parameters', 'PA0435', 1, 2, 0, 6, 0 )
961	ENDIF
962
963	ENDDO
964	!
965	!-- Check if given flight object remains inside model domain if a rate of climb / descent is
966	!-- prescribed.
967	DO l = 1, num_leg
968	IF ( flight_level(l) + max_elev_change(l) > zu(nz) .OR. &
969	flight_level(l) + max_elev_change(l) <= 0.0_wp ) THEN
970	message_string = 'Flight level is outside the model domain '
971	CALL message( 'flight_check_parameters', 'PA0438', 1, 2, 0, 6, 0 )
972	ENDIF
973	ENDDO
974	!
975	!-- Check for appropriate NetCDF format. Definition of more than one unlimited dimension is
976	!-- unfortunately only possible with NetCDF4/HDF5.
977	IF ( netcdf_data_format <= 2 ) THEN
978	message_string = 'netcdf_data_format must be > 2'
979	CALL message( 'flight_check_parameters', 'PA0439', 1, 2, 0, 6, 0 )
980	ENDIF
981
982
983	END SUBROUTINE flight_check_parameters
984
985
986	!--------------------------------------------------------------------------------------------------!
987	! Description:
988	! ------------
989	!> Read module-specific global restart data (Fortran binary format).
990	!--------------------------------------------------------------------------------------------------!
991	SUBROUTINE flight_rrd_global_ftn( found )
992
993
994	USE control_parameters, &
995	ONLY: length, &
996	restart_string
997
998
999	IMPLICIT NONE
1000
1001	LOGICAL, INTENT(OUT) :: found !< flag indicating if a variable string is found
1002
1003
1004	found = .TRUE.
1005
1006
1007	SELECT CASE ( restart_string(1:length) )
1008
1009	CASE ( 'u_agl' )
1010	IF ( .NOT. ALLOCATED( u_agl ) ) ALLOCATE( u_agl(1:num_leg) )
1011	READ ( 13 ) u_agl
1012	CASE ( 'v_agl' )
1013	IF ( .NOT. ALLOCATED( v_agl ) ) ALLOCATE( v_agl(1:num_leg) )
1014	READ ( 13 ) v_agl
1015	CASE ( 'w_agl' )
1016	IF ( .NOT. ALLOCATED( w_agl ) ) ALLOCATE( w_agl(1:num_leg) )
1017	READ ( 13 ) w_agl
1018	CASE ( 'x_pos' )
1019	IF ( .NOT. ALLOCATED( x_pos ) ) ALLOCATE( x_pos(1:num_leg) )
1020	READ ( 13 ) x_pos
1021	CASE ( 'y_pos' )
1022	IF ( .NOT. ALLOCATED( y_pos ) ) ALLOCATE( y_pos(1:num_leg) )
1023	READ ( 13 ) y_pos
1024	CASE ( 'z_pos' )
1025	IF ( .NOT. ALLOCATED( z_pos ) ) ALLOCATE( z_pos(1:num_leg) )
1026	READ ( 13 ) z_pos
1027
1028	CASE DEFAULT
1029
1030	found = .FALSE.
1031
1032	END SELECT
1033
1034
1035	END SUBROUTINE flight_rrd_global_ftn
1036
1037
1038
1039	!------------------------------------------------------------------------------!
1040	! Description:
1041	! ------------
1042	!> Read module-specific global restart data (MPI-IO).
1043	!------------------------------------------------------------------------------!
1044	SUBROUTINE flight_rrd_global_mpi
1045
1046
1047	IMPLICIT NONE
1048
1049	LOGICAL :: array_found !< flag indicating if respective array is found in restart file
1050
1051
1052	CALL rd_mpi_io_check_array( 'u_agl', found = array_found )
1053	IF ( array_found) THEN
1054	IF ( .NOT. ALLOCATED( u_agl ) ) ALLOCATE( u_agl(1:num_leg) )
1055	CALL rrd_mpi_io_global_array( 'u_agl', u_agl )
1056	ENDIF
1057	CALL rd_mpi_io_check_array( 'v_agl', found = array_found )
1058	IF ( array_found) THEN
1059	IF ( .NOT. ALLOCATED( v_agl ) ) ALLOCATE( v_agl(1:num_leg) )
1060	CALL rrd_mpi_io_global_array( 'v_agl', v_agl )
1061	ENDIF
1062	CALL rd_mpi_io_check_array( 'w_agl', found = array_found )
1063	IF ( array_found) THEN
1064	IF ( .NOT. ALLOCATED( w_agl ) ) ALLOCATE( w_agl(1:num_leg) )
1065	CALL rrd_mpi_io_global_array( 'w_agl', w_agl )
1066	ENDIF
1067	CALL rd_mpi_io_check_array( 'x_pos', found = array_found )
1068	IF ( array_found) THEN
1069	IF ( .NOT. ALLOCATED( x_pos ) ) ALLOCATE( x_pos(1:num_leg) )
1070	CALL rrd_mpi_io_global_array( 'x_pos', x_pos )
1071	ENDIF
1072	CALL rd_mpi_io_check_array( 'y_pos', found = array_found )
1073	IF ( array_found) THEN
1074	IF ( .NOT. ALLOCATED( y_pos ) ) ALLOCATE( y_pos(1:num_leg) )
1075	CALL rrd_mpi_io_global_array( 'y_pos', y_pos )
1076	ENDIF
1077	CALL rd_mpi_io_check_array( 'z_pos', found = array_found )
1078	IF ( array_found) THEN
1079	IF ( .NOT. ALLOCATED( z_pos ) ) ALLOCATE( z_pos(1:num_leg) )
1080	CALL rrd_mpi_io_global_array( 'z_pos', z_pos )
1081	ENDIF
1082
1083	END SUBROUTINE flight_rrd_global_mpi
1084
1085
1086	!--------------------------------------------------------------------------------------------------!
1087	! Description:
1088	! ------------
1089	!> This routine writes the respective restart data.
1090	!--------------------------------------------------------------------------------------------------!
1091	SUBROUTINE flight_wrd_global
1092
1093
1094	IMPLICIT NONE
1095
1096
1097	IF ( TRIM( restart_data_format_output ) == 'fortran_binary' ) THEN
1098
1099	CALL wrd_write_string( 'u_agl' )
1100	WRITE ( 14 ) u_agl
1101
1102	CALL wrd_write_string( 'v_agl' )
1103	WRITE ( 14 ) v_agl
1104
1105	CALL wrd_write_string( 'w_agl' )
1106	WRITE ( 14 ) w_agl
1107
1108	CALL wrd_write_string( 'x_pos' )
1109	WRITE ( 14 ) x_pos
1110
1111	CALL wrd_write_string( 'y_pos' )
1112	WRITE ( 14 ) y_pos
1113
1114	CALL wrd_write_string( 'z_pos' )
1115	WRITE ( 14 ) z_pos
1116
1117	ELSEIF ( restart_data_format_output(1:3) == 'mpi' ) THEN
1118
1119	CALL wrd_mpi_io_global_array( 'u_agl', u_agl )
1120	CALL wrd_mpi_io_global_array( 'v_agl', v_agl )
1121	CALL wrd_mpi_io_global_array( 'w_agl', w_agl )
1122	CALL wrd_mpi_io_global_array( 'x_pos', x_pos )
1123	CALL wrd_mpi_io_global_array( 'y_pos', y_pos )
1124	CALL wrd_mpi_io_global_array( 'z_pos', z_pos )
1125
1126	ENDIF
1127
1128	END SUBROUTINE flight_wrd_global
1129
1130
1131	END MODULE flight_mod

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |