Home

Context Navigation

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

Last change on this file since 2000 was 2000, checked in by knoop, 8 years ago
Forced header and separation lines into 80 columns
Property svn:keywords set to `Id`
File size: 37.4 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |