Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |