Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |