Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |