Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |