Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |