Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |