Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |