Home

Context Navigation

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

Last change on this file since 2825 was 2776, checked in by Giersch, 7 years ago
Skipping of module related restart data changed + adapting synthetic turbulence generator to current restart procedure
Property svn:keywords set to `Id`
File size: 38.7 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |