Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |