Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |