Home

Context Navigation

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

Last change on this file since 1957 was 1957, checked in by suehring, 8 years ago
flight module added
Property svn:keywords set to `Id`
File size: 36.0 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |