Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |