subroutine Solve_for_Aurora use ModQset use ModSize use ModAMIE use ModData use ModConductance use ModIndices use ModConstants use ModFiles, only : NCond_Files, lun_dataoutput, cond_file implicit none integer :: i, ith, ix, j, m, mm, n, iFile, iTempCounter real :: x, xm1, x0, xp1, xp2, sp, cp, sza real :: ped, sped, tped, hal, shal, thal, avkev, eflx real :: pedguess, halguess real :: xnorm_avee, xnorm_eflx real, dimension(MAX_DATA) :: obsp, obsh, obsa, obse, & avee_obs, eflx_obs, obs_dummy, xwt_dummy, diff_dummy, & diff_avee, diff_eflx, diff_hal, diff_ped, & dcd_avee, dcd_eflx, dcd_ped real, dimension(MAX_DATA) :: avee_wt, eflx_wt, er_dummy, & cond_wt real, dimension(0:ithmx) :: eflx1d, ekev1d, ped1d, hal1d real, dimension(MAX_DATA, mcoef) :: datfma real, dimension(mcoef) :: solution_avee, solution_eflx, & solution_hal, solution_ped, solution_dummy, & asum_dummy, ao_dummy, ps_dummy, v_dummy real, dimension(0:lonmx,0:ithmx) :: & dummy_err_east, dummy_err_north real :: MagLat, MagLT, DeltaT, TimeErrorAveE, TimeErrorEFlux, & TimeErrorPed, TimeErrorHall n = iteration_number NAuroralObs = 1 do j = 1, NAhn call get_auroral_conductance(AhnLat(j), AhnMLT(j), Indices(n,hpi_norm_), & ped, hal, avkev, eflx) if (solar_model == 'old') then call get_old_solar_cond(AhnLat(j)*ihsoln, AhnMLT(j), sped, shal, sza) else call get_new_solar_cond( Indices(n,f107_), AhnLat(j)*ihsoln, & AhnMLT(j), sped, shal, sza) endif tped = sqrt(ped*ped + sped*sped) thal = sqrt(hal*hal + shal*shal) if (AhnPed(j) > sped .and. AhnHal(j) > shal) then pedguess = sqrt(AhnPed(j)*AhnPed(j) - sped*sped) halguess = sqrt(AhnHal(j)*AhnHal(j) - shal*shal) ! if (ped/tped > 0.25 .and. hal/thal > 0.25) then ! if (/tped > 0.25 .and. hal/thal > 0.25) then ! pedguess = ped/tped * AhnPed(j) ! halguess = hal/thal * AhnHal(j) x = halguess / pedguess avee_obs(NAuroralObs) = exp(alog(x/0.45)/0.85) eflx_obs(NAuroralObs) = & (pedguess*(16.0+avee_obs(NAuroralObs)**2)/(40.0*avee_obs(NAuroralObs)))**2 avee_wt(NAuroralObs) = AhnWeight_AveE(j) eflx_wt(NAuroralObs) = AhnWeight_EFlux(j) cond_wt(NAuroralObs) = AhnWeight(j) if (avee_obs(NAuroralObs) < ConductanceBackground(avee_)) & avee_obs(NAuroralObs) = ConductanceBackground(avee_) if (eflx_obs(NAuroralObs) < ConductanceBackground(eflux_)) & eflx_obs(NAuroralObs) = ConductanceBackground(eflux_) obsp(NAuroralObs) = alog(AhnPed(j)/tped) obsh(NAuroralObs) = alog(AhnHal(j)/thal) obsa(NAuroralObs) = alog(avee_obs(NAuroralObs)/avkev) obse(NAuroralObs) = alog(eflx_obs(NAuroralObs)/eflx) sp = sin(AhnMLT(j)*htor) cp = cos(AhnMLT(j)*htor) MagLat = AhnLat(j) ! load the matrix FF: call fcmp(mmx,1,cp,sp,ff) call Fill_datfma NAuroralObs = NAuroralObs + 1 endif enddo NAuroralObs = NAuroralObs - 1 if (OutputData .and. NAuroralObs > 0) then write(lun_dataoutput, *) "" write(lun_dataoutput, "(a4)") "#AHN" write(lun_dataoutput, *) NAuroralObs do j=1,NAuroralObs avkev = avee_obs(j) / (10**(obsa(j))) eflx = eflx_obs(j) / (10**(obse(j))) write(lun_dataoutput, "(2f6.2,1p,6e9.2)") AhnLat(j), AhnMLT(j), & avee_obs(j), avee_wt(j), avkev, eflx_obs(j), eflx_wt(j), eflx enddo endif ! NAuroralObs = NAhn if (OutputData .and. ncond_files > 0) then write(lun_dataoutput, *) "" write(lun_dataoutput, "(a13)") "#CONDUCTANCES" write(lun_dataoutput, *) ncond_files endif do ifile = 1, ncond_files if (OutputData) then iTempCounter = 0 do j = 1, npts_cond(ifile) DeltaT = abs(ConductanceTimes(ifile, j) - currenttime) if (DeltaT <= window_dt/2.0) iTempCounter = iTempCounter + 1 enddo write(lun_dataoutput, "(a)") cond_file(ifile) write(lun_dataoutput, *) iTempCounter endif do j = 1, npts_cond(ifile) DeltaT = abs(ConductanceTimes(ifile, j) - currenttime) if (DeltaT <= window_dt/2.0) then NAuroralObs = NAuroralObs + 1 MagLat = ConductanceLocations(ifile, j, mag_lat_) MagLT = ConductanceLocations(ifile, j, mlt_) call get_auroral_conductance(MagLat, MagLT, & Indices(n,hpi_norm_), ped, hal, avkev, eflx) sped = 0.0 shal = 0.0 ! if (solar_model == 'old') then ! call get_old_solar_cond(AhnLat(j)*ihsoln, AhnMLT(j), sped, shal, sza) ! else ! call get_new_solar_cond( Indices(n,f107_), AhnLat(j)*ihsoln, & ! AhnMLT(j), sped, shal, sza) ! endif tped = sqrt(ped*ped + sped*sped) thal = sqrt(hal*hal + shal*shal) obsp(NAuroralObs) = alog(ConductanceData(ifile, j, pedersen_)/tped) obsh(NAuroralObs) = alog(ConductanceData(ifile, j, hall_)/thal) obsa(NAuroralObs) = alog(ConductanceData(ifile, j, avee_)/avkev) obse(NAuroralObs) = alog(ConductanceData(ifile, j, eflux_)/eflx) if (LinearDropOffConductance) then TimeErrorAveE = DeltaT / LinearDropOffTimeConductance * LinearDropOffAveE TimeErrorEFlux = DeltaT / LinearDropOffTimeConductance * LinearDropOffEFlux TimeErrorPED = DeltaT / LinearDropOffTimeConductance TimeErrorHALL = DeltaT / LinearDropOffTimeConductance else if (DeltaT > MaxeFoldingDeltaT) DeltaT = MaxeFoldingDeltaT TimeErrorAveE = ConductanceDataError(ifile, j, avee_) * & exp(DeltaT / ExponentialDropOffeFolding) TimeErrorEFlux = ConductanceDataError(ifile, j, eflux_) * & exp(DeltaT / ExponentialDropOffeFolding) endif ! avee_wt(NAuroralObs) = & avee_wt(NAuroralObs) = ConductanceData(ifile, j, avee_)**2* & 1.0/(ConductanceDataError(ifile, j, avee_)**2 + TimeErrorAveE**2) ! eflx_wt(NAuroralObs) = & eflx_wt(NAuroralObs) = ConductanceData(ifile, j, eflux_)**2* & 1.0/(ConductanceDataError(ifile, j, eflux_)**2+ TimeErrorEFlux**2) cond_wt(NAuroralObs) = & 0.5*(ConductanceData(ifile, j, pedersen_)**2/ & (ConductanceDataError(ifile, j, pedersen_)**2 + & TimeErrorPED**2) & + ConductanceData(ifile, j, hall_)**2/ & (ConductanceDataError(ifile, j, hall_)**2 + TimeErrorHALL**2)) if (OutputData) & write(lun_dataoutput, "(2f6.2,1p,7e9.2)") MagLat, MagLT, & ConductanceData(ifile, j, avee_), avee_wt(NAuroralObs),avkev, & ConductanceData(ifile, j, eflux_), eflx_wt(NAuroralObs),eflx, & DeltaT sp = sin(MagLT*htor) cp = cos(MagLT*htor) ! load the matrix FF: call fcmp(mmx,1,cp,sp,ff) call Fill_datfma endif enddo enddo xwt_dummy = 1.0 ao_dummy = 1.0 obs_dummy = 1.0 if (NAuroralObs < 2) NAuroralObs = 0 call matsolv(3, NAuroralObs, datfma, avee_wt, xwt_dummy, obsa, & obs_dummy, er_dummy, asum_dummy, ao_dummy, solution_avee, & solution_dummy, diff_avee, diff_dummy, & dcd_avee, dummy_err_north, dummy_err_east, & rms_avee, err_avee) call matsolv(4, NAuroralObs, datfma, eflx_wt, xwt_dummy, obse, & obs_dummy, er_dummy,asum_dummy,ao_dummy,solution_eflx, & solution_dummy, diff_eflx, diff_dummy, & dcd_eflx, dummy_err_north, dummy_err_east, & rms_eflx, err_eflx) ! Need to figure out rmsaur... Calculated_Average_Energy = 0. Calculated_Energy_Flux = 0. do m=-mmx,mmx ekev1d = 0. eflx1d = 0. mm = iabs(m) do i=ibm(m),iem(m) ix = i - ibm(m) + (ns(mm) - 1)/2 + 1 ! qsclr is only defined to ITHPLT, not ITHMX ! qsclr is defined to ITHMX, but is reasonable to ITHTRNS+~2 (4/94) do ith=0,ithtrns ekev1d(ith) = ekev1d(ith) + solution_avee(i)*qsclr(ix,ith) eflx1d(ith) = eflx1d(ith) + solution_eflx(i)*qsclr(ix,ith) enddo enddo do j=1,lonmx do ith=0,ithtrns Calculated_Average_Energy(j,ith) = & Calculated_Average_Energy(j,ith)+ekev1d(ith)*amie_f(m,j) Calculated_Energy_Flux(j,ith) = & Calculated_Energy_Flux(j,ith)+eflx1d(ith)*amie_f(m,j) enddo enddo enddo ! Convert fit from log form do ith=0,ithmx do j=1,lonmx Calculated_Average_Energy(j,ith) = & Model_Average_Energy(j,ith)*exp(Calculated_Average_Energy(j,ith)) Calculated_Energy_Flux(j,ith) = & Model_Energy_Flux(j,ith)*exp(Calculated_Energy_Flux(j,ith)) enddo Calculated_Average_Energy(0,ith) = Calculated_Average_Energy(lonmx,ith) Calculated_Energy_Flux(0,ith) = Calculated_Energy_Flux(lonmx,ith) enddo HPI_Calculated(n) = 0.0 do i = 0, ithmx do j = 1, lonmx HPI_Calculated(n) = HPI_Calculated(n) + & Calculated_Energy_Flux(j,i) * 1.0e-3 * & st(i) * Ri * Ri * dlon * dth * 1.0e-9 enddo enddo Indices(n, hpi_calc_) = HPI_Calculated(n) Calculated_Auroral_Pedersen = SQRT(Calculated_Energy_Flux) * 40. * & Calculated_Average_Energy / (16. + Calculated_Average_Energy**2) Calculated_Auroral_Hall = & 0.45*(Calculated_Average_Energy**0.85)*Calculated_Auroral_Pedersen where (Calculated_Auroral_Pedersen > 100.0) Calculated_Auroral_Pedersen = 100.0 where (Calculated_Auroral_Hall > 100.0) Calculated_Auroral_Hall = 100.0 !!!! fitting to the conductance data !!! call matsolv(2, NAuroralObs, datfma, cond_wt, xwt_dummy, obsp, & !!! obsh, er_dummy, asum_dummy, ao_dummy, solution_ped, & !!! solution_hal, diff_ped, diff_hal, & !!! dcd_ped, dummy_err_north, dummy_err_east, & !!! rms_avee, err_avee) !!! !!! Calculated_Auroral_Pedersen = 0. !!! Calculated_Auroral_Hall = 0. !!! !!! do m=-mmx,mmx !!! !!! ped1d = 0. !!! hal1d = 0. !!! !!! mm = iabs(m) !!! !!! do i=ibm(m),iem(m) !!! !!! ix = i - ibm(m) + (ns(mm) - 1)/2 + 1 !!! !!! do ith=0,ithtrns !!! ped1d(ith) = ped1d(ith) + solution_ped(i)*qsclr(ix,ith) !!! hal1d(ith) = hal1d(ith) + solution_hal(i)*qsclr(ix,ith) !!! enddo !!! !!! enddo !!! !!! do j=1,lonmx !!! do ith=0,ithtrns !!! Calculated_Auroral_Pedersen(j,ith) = & !!! Calculated_Auroral_Pedersen(j,ith)+ped1d(ith)*amie_f(m,j) !!! Calculated_Auroral_Hall(j,ith) = & !!! Calculated_Auroral_Hall(j,ith)+hal1d(ith)*amie_f(m,j) !!! enddo !!! enddo !!! !!! enddo !!! !!!! Convert fit from log form !!! !!! do ith=0,ithmx !!! !!! do j=1,lonmx !!! Calculated_Auroral_Pedersen(j,ith) = & !!! Model_Auroral_Pedersen(j,ith)*exp(Calculated_Auroral_Pedersen(j,ith)) !!! Calculated_Auroral_Hall(j,ith) = & !!! Model_Auroral_Hall(j,ith)*exp(Calculated_Auroral_Hall(j,ith)) !!! enddo !!! !!! Calculated_Auroral_Pedersen(0,ith) = Calculated_Auroral_Pedersen(lonmx,ith) !!! Calculated_Auroral_Hall(0,ith) = Calculated_Auroral_Hall(lonmx,ith) !!! !!! enddo Calculated_Total_Pedersen = & sqrt(Calculated_Auroral_Pedersen**2 + Model_Solar_Pedersen**2) Calculated_Total_Hall = & sqrt(Calculated_Auroral_Hall**2 + Model_Solar_Hall**2) contains subroutine Fill_datfma x = (90.-abs(MagLat))/(rad*dth) ith = x ith = max0(1,min0(ithtrns-2,ith)) x = x - ith xm1 = x*(-2. + x*(3. - x))/6. x0 = 1. + x*(-.5 + x*(-1. + .5*x)) xp1 = x*(1. + x*(.5 - .5*x)) xp2 = x*(-1. + x*x)/6. do m=-mmx,mmx mm = iabs(m) ! assume that there are equal numbers of symmetric and antisymmetric ! q's for a given value of m. ns(mm) gives the first index value ! for the first subscript of q for wavenumber mm. cubic interpolation. do i=ibm(m),iem(m) ix = i - ibm(m) + nss(mm) datfma(NAuroralObs,i) = (xm1*qsclr(ix,ith-1)+ & x0*qsclr(ix,ith)+ & xp1*qsclr(ix,ith+1)+ & xp2*qsclr(ix,ith+2))*ff(m) enddo enddo end subroutine Fill_datfma end subroutine Solve_for_Aurora subroutine Conductance_Cu_set use ModQset use ModAMIE use ModConductance use ModConstants implicit none real :: xko2, gam2, pcwt, area, x, xn integer :: ithpcb, ith, ix, i, j, k, m, mm real, dimension(ithmx) :: wt real, dimension(mcoef) :: z integer, dimension(mcoef) :: ipvt xko2 = (ri/3.e5)**2 gam2 = 4.0 pcwt = 0.1 ithpcb = ifix(float(ithmx) * 20. / 90. ) area = 0.0 wt = 0.0 ax = 0.0 ccu = 0.0 do ith=1,ithpcb-1 x = ct(ith)-ct(ithpcb) wt(ith) = (st(ith) + .2)*x*x area = area + wt(ith) enddo do ith=1,ithpcb-1 wt(ith) = wt(ith)*pcwt/area enddo area = 0. do m=-mmx,mmx mm = iabs(m) do ith=1,ithpcb-1 do i=ibm(m),iem(m) dqs(i) =qsclr(i-ibm(m)+nss(mm),ith) enddo do i=ibm(m),iem(m) do j=i,iem(m) ccu(i,j) = ccu(i,j) + wt(ith)*(dqs(i)*dqs(j)) enddo enddo enddo do k=ibm(m),iem(m) ix = (k-ibm(m)) +(ns(mm)-1)/2+1 xn = wnt2s(ix) + gam2*wnp2s(ix) ! previous to 9/91, amie assumed a power spectrum of the form 1/k**4. ! x = xko2*xko2 ! ccu(k,k) = (xn*xn + x)/x + ccu(k,k) ! analysis of idm oct 1981 data showed a power spectrum of 1/k**2 x = xko2 ccu(k,k) = (xn + x)/x + ccu(k,k) ccuid(k) = ccu(k,k) do i=1,k ax(i,k) = ccu(i,k) enddo do i=1,k ax(k,i) = ax(i,k) enddo ax(k,k+imx) = 1. enddo enddo call solver(ax,imx,2*imx,mcoef,ipvt,z) do i=1,imx do j=1,i ccu(i,j) = .5*( ax(i,j+imx) + ax(j,i+imx)) enddo enddo ! lower-left part of cu (including diagonal) contains cu. ! test symmetry of cu ! do i=10,imx,10 ! do j=10,i,10 ! write(6,'(1x,2i4,2e10.3)') i,j,ccu(i,j),ax(j,i+imx) ! enddo ! enddo ! test whether ccu is physical do i=1,imx do j=1,i-1 if (ccu(i,j)*ccu(i,j).gt.ccu(i,i)*ccu(j,j)) then write(6,3025) i,j,ccu(i,j),ccu(i,i),ccu(j,j) 3025 format(1x,'i=',i4,'j=',i4,'ccu(i,j)=',e10.3, & 'ccu(i,i)=',e10.3,'ccu(j,j)=',e10.3) endif enddo enddo end subroutine Conductance_Cu_set