Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |