Home

Context Navigation

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

Last change on this file since 3049 was 3049, checked in by Giersch, 6 years ago
Revision history corrected
Property svn:keywords set to `Id`
File size: 39.0 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |