Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |