Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |