Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |