Home

Context Navigation

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

Last change on this file since 3371 was 3274, checked in by knoop, 6 years ago
Modularization of all bulk cloud physics code components
Property svn:keywords set to `Id`
File size: 39.8 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |