subroutine read_inputs( fname, & 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, gamma0, rho_z0, flag_cw, flag_hw, flag_q, flag_sim, & i_unstable, tempz, temp, theta, n_temp, hwz, hw, n_hw, & cwz, cw, n_cw, qz, q, n_q, 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 ) implicit none include 'ige.h' real tempz(nzmax), temp(nzmax), theta(nzmax), bvfsq(nzmax), & hwz(nzmax), hw(nzmax), cwz(nzmax), cw(nzmax), & qz(nzmax), q(nzmax), rho(nzmax), press(nzmax), & time_meas_p(nzmax), y_meas_p(nzmax), z_meas_p(nzmax), & time_meas_s(nzmax), y_meas_s(nzmax), z_meas_s(nzmax), & time_meas_cir_p(nzmax), cir_meas_p(nzmax), & time_meas_cir_s(nzmax), cir_meas_s(nzmax), & ac_weight, ac_span, ac_speed, & b0, R0, y0, z0, heading, runway_elev, & surf_wind_dir, surf_wind_spd, surf_temp, & msl_pressure, stability, turbulence, v0, gamma0, & rho_z0, grav, R_gas, gamma_gas, hw_avg integer units_flag, flag_pot_temp, flag_cw, flag_hw, flag_q, & flag_sim, i_unstable, n_temp, n_hw, n_cw, n_q, & nt_meas_p, nt_meas_s, nt_meas_cir_p, & nt_meas_cir_s, nzmax character*5 ac_type character*7 runway character*20 ac_descriptor character*80 base_name character(*) fname logical no_prof_t, no_prof_hw, no_prof_cw, no_prof_q real pi, c13, cp_gas, kt2ms, ft2m, inHg2Pa, lb2kg, & slugft32kgm3, ftsec2msec, beta, & T_standard_0, expon, arg, surf_rho, theta_0, & surf_press, theta_z0, T_z0, c1, e1, p_z0, z, dz, & dTdz, z_top, del_z, del_T, dT, & fac, angle, head_wind, cross_wind, weight, gamma, & gamma1, gamma2, edr, edr_star, tlag_fac,tige1_fac, & ff2, ff3, gg2, gg3, yp0, ys0, zp0, zs0, vyp, vys, & vzp, vzs, cwp0, cws0, vel_factor, sum2, z_min integer i, L, ic, ierr, len_fname, len_label, i_lab_beg, & i_lab_end, i_val_beg, i_val_end, & i_base_beg, i_base_end, isamp & logical found_units_flag, found_ac_weight, found_ac_span, & found_ac_speed, found_b0, found_y0, & found_z0, found_runway, found_heading, & found_runway_elev, found_surf_wind_dir, & found_surf_wind_spd, found_surf_temp, found_surf_rho, & found_msl_pressure, found_stability, & found_turbulence, found_v0, found_gamma0, & found_flag_pot_temp character*3 gamma_bool character*20 label character*80 line pi = acos(-1.0) c13 = 1.0/3.0 cp_gas = gamma_gas/(gamma_gas-1.0)*R_gas kt2ms = 1852.0/3600.0 ft2m = (2.54*12.0)/100.0 inHg2Pa = 101325.0/29.92 lb2kg = 0.453592 slugft32kgm3 = 515.37882 ftsec2msec = ft2m ierr = 0 no_prof_t = .false. no_prof_hw = .false. no_prof_cw = .false. no_prof_q = .false. found_units_flag = .false. found_ac_weight = .false. found_ac_span = .false. found_ac_speed = .false. found_b0 = .false. found_y0 = .false. found_z0 = .false. found_runway = .false. found_heading = .false. found_runway_elev = .false. found_surf_wind_dir = .false. found_surf_wind_spd = .false. found_surf_temp = .false. found_surf_rho = .false. found_msl_pressure = .false. found_stability = .false. found_turbulence = .false. found_v0 = .false. found_gamma0 = .false. found_flag_pot_temp = .false. n_temp=0; n_hw=0; n_cw=0; n_q=0; nt_meas_p=0; nt_meas_s=0 nt_meas_cir_p=0; nt_meas_cir_s=0; flag_pot_temp=0; flag_sim=0 time_meas_p(1)=0.0; time_meas_s(1)=0.0; time_meas_cir_p(1)=0.0; time_meas_cir_s(1)=0.0; open(unit=1,file=fname,status='old') c *** Parse the input file. For each line, identify the label, read c *** the input and assign it to the corresponding variable. do L=1, 20000 c *** Read a line and extract the label. Comments are flagged with # c *** Blank lines are skipped. The string 'eof' can signify the end c *** of file. read(1,'(a80)',end=60) line if( line(1:1) .eq. '#' ) goto 40 i_lab_beg = scan(line,'ABCDEFGHIJKLMNOPQRSTUVWXYZ'// & 'abcdefghijklmnopqrstuvwxyz' ) if( i_lab_beg .eq. 0 ) goto 40 i_lab_end = index(line(i_lab_beg+1:80),' ')-1 + i_lab_beg len_label = i_lab_end-i_lab_beg+1 label = line(i_lab_beg:i_lab_end) c *** Check for non-alpha characters and convert to lower case do i=1, len_label ic = ichar(label(i:i)) if( ic .lt. 48 .or. & (ic .gt. 57 .and. ic .lt. 65) .or. & (ic .gt. 90 .and. ic .lt. 95) .or. & ic .eq. 96 .or. & ic .gt. 122 ) then print *, 'ERROR: while reading file ', fname print *, 'Invalid character in label = ',label ierr = 4 print *, 'line = ',line return end if if( ic .ge. 65 .and. ic .le. 90 ) then label(i:i) = achar(ic+32) end if end do if( label(1:len_label) .eq. 'eof' ) goto 60 c *** Determine the position of the parameter value do i=i_lab_end+1, 80 if( line(i:i) .ne. ' ' ) goto 15 end do 15 i_val_beg = i i_val_end = index(line(i_val_beg+1:80),' ')-1 + i_val_beg c *** Match the label with the corresponding parameter and set its value. select case( label(1:len_label) ) case( 'units_flag' ) read(line(i_val_beg:i_val_end),*,err=50) units_flag found_units_flag = .true. case( 'ac_type' ) ac_type = line(i_val_beg:i_val_end) case( 'ac_descriptor' ) c ac_descriptor = line(i_val_beg:i_val_end) case( 'ac_weight' ) read(line(i_val_beg:i_val_end),*,err=50) ac_weight found_ac_weight = .true. case( 'ac_span' ) read(line(i_val_beg:i_val_end),*,err=50) ac_span found_ac_span = .true. case( 'ac_speed' ) read(line(i_val_beg:i_val_end),*,err=50) ac_speed found_ac_speed = .true. case( 'ac_b0' ) read(line(i_val_beg:i_val_end),*,err=50) b0 found_b0 = .true. case( 'ac_y0' ) read(line(i_val_beg:i_val_end),*,err=50) y0 found_y0 = .true. case( 'ac_z0' ) read(line(i_val_beg:i_val_end),*,err=50) z0 found_z0 = .true. case( 'runway' ) runway = line(i_val_beg:i_val_end) found_runway = .true. case( 'heading' ) read(line(i_val_beg:i_val_end),*,err=50) heading found_heading = .true. case( 'runway_elev' ) read(line(i_val_beg:i_val_end),*,err=50) runway_elev found_runway_elev = .true. case( 'surf_wind_dir' ) read(line(i_val_beg:i_val_end),*,err=50) surf_wind_dir found_surf_wind_dir = .true. case( 'surf_wind_spd' ) read(line(i_val_beg:i_val_end),*,err=50) surf_wind_spd found_surf_wind_spd = .true. case( 'surf_temp' ) read(line(i_val_beg:i_val_end),*,err=50) surf_temp found_surf_temp = .true. case( 'surf_rho' ) read(line(i_val_beg:i_val_end),*,err=50) surf_rho found_surf_rho = .true. case( 'msl_pressure' ) read(line(i_val_beg:i_val_end),*,err=50) msl_pressure found_msl_pressure = .true. case( 'stability' ) read(line(i_val_beg:i_val_end),*,err=50) stability found_stability = .true. case( 'turbulence' ) read(line(i_val_beg:i_val_end),*,err=50) turbulence found_turbulence = .true. case( 'v0' ) read(line(i_val_beg:i_val_end),*,err=50) v0 found_v0 = .true. case( 'gamma0' ) read(line(i_val_beg:i_val_end),*,err=50) gamma0 found_gamma0 = .true. case( 'flag_sim' ) read(line(i_val_beg:i_val_end),*,err=50) flag_sim case( 'flag_pot_temp' ) read(line(i_val_beg:i_val_end),*,err=50) flag_pot_temp found_flag_pot_temp = .true. case( 'n_temp' ) read(line(i_val_beg:i_val_end),*,err=50) n_temp do i=1, n_temp read(1,*) tempz(i), temp(i) end do case( 'flag_hw' ) read(line(i_val_beg:i_val_end),*,err=50) flag_hw case( 'n_hw' ) read(line(i_val_beg:i_val_end),*,err=50) n_hw do i=1, n_hw read(1,*) hwz(i), hw(i) end do case( 'flag_cw' ) read(line(i_val_beg:i_val_end),*,err=50) flag_cw case( 'n_cw' ) read(line(i_val_beg:i_val_end),*,err=50) n_cw do i=1, n_cw read(1,*) cwz(i), cw(i) end do case( 'flag_q' ) read(line(i_val_beg:i_val_end),*,err=50) flag_q case( 'n_q' ) read(line(i_val_beg:i_val_end),*,err=50) n_q do i=1, n_q read(1,*) qz(i), q(i) end do case( 'nt_meas_p' ) read(line(i_val_beg:i_val_end),*,err=50) nt_meas_p do i=1, nt_meas_p read(1,*) time_meas_p(i), y_meas_p(i), z_meas_p(i) end do case( 'nt_meas_s' ) read(line(i_val_beg:i_val_end),*,err=50) nt_meas_s do i=1, nt_meas_s read(1,*) time_meas_s(i), y_meas_s(i), z_meas_s(i) end do case( 'nt_meas_cir_p' ) read(line(i_val_beg:i_val_end),*,err=50) nt_meas_cir_p do i=1, nt_meas_cir_p read(1,*) time_meas_cir_p(i), cir_meas_p(i) end do case( 'nt_meas_cir_s' ) read(line(i_val_beg:i_val_end),*,err=50) nt_meas_cir_s do i=1, nt_meas_cir_s read(1,*) time_meas_cir_s(i), cir_meas_s(i) end do case default print *, 'ERROR: Unrecognized label = ', label ierr = 4 return end select 40 continue end do 50 ierr = 2 print *, 'ERROR: Malformed numeric entry for ', label(1:len_label) return 60 continue close(1) c *** Convert all quantities to SI units if( units_flag .eq. 1 ) then if( found_ac_weight ) then ac_weight = ac_weight *lb2kg end if if( found_ac_span ) then ac_span = ac_span *ft2m end if if( found_ac_speed ) then ac_speed = ac_speed *kt2ms end if if( found_b0 ) then b0 = b0 *ft2m end if if( found_y0 ) then y0 = y0 *ft2m end if if( found_z0 ) then z0 = z0 *ft2m end if if( found_msl_pressure ) then msl_pressure = msl_pressure *inHg2Pa end if if( found_runway_elev ) then runway_elev = runway_elev *ft2m end if if( found_surf_rho ) then surf_rho = surf_rho *slugft32kgm3 end if if( found_surf_wind_spd ) then surf_wind_spd = surf_wind_spd *kt2ms end if if( found_turbulence ) then turbulence = turbulence *ft2m**2 end if if( found_v0 ) then v0 = v0 *ftsec2msec end if if( found_gamma0 ) then gamma0 = gamma0 *ft2m**2 end if do i=1, n_temp tempz(i) = tempz(i)*ft2m end do do i=1, n_hw hwz(i) = hwz(i)*ft2m hw(i) = hw(i)*kt2ms end do do i=1, n_cw cwz(i) = cwz(i)*ft2m cw(i) = cw(i)*kt2ms end do do i=1, n_q qz(i) = qz(i)*ft2m q(i) = q(i)*ft2m**2 end do do i=1, nt_meas_p y_meas_p(i) = y_meas_p(i)*ft2m z_meas_p(i) = z_meas_p(i)*ft2m end do do i=1, nt_meas_s y_meas_s(i) = y_meas_s(i)*ft2m z_meas_s(i) = z_meas_s(i)*ft2m end do do i=1, nt_meas_cir_p cir_meas_p(i) = cir_meas_p(i)*ft2m**2 end do do i=1, nt_meas_cir_s cir_meas_s(i) = cir_meas_s(i)*ft2m**2 end do end if if( found_surf_temp .and. surf_temp .lt. 150.0 ) then surf_temp = surf_temp + 273.15 end if if( n_temp .gt. 0 .and. temp(1) .lt. 150.0 ) then do i=1, n_temp temp(i) = temp(i) + 273.15 end do end if if( n_temp .gt. 0 ) surf_temp = temp(1) c *** Check for missing required inputs if( .not. found_units_flag ) then print *, 'ERROR: the required input units_flag is missing' stop end if if( .not. found_ac_speed ) then print *, 'ERROR: the required input ac_speed is missing' stop end if if( .not. found_z0 ) then print *, 'ERROR: the required input ac_z0 is missing' stop end if if( n_temp .gt. 0 .and. (.not. found_flag_pot_temp) ) then print *, 'ERROR: you must supply flag_pot_temp for a ', & 'specified temperature profile' stop end if if( .not. found_b0 ) then if( .not. found_ac_span ) then print *, 'ERROR: Neither ac_b0 nor ac_span are given' stop else print *, 'Warning: ac_b0 is not given, estimating this ', & 'from ac_span' b0 = ac_span*pi/4.0 end if end if if( .not. found_ac_span ) ac_span = b0*4.0/pi if( .not. found_y0 ) then print *, 'Warning: ac_y0 is not given.' print *, ' Assuming ac_y0 = 0' end if c *** Define the upper limit for profiles that we create (or alter). c *** Here we choose the first multiple of 10 greater than z0 z_top = float( ( int(z0*0.1) + 1 )*10 ) c *** If the read-in temperature profile contains only two points, c *** interpolate it to more. if( n_temp .eq. 2 ) then if( abs(temp(2)-temp(1)) .lt. 1.0e-4 ) no_prof_t = .true. dTdz = ( temp(2) - temp(1) )/( tempz(2) - tempz(1) ) del_z = z_top - tempz(1) del_T = del_z*dTdz n_temp = 51 dz = del_z/float(n_temp-1) dT = del_T/float(n_temp-1) do i=2, n_temp tempz(i) = tempz(1) + float(i-1)*dz temp(i) = temp(1) + float(i-1)*dT end do end if c *** If we read temperature and flag_pot_temp=1, then we actually c *** read potential temperature. Copy the values into the theta array c *** in this case. if( n_temp .gt. 0 .and. flag_pot_temp .eq. 1 ) then do i=1, n_temp theta(i) = temp(i) end do end if c *** Construct the potential temperature profile if no temperature profile c *** data was given if( n_temp .eq. 0 ) then no_prof_t = .true. if( .not. found_stability ) then print *, 'Warning: attempting to compute temperature ', & 'profile but stability is not given' print *, ' Assuming stability = 0.0' stability = 0.0 end if if( .not. found_surf_temp ) then print *, 'Warning: attempting to compute temperature ', & 'profile but surf_temp is not given' print *, ' Assuming surf_temp = 15.0 C' surf_temp = 15.0 + 273.15 found_surf_temp = .true. end if theta_0 = surf_temp n_temp = 51 dz = z_top/float(n_temp-1) fac = stability**2/grav do i=1, n_temp z = float(i-1)*dz tempz(i) = z theta(i) = theta_0*exp(fac*z) end do flag_pot_temp = 1 end if c *** Make sure that the first profile point is at the surface if( abs(tempz(1)) .gt. 1.0e-3 ) then print *, 'The first point in your temperature profile is at ', & tempz(1), ' AGL' print *, 'This is a problem since it is assumed to be 0 AGL' print *, 'Either set the first point at 0, or change the ', & 'code to compute the surface density properly in ', & 'this case' stop end if c *** Determine the surface air density if it was not given if( .not. found_surf_rho ) then if( .not. found_runway_elev ) then print *, 'Warning: attempting to compute surf_rho but ', & 'runway_elev is not given.' print *, ' Assuming runway_elev = 0.0' runway_elev = 0.0 end if if( .not. found_msl_pressure ) then print *, 'Warning: attempting to compute surf_rho but ', & 'msl_pressure is not given.' print *, ' Assuming msl_pressure = 29.92 in Hg' msl_pressure = 29.92*inHg2Pa end if if( .not. found_surf_temp ) then print *, 'Warning: attempting to compute surf_rho but ', & 'surf_temp is not given.' print *, ' Assuming surf_temp = 15.0 C' surf_temp = 15.0 + 273.15 found_surf_temp = .true. end if beta = -0.0065 T_standard_0 = 288.15 expon = -grav/(beta*R_gas) arg = 1.0 + beta*runway_elev/T_standard_0 surf_press = msl_pressure*arg**expon surf_rho = surf_press/(R_gas*surf_temp) end if c *** Get profiles for all of the thermodynamic variables rho(1) = surf_rho call compute_state( tempz, rho, press, temp, theta, bvfsq, grav, & R_gas, gamma_gas, flag_pot_temp, i_unstable, n_temp ) c *** Get the air density at the initial aircraft altitude do i=1, n_temp-1 if( z0 .gt. tempz(i) .and. z0 .le. tempz(i+1) ) goto 55 end do 55 rho_z0 = rho(i) + (rho(i+1)-rho(i))/(tempz(i+1)-tempz(i))* & (z0-tempz(i)) c *** Adjust the aircraft speed for the air density ac_speed = ac_speed*sqrt(1.225/rho_z0) c *** Further checks for missing inputs gamma_bool = '000' if( found_gamma0 ) gamma_bool(1:1) = '1' if( found_v0 ) gamma_bool(2:2) = '1' if( found_ac_weight ) gamma_bool(3:3) = '1' select case( gamma_bool ) case( '000' ) print *, 'ERROR: None of ac_weight, gamma0, or v0 are given' stop case( '100' ) ac_weight = gamma0*rho_z0*ac_speed*b0/grav v0 = gamma0/(2.0*pi*b0) case( '010' ) gamma0 = 2.0*pi*b0*v0 ac_weight = gamma0*rho_z0*ac_speed*b0/grav case( '001' ) gamma0 = ac_weight*grav/(rho_z0*ac_speed*b0) v0 = gamma0/(2.0*pi*b0) c print *, 'ac_weight, rho_z0, ac_speed, b0, gamma0 = ' c print *, ac_weight, rho_z0, ac_speed, b0, gamma0 case( '110' ) gamma = 2.0*pi*b0*v0 if( abs(1.0-gamma/gamma0) .gt. 1.0e-3 ) then print *, 'Warning: both gamma0 and v0 are ', & 'given but they appear inconsistent' end if ac_weight = gamma0*rho_z0*ac_speed*b0/grav case( '101' ) gamma = ac_weight*grav/(rho_z0*ac_speed*b0) v0 = gamma0/(2.0*pi*b0) if( abs(1.0-gamma/gamma0) .gt. 1.0e-3 ) then print *, 'Warning: both gamma0 and ac_weight are ', & 'given but they appear inconsistent' end if case( '011' ) gamma0 = 2.0*pi*b0*v0 weight = gamma0*rho_z0*ac_speed*b0/grav if( abs(1.0-weight/ac_weight) .gt. 1.0e-3 ) then print *, 'Warning: both v0 and ac_weight are given ', & 'but they appear inconsistent' end if case( '111' ) gamma1 = 2.0*pi*b0*v0 gamma2 = ac_weight*grav/(rho_z0*ac_speed*b0) if( abs(1.0-gamma1/gamma0)+abs(1.0-gamma2/gamma0) & .gt. 1.0e-3 ) then print *, 'Warning: all three of gamma0, v0 and ac_weight ', & 'are given but they appear inconsistent' end if end select c *** Construct the headwind profile if( n_hw .eq. 0 ) then if( found_surf_wind_spd .and. surf_wind_spd .eq. 0.0 ) then head_wind = 0.0 else if( (.not. found_surf_wind_dir) .or. (.not. found_heading) & .or. (.not. found_surf_wind_spd) ) then print *, 'Warning: attempting to compute the headwind ', & 'but at least one of' print *, 'surf_wind_dir, surf_wind_spd, heading is missing' print *, ' Headwind and crosswind will assumed to be zero' head_wind = 0.0 else angle = ( surf_wind_dir - heading )/180.0*pi head_wind = surf_wind_spd*cos(angle) c print *, c & 'surf_wind_spd, surf_wind_dir, heading, angle, head_wind = ', c & surf_wind_spd, surf_wind_dir, heading, angle, head_wind end if n_hw = 2 hwz(1) = 0.0 hwz(2) = z_top hw(1) = head_wind hw(2) = head_wind end if if( n_hw .eq. 2 .and. abs(hw(2)-hw(1)) .lt. 1.e-6 ) then no_prof_hw = .true. end if c *** Construct the crosswind profile if( n_cw .le. 1 ) then if( n_cw .eq. 0 ) then if( (.not. found_surf_wind_dir) .or. (.not. found_heading) & .or. (.not. found_surf_wind_spd) ) then cross_wind = 0.0 else angle = ( surf_wind_dir - heading )/180.0*pi cross_wind = surf_wind_spd*sin(angle) end if else cross_wind = cw(1) end if n_cw = 2 cwz(1) = 0.0 cwz(2) = z_top cw(1) = cross_wind cw(2) = cross_wind end if if( n_cw .eq. 2 .and. abs(cw(2)-cw(1)) .lt. 1.e-6 ) then no_prof_cw = .true. end if c *** Delete a crosswind value at the surface if it is simply a repeat c *** of the first measurement station if( n_cw .ge. 3 .and. cw(1) .eq. cw(2) ) then do i=2, n_cw+1 cw( i-1) = cw( i) cwz(i-1) = cwz(i) end do n_cw = n_cw - 1 end if c *** construct the eddy dissipation rate profile if( n_q .eq. 0 ) then if( .not. found_turbulence ) then print *, 'Warning: attempting to compute turbulence ', & 'profile but turbulence is missing' print *, ' Assuming turbulence = 0.0' turbulence = 0.0 end if n_q = 2 qz(1) = 0.0 qz(2) = z_top q(1) = turbulence q(2) = turbulence end if if( n_q .eq. 2 .and. abs(q(2)-q(1)) .lt. 1.e-6 ) then no_prof_q = .true. end if c print *, 'rho_z0, ac_speed, b0, ac_weight, grav, v0, gamma0 = ' c print *, rho_z0, ac_speed, b0, ac_weight, grav, v0, gamma0 c *** Correct the edr profile for Frankfurt OGE cases if( index(fname,'Fra04_OGE') .ne. 0 ) then print *, 'Fixing up edr for Fra04_OGE' do i=1, n_q edr_star = ((q(i)*b0)**c13)/v0 edr = (max((edr_star-0.10),1.0e-3)*v0)**3/b0 q(i) = edr end do end if c *** extrapolate surface wind from 10m height for Frankfurt IGE cases if( index(fname,'Fra04_IGE') .ne. 0 ) then if ( n_cw .ge. 4 ) then gg2 = (cwz(2)-cwz(1))/(cwz(3)-cwz(1))-(cwz(4)-cwz(1)) & *(cwz(4)-cwz(1))/(cwz(2)-cwz(1))/(cwz(3)-cwz(1)) gg3 = (cwz(2)-cwz(1))/(cwz(4)-cwz(1))-(cwz(3)-cwz(1)) & *(cwz(3)-cwz(1))/(cwz(2)-cwz(1))/(cwz(4)-cwz(1)) ff2 = gg2*(cwz(2)-cwz(3)) ff3 = gg3*(cwz(2)-cwz(4)) cw(1) = ( (cw(2)*(cwz(4)-cwz(1))-cw(4)*(cwz(2)-cwz(1)))*gg3 & - (cw(2)*(cwz(3)-cwz(1))-cw(3)*(cwz(2)-cwz(1)))*gg2 & )/(ff2-ff3) else if ( n_cw .eq. 3 ) then cw(1) = cw(2) - (cw(3)-cw(2))/(cwz(3)-cwz(2))*cwz(2) endif endif ! if ( n_cw .ge. 3 ) then ! cw(1) = cw(2) - (cw(3)-cw(2))/(cwz(3)-cwz(2))*cwz(2) ! endif if (cw(1)*cw(2) .lt. 0.) then cw(1) = 0.0 endif ! calibrate the crosswind with an initial proxy yp0 = y_meas_p(1) + (y_meas_p(2)-y_meas_p(1))* & (0.-time_meas_p(1))/(time_meas_p(2)-time_meas_p(1)) ys0 = y_meas_s(1) + (y_meas_s(2)-y_meas_s(1))* & (0.-time_meas_s(1))/(time_meas_s(2)-time_meas_s(1)) zp0 = z_meas_p(1) + (z_meas_p(2)-z_meas_p(1))* & (0.-time_meas_p(1))/(time_meas_p(2)-time_meas_p(1)) zs0 = z_meas_s(1) + (z_meas_s(2)-z_meas_s(1))* & (0.-time_meas_s(1))/(time_meas_s(2)-time_meas_s(1)) if (nt_meas_p .ge. 3) then vyp=(y_meas_p(3)-y_meas_p(1))/(time_meas_p(3)-time_meas_p(1)) vzp=(z_meas_p(3)-z_meas_p(1))/(time_meas_p(3)-time_meas_p(1)) else vyp=(y_meas_p(2)-y_meas_p(1))/(time_meas_p(2)-time_meas_p(1)) vzp=(z_meas_p(2)-z_meas_p(1))/(time_meas_p(2)-time_meas_p(1)) endif if (nt_meas_s .ge. 3) then vys=(y_meas_s(3)-y_meas_s(1))/(time_meas_s(3)-time_meas_s(1)) vzs=(z_meas_s(3)-z_meas_s(1))/(time_meas_s(3)-time_meas_s(1)) else vys=(y_meas_s(2)-y_meas_s(1))/(time_meas_s(2)-time_meas_s(1)) vzs=(z_meas_s(2)-z_meas_s(1))/(time_meas_s(2)-time_meas_s(1)) endif do i=1,n_cw-1 if ( cwz(i).lt.zp0 .and. cwz(i+1).gt.zp0 ) then cwp0 = cw(i) + (cw(i+1)-cw(i)) & *(zp0-cwz(i))/(cwz(i+1)-cwz(i)) endif if ( cwz(i).lt.zs0 .and. cwz(i+1).gt.zs0 ) then cws0 = cw(i) + (cw(i+1)-cw(i)) & *(zs0-cwz(i))/(cwz(i+1)-cwz(i)) endif enddo i_base_beg = scan( fname, '/\', .true. ) + 1 i_base_end = scan( fname, '.', .true. ) - 1 if( i_base_end .eq. -1 ) i_base_end = index(fname,' ') - 1 base_name = fname(i_base_beg:i_base_end) open(unit=25,file='Output/'//trim(base_name)//'.meas_init') write(25,126) write(25,57) 'y0p = ',yp0 write(25,57) 'y0s = ',ys0 write(25,57) 'z0p = ',zp0 write(25,57) 'z0s = ',zs0 write(25,57) 'w0p = ',vzp write(25,57) 'w0s = ',vzs write(25,57) 'v0p = ',vyp write(25,57) 'v0s = ',vys write(25,57) 'CW0p = ',cwp0 write(25,57) 'CW0s = ',cws0 write(25,58) 'wind ratio p = ',vyp/cwp0 write(25,58) 'wind ratio s = ',vys/cws0 57 format(a7,f8.3) 58 format(a15,f8.3) 126 format('# measured data initial conditions') close(25) vel_factor = (vyp/cwp0+vys/cws0)*0.5 vel_factor = (vyp+vys)/(cwp0+cws0) if (abs(cwp0+cws0) < 0.10) then vel_factor = 1.0 endif c write(6,*)"vel_factor=",vel_factor do i=1,n_cw cw(i) = cw(i)*vel_factor enddo do i=1,n_hw ! hw(i) = hw(i)*sign(1.,vel_factor) hw(i) = hw(i)*(vel_factor) enddo end if c *** Calculate the average headwind for use in normalizing the data z_min = 0.0 + float(nt_meas_p+nt_meas_s)*1.0e+06 do i=1, nt_meas_p z_min = min( z_min, z_meas_p(i) ) end do do i=1, nt_meas_s z_min = min( z_min, z_meas_s(i) ) end do sum2 = 0.0 isamp = 0 do i=1, n_hw if( hwz(i) .ge. z_min .and. hwz(i) .le. z0 ) then sum2 = sum2 + hw(i) isamp = isamp + 1 end if end do if( isamp .eq. 0 ) then hw_avg = hw(1) else hw_avg = sum2/float(isamp) end if c *** The initial vortex core radii are not currently in the ascii data c *** file. These will be added once I figure out how to estimate their c *** values for cases where they were not measured. For now we set them c *** to a crude guess. the current estimate is from memphis field data c *** for 727's and from NWRA lab data emailed by don 4 april 2013. c r0p = 0.0625*ac_span ! in meters c r0s = 0.0625*ac_span ! in meters c *** The vortex core radius is estimated here. We use 3.0% of the wingspan. R0 = 0.03*ac_span R0p = R0 R0s = R0 return end