Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |