program run_viper implicit none integer, parameter :: mx=2000, nzmax=401 real y_meas_p(nzmax), z_meas_p(nzmax), cir_meas_p(nzmax), & y_meas_s(nzmax), z_meas_s(nzmax), cir_meas_s(nzmax), & time_meas_p( nzmax), time_meas_s( nzmax), & time_meas_cir_p(nzmax), time_meas_cir_s(nzmax) integer nt_meas_p, nt_meas_s, nt_meas_cir_p, nt_meas_cir_s real time_model(mx), & y_model_p(mx), z_model_p(mx), cir_model_p(mx), & y_model_s(mx), z_model_s(mx), cir_model_s(mx), & time_tmp(mx), hw_tmp_p(mx), hw_tmp_s(mx), & x_tmp_p(mx), y_tmp_p(mx), z_tmp_p(mx), cir_tmp_p(mx), & x_tmp_s(mx), y_tmp_s(mx), z_tmp_s(mx), cir_tmp_s(mx), & work(nzmax*12), cw_orig(nzmax), cw_spln(nzmax) real ac_weight, ac_span, ac_speed, b0, R0, y0, z0, & heading, runway_elev, surf_wind_dir, surf_wind_spd, & surf_temp, surf_rho, msl_pressure, stability, turbulence, & v0, Gamma_oo, rho_z0, hw_avg integer units_flag, flag_cw, flag_hw, flag_q, flag_sim, n_steps, & i_unstable, ierr, nt_model, i_hw real t1, t2, dt, time_max, data_time_max, default_time_max, & max_t_valu, max_t_step, t_scale, tld, tssd, zmf, maxsub real pi, tan3, time_now, time_release, z_now, z_est, & x_release, z_release, z_release_sav, z_release_pre, & z_rel_est, hw_z0, hw_est, hw_release, hw_rel_est, w1, w2, & hw_vort_avg, hw_ac_avg, hw_sav, hw_err, hw_add, hw_interp, & z0_tmp, t1_tmp, t2_tmp, del_t, deviation, & x_err, mu, mu_ref, T_ref, T_s, dev_cw integer L, khw, khw_sav, nt_tmp, n_steps_tmp integer i, ii, j, n, i_beg, i_end, & i_base_beg, i_base_end, len_base, i_ext_beg, len_arg logical no_prof_t, no_prof_hw, no_prof_cw, no_prof_q, & write_lin, write_star, model_unknown, input_unknown character( 5) ac_type, output_extension character( 7) runway character( 10) max_t_mthd character( 12) model character( 20) ac_descriptor, extension character( 80) base_name character(200) arg, fname real grav, tzero, adlaps, c1, c2, c3, zz, zim, zge, yover, zdown, & gamma_gas, R_gas, nu common /cvalb/grav,tzero,adlaps,c1,c2,c3,zz,zim,zge,yover,zdown, & gamma_gas, R_gas, nu real tempz, temp, cwz, cw, hwz, hw, qz, q, & rho, theta, press, bvfsq integer ntemp, nq, ncw, nhw common /atmb/tempz(nzmax), temp(nzmax), qz(nzmax), q(nzmax), & cwz(nzmax), cw(nzmax), hwz(nzmax), hw(nzmax), & press(nzmax), theta(nzmax), bvfsq(nzmax), rho(nzmax), & ntemp, nq, ncw, nhw grav = 9.81 gamma_gas = 1.4 R_gas = 287.0 pi = acos(-1.0) tan3 = tan(3.0/180.0*pi) write_lin = .false. write_star = .false. model_unknown = .true. input_unknown = .true. c *** Parse the command line inputs to get the model name and input file. n = 1 do i=1, 20 call get_command_argument( n, arg, len_arg ) if( len_arg .eq. 0 ) goto 10 select case( arg(1:2) ) case( '-m' ) n = n + 1 call get_command_argument( n, arg, len_arg ) if( len_arg .ne. 0 ) then model = trim(arg) model_unknown = .false. end if case( '-i' ) n = n + 1 call get_command_argument( n, arg, len_arg ) if( len_arg .ne. 0 ) then fname = trim(arg) input_unknown = .false. end if end select n = n + 1 end do 10 continue c *** Set i_unstable=0 to convert unstable regions to neutral i_unstable = 0 c *** Ask the user to input the model name and/or input file if they were c *** not supplied on the command line. if( model_unknown ) then print *, 'Enter the name of the model to run (e.g viper20)' read(5,'(a12)') model end if if( input_unknown ) then print *, 'Enter the name of the input file' read(5,'(a200)') fname end if c *** Set the output file extension based on the selected model. select case( model ) case('viper30' ) output_extension = 'vpr30' case('viper20' ) output_extension = 'vpr20' case default print *, trim(model), ' is not a valid model' stop end select c *** Parse the input file path to extract the case name. i_beg = verify(fname,' ') i_end = verify(fname,' ',.true.) i_base_beg = scan( fname, '/\', .true. ) + 1 i_base_beg = max( i_base_beg, i_beg ) i_base_end = scan( fname, '.', .true. ) - 1 if( i_base_end .eq. -1 ) i_base_end = i_end len_base = i_base_end - i_base_beg + 1 base_name = fname(i_base_beg:i_base_end) print *, ' ' print *, 'Running ', trim(base_name), ' with ', trim(model) c *** Read the input file call read_inputs(fname(i_beg:i_end), & units_flag, ac_type, ac_descriptor, ac_weight, ac_span, & ac_speed, b0, R0, y0, z0, runway, heading, & runway_elev, surf_wind_dir, surf_wind_spd, surf_temp, & surf_rho, msl_pressure, stability, turbulence, & v0, Gamma_oo, rho_z0, flag_cw, flag_hw, flag_q, flag_sim, & i_unstable, tempz, temp, theta, ntemp, hwz, hw, nhw, & cwz, cw, ncw, qz, q, nq, rho, press, bvfsq, hw_avg, & time_meas_p, y_meas_p, z_meas_p, nt_meas_p, & time_meas_s, y_meas_s, z_meas_s, nt_meas_s, & time_meas_cir_p, cir_meas_p, nt_meas_cir_p, & time_meas_cir_s, cir_meas_s, nt_meas_cir_s, nzmax, & grav, R_gas, gamma_gas, & no_prof_t, no_prof_hw, no_prof_cw, no_prof_q, ierr ) if( ierr .ne. 0 ) stop cw_orig = cw t_scale = b0/v0 c *** Read the preferences data call read_prefs( max_t_mthd, max_t_valu, max_t_step ) c *** Determine the maximum time for the vortex measurements if this c *** is to be used as the maximum model time. default_time_max = 240.0 if( trim(max_t_mthd) .eq. 'measured' ) then data_time_max = max( maxsub(time_meas_p,nt_meas_p), & maxsub(time_meas_s,nt_meas_s), & maxsub(time_meas_cir_p,nt_meas_cir_p), & maxsub(time_meas_cir_s,nt_meas_cir_s), & 0.0 ) if( data_time_max .gt. 0.0 ) then time_max = ( int(data_time_max/max_t_step) + 1 )*max_t_step else time_max = default_time_max end if else if( trim(max_t_mthd) .eq. 'model' ) then time_max = default_time_max else time_max = max_t_valu end if c *** Set the model parameter values call set_vpr20_params( base_name, zmf ) call set_vpr30_params( zmf ) t1 = 0.0 t2 = time_max dt = 0.10*t_scale n_steps = nint((t2-t1)/dt) zim = zmf*b0 zz = z0 c *** Smooth and spline fit the crosswind data. if( trim(model) .eq. 'viper30' .and. ncw .ge. 3 ) then dev_cw = 0.1 call auto_smooth( cwz, cw, work, ncw, 1, dev_cw, 0.05 ) do i=1, ncw cw(i) = work(i) cw_spln(i) = work(2*ncw+i) end do end if c *** Compute the kinematic viscosity at the surface (needed for viper_ige). mu_ref = 1.7885e-5 T_ref = 288.15 T_s = 110.0 mu = mu_ref*(surf_temp/T_ref)**1.5*(T_ref+T_s)/(surf_temp+T_s) nu = mu/surf_rho c *** Check to see if we need to use the exact procedure for headwind c *** effects. To do this, we estimate the altitude deviation from z0 c *** where the vortices must be released in order to be advected to the c *** measurement plane at time t2. time_release = hw_avg/ac_speed*t2 x_release = ( ac_speed - hw_avg )*time_release z_release = z0 - x_release*tan3 deviation = abs(z0-z_release)/z0 i_hw = 1 if( deviation .lt. 5.0e-2 ) then c *** The headwind is either small or zero. We can use the simple c *** built-in headwind correction scheme. call execute_model( model, & y0, z0, v0, b0, R0, Gamma_oo, ac_speed, & runway_elev, surf_rho, i_hw, hw_avg, & ntemp, tempz, temp, theta, rho, bvfsq, & ncw, cwz, cw, cw_spln, nhw, hwz, hw, nq, qz, q, & t1, t2, n_steps, & nt_model, time_model, & y_model_p, z_model_p, cir_model_p, & y_model_s, z_model_s, cir_model_s, & tld, tssd ) else c *** The headwinds are strong enough where we need to use the exact c *** headwind scheme. Initialize some stuff here. i_hw = 0 dt = ( t2 - t1 )/float(n_steps) khw = nhw - 1 hw_z0 = hw_interp( hw, hwz, nhw, khw, z0 ) z_release = z0 c *** The extrapolation scheme used below requires the release point c *** on the previous time step. Estimate this here by taking a negative c *** time step away from t=0. time_release = -(hw_z0/ac_speed)*dt z_release_pre = z0 - (ac_speed-hw_z0)*time_release*tan3 c *** We need to estimate the average headwind computed over the c *** vortex trajectory for the first time step. Compute this by c *** taking an average of the headwind at the estimated release c *** point and at the estimated vortex location at the measurement plane. z_est = z0 - & ( ((ac_speed-hw_z0)*tan3-v0)*hw_z0/ac_speed + v0 )*dt hw_est = hw_interp( hw, hwz, nhw, khw, z_est ) z_rel_est = 2.0*z_release - z_release_pre hw_rel_est = hw_interp( hw, hwz, nhw, khw, z_rel_est ) hw_vort_avg = 0.5*( hw_rel_est + hw_est ) hw_add = 0.0 c *** Time advancement loop. The effects of a headwind are accounted for c *** explicitly by computing a complete vortex trajectory with unique c *** starting location (release point) for each output time step. c *** The vortex x-coordinate time histories are computed fairly c *** accurately by integrating dx/dt=-hw over the vortex trajectory. c *** For each time step, we compute the length of time to run the model c *** so that the vortices are advected to the measurement plane c *** at the end of the run. An error correction feature is included c *** and, with it, the vortices end up very close to the measurement c *** plane. Vortex tilt in conjuction with a headwind c *** complicates things since the port and starboard vortices will c *** reach the measurement plane at different times (since they experience c *** different headwind histories owing to their different altitude c *** histories). Right now we handle this by choosing a time where c *** the vortices are equi-distant to the measurement plane (one ahead c *** and one behind). do i=2, n_steps+1 time_now = float(i-1)*dt c *** Estimate the vortex altitude when it is released from the c *** aircraft (z_release) via extrapolation. z_release_sav = z_release z_release = 2.0*z_release - z_release_pre z_release_pre = z_release_sav c *** Find the average headwind over the estimated aircraft altitude range. c *** The average headwind over the vortex altitude range will be taken c *** from the previous model run. hw_release = hw_interp( hw, hwz, nhw, khw, z_release ) hw_ac_avg = 0.5*( hw_release + hw_z0 ) c *** Now compute the time and altitude for the vortex release more c *** accurately using the estimated average headwinds. time_release = hw_vort_avg/(ac_speed+hw_vort_avg-hw_ac_avg)* & time_now x_release = ( ac_speed - hw_ac_avg )*time_release z_release = z0 - x_release*tan3 z0_tmp = z_release t2_tmp = time_now t1_tmp = time_release n_steps_tmp = 5*(i-1) call execute_model( model, & y0, z0_tmp, v0, b0, R0, gamma_oo, ac_speed, & runway_elev, surf_rho, i_hw, hw_avg, & ntemp, tempz, temp, theta, rho, bvfsq, & ncw, cwz, cw, cw_spln, nhw, hwz, hw, nq, qz, q, & t1_tmp, t2_tmp, n_steps_tmp, & nt_tmp, time_tmp, & y_tmp_p, z_tmp_p, cir_tmp_p, & y_tmp_s, z_tmp_s, cir_tmp_s, & tld, tssd ) z_now = 0.5*( z_tmp_p(nt_tmp ) + z_tmp_s(nt_tmp ) ) c *** Determine the x position of the vortices as a function of time. khw = nhw-1 do j=1, nt_tmp hw_tmp_p(j) = hw_interp( hw, hwz, nhw, khw, z_tmp_p(j) ) end do khw = nhw-1 do j=1, nt_tmp hw_tmp_s(j) = hw_interp( hw, hwz, nhw, khw, z_tmp_s(j) ) end do x_tmp_p(1) = x_release x_tmp_s(1) = x_release do j=2, nt_tmp del_t = time_tmp(j) - time_tmp(j-1) x_tmp_p(j) = x_tmp_p(j-1) - & 0.5*( hw_tmp_p(j) + hw_tmp_p(j-1) )*del_t x_tmp_s(j) = x_tmp_s(j-1) - & 0.5*( hw_tmp_s(j) + hw_tmp_s(j-1) )*del_t end do c *** Compute the effective average headwind over the vortex trajectory. c *** We will use this value to estimate the length of time to run the c *** model for the next time step. We also apply a bit of error c *** correction since we have access to the estimated and exact c *** headwinds for the current time step at this point. The error c *** correction seems to work well but we should keep an eye on it. hw_sav = hw_vort_avg hw_vort_avg = ( x_release - & 0.5*( x_tmp_p(nt_tmp) + x_tmp_s(nt_tmp) ) )/ & ( t2_tmp - t1_tmp ) hw_err = hw_vort_avg - hw_sav hw_add = hw_add + hw_err hw_vort_avg = hw_vort_avg + hw_add c *** Abort the run if the model did not stop at the expected time if( abs(time_tmp(nt_tmp)-time_now) .gt. 1.0e-4 ) goto 90 c *** The initial condition is contained in the solution for the c *** first run (at i=2) if( i .eq. 2 ) then time_model( 1) = time_tmp( 1) y_model_p( 1) = y_tmp_p( 1) y_model_s( 1) = y_tmp_s( 1) z_model_p( 1) = z_tmp_p( 1) z_model_s( 1) = z_tmp_s( 1) cir_model_p(1) = cir_tmp_p(1) cir_model_s(1) = cir_tmp_s(1) end if c *** The last point in the solution arrays should have the vortices at c *** the measurement plane. Save just these values to the global solution c *** arrays. time_model( i) = time_tmp( nt_tmp) y_model_p( i) = y_tmp_p( nt_tmp) y_model_s( i) = y_tmp_s( nt_tmp) z_model_p( i) = z_tmp_p( nt_tmp) z_model_s( i) = z_tmp_s( nt_tmp) cir_model_p(i) = cir_tmp_p(nt_tmp) cir_model_s(i) = cir_tmp_s(nt_tmp) end do 90 nt_model = i - 1 end if c *** Write the measured data for plotting purposes. call write_inputs( base_name, nzmax, mx, c & tempz, temp, cwz, cw_orig, hwz, hw, qz, q, & tempz, temp, cwz, cw, hwz, hw, qz, q, & rho, press, theta, bvfsq, & time_meas_p, y_meas_p, z_meas_p, & time_meas_s, y_meas_s, z_meas_s, & time_meas_cir_p, cir_meas_p, & time_meas_cir_s, cir_meas_s, & time_model, & y_model_p, z_model_p, cir_model_p, & y_model_s, z_model_s, cir_model_s, & y0, z0, v0, b0, Gamma_oo, zmf, & ac_speed, hw_avg, ntemp, ncw, nhw, nq, & nt_meas_p, nt_meas_s, & nt_meas_cir_p, nt_meas_cir_s, & nt_model, tld, tssd, & write_lin, write_star ) c *** Write the model results open(unit=25, & file='Output/'//trim(base_name)//'.'//output_extension) write(25,125) 125 format('# time y_p z_p cir_p y_s', & ' z_s cir_s') do i=1, nt_model write(25,45) time_model(i), & y_model_p(i), z_model_p(i), cir_model_p(i), & y_model_s(i), z_model_s(i), cir_model_s(i) 45 format(f8.4,1p,6e12.4) end do close(25) print *, ' ' print *, model//' is done' print *, ' ' stop end function hw_interp( hw, hwz, nhw, khw, z ) real hw(nhw), hwz(nhw) if( z .gt. hwz(nhw) ) then print *, 'ERROR in hw_interp, z is out of range. z = ', z stop end if khw = min( khw, nhw-1 ) k = khw if( z .gt. hwz(k) .and. z .le. hwz(k+1) ) goto 5 if( z .gt. hwz(k+1) ) khw = nhw-1 do k=khw, 1, -1 if( hwz(k) .le. z ) goto 5 end do print *, 'ERROR in hw_interp, z is out of range. z = ', z stop 5 khw = k wp = ( z - hwz(khw) )/( hwz(khw+1) - hwz(khw) ) wm = 1.0 - wp hw_interp = wm*hw(khw) + wp*hw(khw+1) return end real function maxsub(a,n) real a(n) maxsub = -1.0e-15 do i=1, n maxsub = max(maxsub,a(i)) end do return end