Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |