Home

Context Navigation

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

Last change on this file since 1958 was 1958, checked in by suehring, 8 years ago
last commit documented
Property svn:keywords set to `Id`
File size: 36.0 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |