Home

Context Navigation

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

Last change on this file since 3909 was 3885, checked in by kanani, 6 years ago
restructure/add location/debug messages
Property svn:keywords set to `Id`
File size: 40.3 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |