subroutine Solve_for_EField use ModQset use ModAMIE use ModData use ModConductance, only : Calculated_Total_Pedersen, Calculated_Total_Hall use ModEField use ModIndices use ModConstants use ModFiles implicit none real, dimension(mcoef) :: ao, a, asum, ah, qa, dqa, hmg, dmg, zmg, ymg real, dimension(mcoef,mcoef) :: tfma, psfma, zfma, vfma real, dimension(-mmx:mmx) :: curtm, curpm real, dimension(lonmx) :: curp, curt real, dimension(MAX_DATA, mcoef) :: datfma, efma integer :: ith, m, mm, mii, k, ix, lon, i, j, ibmm1, ip, n, ifile integer :: NEfield, ixx, iTempCounter, idel, idel2, nE, nBad, nGood, iDir integer :: nInversions real :: rst, rdt, rdtst, curtx, curpx, et, ep, sth real :: x, xm1, x0, xp1, xp2, sp, cp, sza, dist real :: fint, dfint, dpotdp, dpotdt real :: cos_dec, sin_dec, colat, sin_se, cos_se real :: eflx, avkev, thal_sqrd, hal, shal, ped, sped real :: ztemp, ztempwt real :: MagError(3) real, dimension(MAX_DATA) :: xwt_dummy, efield_wt, efield_obs, obsh real, dimension(0:ithmx,NEfield_Solutions) :: EField_Coef real, dimension(mcoef,NEfield_Solutions) :: Magnetic_Coef real, dimension(MAX_DATA) :: diff_dummy, diff_pot, dcd_pot, er real, dimension(0:lonmx,0:ithmx) :: err_east, err_north real :: MagLat, MagLT, DeltaT, TimeErrorEfield, TimeErrorCurrent, TotalError real :: clatsap, dclat, dxlat, dxlon, dxmlt, signyz, dalat, dalon logical :: IsFound, IsWorse, IsDone real :: OldError, NewError integer :: iMaxErrorFlag iMaxErrorFlag = 1 n = iteration_number nBad = 100 nInversions = 0 IsDone = .false. IsWorse = .false. OldError = 1.0e32 do while (.not.IsDone) if (DebugLevel > 1) & write(*,*) "==> Starting an Electric Field Inversion!" ! This condition is only true if we have completed 2 inversions already: ! Basically, we want to go back to the last inversion if the current ! inversion is worse. This means that we have tightened up the ! magnetometer error flag checking and we are going to just call it ! done. IsDone = .true. if (IsWorse) IsDone = .true. xwt_dummy = 1.0 ! COMPUTE COEFFICIENTS FOR MODEL ELECTRIC FIELD, AO ao(1:imx) = 0.0 ! ENSMDL is the negative theta component, and EESMDL is the phi component ! of the model E field. ! LOOP IN COLATITUDE NEField = 0 do ith=1,ithmx rst = ri*st(ith) rdt = ri*dth if (ith == ithmx) rdt = .5*rdt rdtst = rdt*st(ith) do m=-mmx,mmx mm = iabs(m) ! ASSUME THAT THERE ARE EQUAL NUMBERS OF SYMMETRIC AND ANTISYMMETRIC ! QS FOR A GIVEN VALUE OF M. NS(MM) GIVES THE FIRST INDEX VALUE ! FOR THE FIRST SUBSCRIPT OF Q FOR WAVENUMBER MM. do k=ibm(m),iem(m) ix = k - ibm(m) dqs(k) = dq(2*ix+ns(mm),ith) qs(k) = q(2*ix+ns(mm),ith) enddo enddo ! CALCULATE FOURIER COEFFICIENTS OF E FIELD WITH RESPECT TO LONGITUDE. curtm(-mmx:mmx) = 0.0 curpm(-mmx:mmx) = 0.0 do lon=1,lonmx curtx = -EField_BGModel(LON,ITH, efield_north_)/lonmx curpx = EField_BGModel(LON,ITH, efield_east_)/lonmx curtm(-mmx:mmx) = curtm(-mmx:mmx) + curtx*amie_f(-mmx:mmx,lon) curpm(-mmx:mmx) = curpm(-mmx:mmx) + curpx*amie_f(-mmx:mmx,lon) enddo do m=-mmx,mmx ao(ibm(m):iem(m)) = ao(ibm(m):iem(m)) - & dqs(ibm(m):iem(m))*curtm(m)*rdtst - & curpm(-m)*qs(ibm(m):iem(m))*m*rdt enddo enddo ! Percent_Error = 100.0 a(1:imx) = ao(1:imx) ! ! ********************* Solve for electric field ******************* ! WRITE(6,"(1X,'CALCULATING PSFMA',)") ! CALCULATE MATRIX PSFMA, WHICH CONVERTS ELECTRIC POTENTIAL ! COEFFICIENTS (A) TO TOROIDAL IONOSPHERIC CURRENT FUNCTION ! COEFFICIENTS (PS), AND MATRIX TFMA, WHICH CONVERTS A TO POLOIDAL ! IONOSPHERIC CURRENT FUNCTION COEFFICIENTS (T). tfma(1:imx,1:imx) = 0.0 psfma(1:imx,1:imx) = 0.0 ! LOOP IN COLATITUDE DO ITH=1,ITHMX rst = ri*st(ith) rdt = ri*dth if (ith == ithmx) rdt = .5*rdt rdtst = rdt*st(ith) 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. do k=ibm(m),iem(m) ix = k - ibm(m) dqs(k) = dq(2*ix+ns(mm),ith) qs(k) = q(2*ix+ns(mm),ith) dqa(k) = dq(2*ix+ns(mm)+1,ith) qa(k) = q(2*ix+ns(mm)+1,ith) enddo enddo DO I=1,IMX MII = MI(I) ! USE SIN(INCLINATION)=1 HERE. curt(1:lonmx) = -Calculated_Total_Pedersen(1:lonmx,ITH)* & DQS(I)*amie_F(MII,1:lonmx)/RI - & Calculated_Total_Hall(1:lonmx,ITH)* & MII*QS(I)*amie_F(-MII,1:lonmx)/RST curp(1:lonmx) = -Calculated_Total_Pedersen(1:lonmx,ITH)* & MII*QS(I)*amie_F(-MII,1:lonmx)/RST + & Calculated_Total_Hall(1:lonmx,ITH)* & DQS(I)*amie_F(MII,1:lonmx)/RI ! CALCULATE FOURIER COEFFICIENTS OF CURRENT WITH RESPECT TO LONGITUDE. curtm(-mmx:mmx) = 0.0 curpm(-mmx:mmx) = 0.0 do lon=1,lonmx curtx = curt(lon)/lonmx curpx = curp(lon)/lonmx curtm(-mmx:mmx) = curtm(-mmx:mmx) + curtx*amie_f(-mmx:mmx,lon) curpm(-mmx:mmx) = curpm(-mmx:mmx) + curpx*amie_f(-mmx:mmx,lon) enddo do m=-mmx,mmx tfma(ibm(m):iem(m),i) = tfma(ibm(m):iem(m),i) - & dqs(ibm(m):iem(m))*curtm(m)*rdtst - & curpm(-m)*qs(ibm(m):iem(m))*m*rdt psfma(ibm(m):iem(m),i) = psfma(ibm(m):iem(m),i) + & dqa(ibm(m):iem(m))*curpm(m)*rdtst - & curtm(-m)*qa(ibm(m):iem(m))*m*rdt enddo enddo enddo ! CALCULATE MATRIX VFMA, WHICH CONVERTS ELECTRIC POTENTIAL ! COEFFICIENTS (A) TO GROUND MAGNETIC POTENTIAL COEFFICIENTS (V). do m=-mmx,mmx ibmm1 = ibm(m) - 1 mm = iabs(m) do i=ibm(m),iem(m) ip = i - ibm(m) + nss(mm) ! nss(mm) is starting row in vfmps. zfma(i,1:imx) = 0. vfma(i,1:imx) = 0. do k=ibm(m),iem(m) zfma(i,1:imx) = zfma(i,1:imx) + zfmps(ip,k-ibmm1)*psfma(k,1:imx) vfma(i,1:imx) = vfma(i,1:imx) + vfmps(ip,k-ibmm1)*psfma(k,1:imx) enddo enddo enddo if (Use_Magnetometers) then NEField = 1 do j=1,nmagnetometers ! IsFound = .false. ! if (mglat(j)*ihsoln < 0 .or. hgm(n,j) > 9998.0) IsFound = .true. ! do iDir = 1, 3 ! if (iMagErrorFlag(n,j,iDir) > iMaxErrorFlag) IsFound = .true. ! enddo ! if (.not. IsFound) then if (mglat(j)*ihsoln > 0) then Cos_Dec = cos(decmg(j)*dtor) Sin_Dec = sin(decmg(j)*dtor) Colat = (90.0-abs(mglat(j)))*dtor call get_auroral_conductance(mglat(j), amltmg(j), & Indices(n,hpi_norm_), ped, hal, avkev, eflx) if (solar_model == 'old') then call get_old_solar_cond(mglat(j)*ihsoln, amltmg(j), & sped, shal, sza) else call get_new_solar_cond(Indices(n,f107_), & mglat(j)*ihsoln, amltmg(j), sped, shal, sza) endif thal_sqrd = hal*hal + shal*shal if ( iMagErrorFlag(n,j,1) <= iMaxErrorFlag .and. & iMagErrorFlag(n,j,2) <= iMaxErrorFlag) then efield_wt(NEfield) = & wth(n,j) * xwt_dummy(j)/(1.0 + wth(n,j)/thal_sqrd) efield_wt(NEfield+1) = & wtd(n,j) * xwt_dummy(j)/(1.0 + wtd(n,j)/thal_sqrd) else efield_wt(NEfield) = 0.0 efield_wt(NEfield+1) = 0.0 endif efield_obs(NEfield) = hgm(n,j) efield_obs(NEfield+1) = dgm(n,j) datfma(NEfield,1:imx) = 0.0 datfma(NEfield+1,1:imx) = 0.0 sp = sin(amltmg(j)*htor) cp = cos(amltmg(j)*htor) call fcmp(mmx,1,cp,sp,ff) x = colat/dth ith = x ith = max(1,min(ithmx-2,ith)) x = x - float(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. sth = sin(colat) do m=-mmx,mmx fint = ff(m)/ri dfint = m*ff(-m)/sth/ri ixx = ns(iabs(m)) + 1 do i=ibm(m),iem(m) ix = i - ibm(m) x = xm1*q(2*ix+ixx,ith-1)+ & x0*q(2*ix+ixx,ith)+ & xp1*q(2*ix+ixx,ith+1)+ & xp2*q(2*ix+ixx,ith+2) dpotdp = x*dfint zmg(i) = x*ff(m) dpotdt = (xm1*dq(2*ix+ixx,ith-1)+ & x0*dq(2*ix+ixx,ith)+ & xp1*dq(2*ix+ixx,ith+1)+ & xp2*dq(2*ix+ixx,ith+2))*fint ! cd = COS(ANGLE BETWEEN NORTH AND OBSERVATION DIRECTION, ! POS. TOWARDS EAST). SD = SIN OF THIS ANGLE. ! FACTOR 1.E9 IS TO CONVERT TESLA TO NANOTESLA. hmg(i) = (Cos_Dec*dpotdt - Sin_Dec*dpotdp)*1.E9 dmg(i) = (-Cos_Dec*dpotdp - Sin_Dec*dpotdt)*1.E9 enddo enddo DO I=1,IMX datfma(NEfield,1:imx) = datfma(NEfield,1:imx) + & hmg(i)*vfma(i,1:imx) datfma(NEField+1,1:imx) = datfma(NEField+1,1:imx) + & dmg(i)*vfma(i,1:imx) enddo NEField = NEField + 2 if (use_zed) then if (iMagErrorFlag(n,j,3) <= iMaxErrorFlag) then Efield_WT(NEField) = wtzed(n,j)*xwt_dummy(J) else Efield_WT(NEField) = 0.0 endif Efield_obs(NEField) = zgm(n,J) datfma(NEField,1:IMX) = 0.0 do i=1,imx datfma(NEField,1:imx) = datfma(NEField,1:imx) + & zmg(i)*zfma(i,1:imx) enddo NEField = NEField + 1 endif endif enddo NEField = NEField - 1 if (NEField > Max_Data) then write(*,*) "Too many Data points in efield_solver!" call stop_mpi("Yikes!") endif endif if (OutputData .and. nefield_files > 0) then write(lun_dataoutput, *) "" write(lun_dataoutput, "(a8)") "#EFIELDS" write(lun_dataoutput, *) nefield_files endif do ifile = 1, nefield_files if (dt_window == 0) window_dt = Efield_Weightings(ifile,htwindow_) if (OutputData) then iTempCounter = 0 do j = 1, npts_efield(ifile) DeltaT = abs(EfieldTimes(ifile, j) - currenttime) if (DeltaT <= window_dt) iTempCounter = iTempCounter + 1 enddo write(lun_dataoutput, "(a)") efield_file(ifile) write(lun_dataoutput, "(i4,1x,a1)") iTempCounter,efield_flag(ifile) endif if (efield_flag(ifile) == 'S') then do j = 1, npts_efield(ifile) DeltaT = abs(EfieldTimes(ifile, j) - currenttime) if (DeltaT <= window_dt) then NEField = NEField + 1 if (NEField > Max_Data) then write(*,*) "Too many Data points in efield_solver : EFields!" call stop_mpi("Yikes!") endif MagLat = EfieldLocations(ifile, j, mag_lat_) MagLT = EfieldLocations(ifile, j, mlt_) ! TimeErrorEfield = DeltaT / LinearDropOffTimeEfield * LinearDropOffEfield ! efield_wt(NEField) = & ! 1.0 / (EfieldDataError(ifile, j, 1)**2 + TimeErrorEfield**2) TotalError = ((1.+Efield_Weightings(ifile,errorper1_)) * & EfieldDataError(ifile, j, 1))**2 + & (Efield_Weightings(ifile,valueper1_) * & EfieldData(ifile, j, 1))**2 + & (DeltaT * Efield_Weightings(ifile,timeerror1_)/window_dt)**2 + & Efield_Weightings(ifile,minerror1_)**2 efield_wt(NEField) = Efield_Weightings(ifile,xwt1_) / TotalError efield_obs(NEField) = EfieldData(ifile, j, 1) if (OutputData) & write(lun_dataoutput, "(2f6.2,1p,4e10.2)") MagLat, MagLT, & EfieldData(ifile, j, 1), EfieldData(ifile, j, 2), & Efieldangle(ifile, j), DeltaT sp = sin(MagLT*htor) cp = cos(MagLT*htor) call fcmp(mmx,1,cp,sp,ff) Cos_se = cos(EfieldAngle(ifile, j)*dtor) Sin_se = sin(EfieldAngle(ifile, j)*dtor) Colat = (90.0-abs(MagLat))*dtor x = colat/dth ith = x ith = max(1,min(ithmx-2,ith)) x = x - float(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. sth = sin(colat) do m=-mmx,mmx fint = ff(m)/ri dfint = m*ff(-m)/sth/ri mm = iabs(m) do i=ibm(m),iem(m) ix = i - ibm(m) dpotdp = (xm1*q(2*ix+ns(mm),ith-1)+ & x0*q(2*ix+ns(mm),ith)+ & xp1*q(2*ix+ns(mm),ith+1)+ & xp2*q(2*ix+ns(mm),ith+2))*dfint dpotdt = (xm1*dq(2*ix+ns(mm),ith-1)+ & x0*dq(2*ix+ns(mm),ith)+ & xp1*dq(2*ix+ns(mm),ith+1)+ & xp2*dq(2*ix+ns(mm),ith+2))*fint ! CD = COS(ANGLE BETWEEN NORTH AND OBSERVATION DIRECTION, POS. TOWARDS ! EAST). SD = SIN OF THIS ANGLE. FACTOR 1.E3 CONVERTS V/M TO MV/M. efma(j,i) = (-Cos_se*dpotdp - Sin_se*dpotdt)*1.E3 datfma(NEField,i) = (Cos_se*dpotdt - Sin_se*dpotdp)*1.E3 enddo enddo endif enddo else do j = 1, npts_efield(ifile) DeltaT = abs(EfieldTimes(ifile, j) - currenttime) if (DeltaT <= window_dt) then NEField = NEField + 2 if (NEField > Max_Data) then write(*,*) "Too many Data points in efield_solver: EFields, 2 comps!" call stop_mpi("Yikes!") endif MagLat = EfieldLocations(ifile, j, mag_lat_) MagLT = EfieldLocations(ifile, j, mlt_) ! TimeErrorEfield = DeltaT / LinearDropOffTimeEfield * LinearDropOffEfield ! efield_wt(NEField-1) = & ! 1.0 / (EfieldDataError(ifile, j, 1)**2 + TimeErrorEfield**2) ! efield_wt(NEField) = & ! 1.0 / (EfieldDataError(ifile, j, 2)**2 + TimeErrorEfield**2) TotalError = ((1.+Efield_Weightings(ifile,errorper1_)) * & EfieldDataError(ifile, j, 1))**2 + & (Efield_Weightings(ifile,valueper1_) * & EfieldData(ifile, j, 1))**2 + & (DeltaT * Efield_Weightings(ifile,timeerror1_)/window_dt)**2 + & Efield_Weightings(ifile,minerror1_)**2 efield_wt(NEField-1) = Efield_Weightings(ifile,xwt1_) / TotalError efield_obs(NEField-1) = EfieldData(ifile, j, 1) TotalError = ((1.+Efield_Weightings(ifile,errorper2_)) * & EfieldDataError(ifile, j, 2))**2 + & (Efield_Weightings(ifile,valueper2_) * & EfieldData(ifile, j, 2))**2 + & (DeltaT * Efield_Weightings(ifile,timeerror2_)/window_dt)**2 + & Efield_Weightings(ifile,minerror2_)**2 efield_wt(NEField) = Efield_Weightings(ifile,xwt2_) / TotalError efield_obs(NEField) = EfieldData(ifile, j, 2) if (OutputData) & write(lun_dataoutput, "(2f6.2,1p,4e10.2)") MagLat, MagLT, & EfieldData(ifile, j, 1), EfieldData(ifile, j, 2), & Efieldangle(ifile, j), DeltaT sp = sin(MagLT*htor) cp = cos(MagLT*htor) call fcmp(mmx,1,cp,sp,ff) Cos_se = cos(EfieldAngle(ifile, j)*dtor) Sin_se = sin(EfieldAngle(ifile, j)*dtor) Colat = (90.0-abs(MagLat))*dtor x = colat/dth ith = x ith = max(1,min(ithmx-2,ith)) x = x - float(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. sth = sin(colat) do m=-mmx,mmx fint = ff(m)/ri dfint = m*ff(-m)/sth/ri mm = iabs(m) do i=ibm(m),iem(m) ix = i - ibm(m) dpotdp = (xm1*q(2*ix+ns(mm),ith-1)+ & x0*q(2*ix+ns(mm),ith)+ & xp1*q(2*ix+ns(mm),ith+1)+ & xp2*q(2*ix+ns(mm),ith+2))*dfint dpotdt = (xm1*dq(2*ix+ns(mm),ith-1)+ & x0*dq(2*ix+ns(mm),ith)+ & xp1*dq(2*ix+ns(mm),ith+1)+ & xp2*dq(2*ix+ns(mm),ith+2))*fint ! CD = COS(ANGLE BETWEEN NORTH AND OBSERVATION DIRECTION, POS. TOWARDS ! EAST). SD = SIN OF THIS ANGLE. FACTOR 1.E3 CONVERTS V/M TO MV/M. datfma(NEField-1,i) = (Cos_se*dpotdt - Sin_se*dpotdp)*1.E3 datfma(NEField,i) = (-Cos_se*dpotdp - Sin_se*dpotdt)*1.E3 enddo enddo endif enddo endif enddo if (OutputData .and. ncurrent_files > 0) then write(lun_dataoutput, *) "" write(lun_dataoutput, "(a8)") "#CURRENTS" write(lun_dataoutput, *) ncurrent_files endif do ifile = 1, ncurrent_files if (dt_window == 0) window_dt = Current_Weightings(ifile,htwindow_) if (OutputData) then iTempCounter = 0 do j = 1, npts_current(ifile) DeltaT = abs(CurrentTimes(ifile, j) - currenttime) if (DeltaT <= window_dt) iTempCounter = iTempCounter + 1 enddo write(lun_dataoutput, "(a)") current_file(ifile) write(lun_dataoutput, "(i4,1x,a2)") iTempCounter,current_flag(ifile) endif do j = 1, npts_current(ifile) DeltaT = abs(CurrentTimes(ifile, j) - currenttime) if (DeltaT <= window_dt) then NEField = NEField + 2 if (NEField > Max_Data) then write(*,*) "Too many Data points in efield_solver: Currents!" call stop_mpi("Yikes!") endif do k=1,imx ymg(k) = 0. zmg(k) = 0. datfma(NEField-1,k) = 0. datfma(NEField,k) = 0. enddo MagLat = CurrentLocations(ifile, j, mag_lat_) MagLT = CurrentLocations(ifile, j, mlt_) Cos_se = cos(CurrentAngle(ifile, j)*dtor) Sin_se = sin(CurrentAngle(ifile, j)*dtor) if (OutputData) & write(lun_dataoutput, "(2f6.2,1p,8e10.2)") MagLat, MagLT, & CurrentDalat(ifile, j), CurrentDalon(ifile, j), & CurrentData(ifile, j, ddby_), CurrentData(ifile, j, ddbz_), & CurrentData(ifile, j, dby_), CurrentData(ifile, j, dbz_), & CurrentAngle(ifile, j), DeltaT if (current_flag(ifile) == 'DD') then idel2 = 2 x = sin((90.-abs(MagLat))*dtor) dalat = CurrentDalat(ifile,j) dalon = CurrentDalon(ifile,j) dist = ri * sqrt(dalat**2 + x*x*dalon**2) / rad x = 1.E6/dist CurrentDalat(ifile,j) = dalat * x CurrentDalon(ifile,j) = dalon * x CurrentData(ifile, j, ddby_) = CurrentData(ifile, j, ddby_) * x CurrentData(ifile, j, ddbz_) = CurrentData(ifile, j, ddbz_) * x CurrentDataError(ifile, j, ddby_) = CurrentDataError(ifile, j, ddby_) * x CurrentDataError(ifile, j, ddbz_) = CurrentDataError(ifile, j, ddbz_) * x dclat = -.005*CurrentDalat(ifile, j) / rad dxlon = .005*CurrentDalon(ifile, j) / rad dxmlt = .005*CurrentDalon(ifile, j)*12./180. signyz = 100. else idel2 = 1 dclat = 0. dxlon = 0. dxmlt = 0. signyz = 1. x = 1. CurrentData(ifile, j, ddby_) = CurrentData(ifile, j, dby_) CurrentData(ifile, j, ddbz_) = CurrentData(ifile, j, dbz_) !Increase standard deviation by 3 times CurrentDataError(ifile, j, ddby_) = CurrentDataError(ifile, j, dby_) * 3. CurrentDataError(ifile, j, ddbz_) = CurrentDataError(ifile, j, dbz_) * 3. endif ! TimeErrorCurrent = DeltaT / LinearDropOffTimeCurrent * LinearDropOffCurrent + x ! efield_wt(NEField-1) = & ! 1.0 / (CurrentDataError(ifile, j, ddby_)**2 + TimeErrorCurrent**2) ! efield_wt(NEField) = & ! 1.0 / (CurrentDataError(ifile, j, ddbz_)**2 + TimeErrorCurrent**2) ! efield_obs(NEField-1) = CurrentData(ifile, j, ddby_) ! efield_obs(NEField) = CurrentData(ifile, j, ddbz_) TotalError = ((1.+Current_Weightings(ifile,errorper1_)) * & CurrentDataError(ifile, j, ddby_))**2 + & (Current_Weightings(ifile,valueper1_) * & CurrentData(ifile, j, ddby_))**2 + & (DeltaT * Current_Weightings(ifile,timeerror1_)/window_dt)**2 + & Current_Weightings(ifile,minerror1_)**2 efield_wt(NEField-1) = Current_Weightings(ifile,xwt1_) / TotalError TotalError = ((1.+Current_Weightings(ifile,errorper2_)) * & CurrentDataError(ifile, j, ddbz_))**2 + & (Current_Weightings(ifile,valueper2_) * & CurrentData(ifile, j, ddbz_))**2 + & (DeltaT * Current_Weightings(ifile,timeerror2_)/window_dt)**2 + & Current_Weightings(ifile,minerror2_)**2 efield_wt(NEField) = Current_Weightings(ifile,xwt2_) / TotalError efield_obs(NEField-1) = CurrentData(ifile, j, ddby_) efield_obs(NEField) = CurrentData(ifile, j, ddbz_) do idel=1,idel2 clatsap = (90. - abs(Maglat)) / rad + dclat sp = sin(MagLT*htor+dxlon) cp = cos(MagLT*htor+dxlon) call fcmp(mmx,1,cp,sp,ff) x = clatsap / dth ith = x ith = max0(1,min0(ithmx-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. sth = sin(clatsap) do m=-mmx,mmx fint = ff(m)/ri dfint = m*ff(-m)/sth/ri ixx = ns(iabs(m)) + 1 do i=ibm(m),iem(m) ix = i - ibm(m) dpotdp = (xm1*q(2*ix+ixx,ith-1)+x0*q(2*ix+ixx,ith) & +xp1*q(2*ix+ixx,ith+1)+xp2*q(2*ix+ixx,ith+2))*dfint dpotdt =(xm1*dq(2*ix+ixx,ith-1)+x0*dq(2*ix+ixx,ith) & +xp1*dq(2*ix+ixx,ith+1)+xp2*dq(2*ix+ixx,ith+2))*fint dmg(i) = (Cos_se*dpotdt - Sin_se*dpotdp)*200.*twopi hmg(i) = (Cos_se*dpotdp + Sin_se*dpotdt)*200.*twopi enddo enddo do I=1,IMX ymg(I) = ymg(I) + hmg(I)*signyz zmg(I) = zmg(I) + dmg(I)*signyz enddo dclat = .005*CurrentDalat(ifile, j)/RAD dxlon = -.005*CurrentDalon(ifile, j)/RAD dxmlt = .005*CurrentDalon(ifile, j)*12./180. signyz = -100. enddo do i=1,imx do k=1,imx datfma(NEField-1,k) = datfma(NEField-1,k) + ymg(i)*tfma(i,k) datfma(NEField,k) = datfma(NEField,k) + zmg(i)*tfma(i,k) enddo enddo endif enddo enddo ! COMPUTE "MODEL" OBSERVATIONS, OBSH, AND DIFFERENCES BETWEEN ACTUAL AND ! MODEL OBSERVATIONS, OBS. obsh(1:NEField) = 0.0 DO I=1,IMX OBSH(1:NEField) = OBSH(1:NEField) + AO(I)*DATFMA(1:NEField,I) enddo efield_obs(1:NEField) = efield_obs(1:NEField) - OBSH(1:NEField) ! COMPUTE EXPECTED RMS VARIATION IN E FIELD AND EXPECTED ERROR IN FIT. ! EVAR HOLDS VARIATION AND EFER HOLDS ERROR. ! 88/11/7 EEERR AND ENERR HOLD EAST AND NORTH COMPONENTS OF ERROR IN E. if (DebugLevel > 1) & write(*,*) 'Going into matsolv for e-fields, NEField = ',NEField call matsolv(1, NEField, datfma, efield_wt, xwt_dummy, efield_obs, & obsh, er, asum, ao, a, ah, diff_pot, diff_dummy, dcd_pot, & err_north, err_east, rms_pot, err_pot) ! CALCULATE COEFFICIENTS FOR CURRENT FUNCTION (PS), MAGNETIC ! POTENTIAL (V), and current potential (TF) a(1:imx) = asum(1:imx) ! Let's get the predicted values at the station locations: nE = 1 nBad = 0 nGood = 0 NewError = 0 do j=1,nmagnetometers ! IsFound = .false. ! if (mglat(j)*ihsoln < 0 .or. hgm(n,j) > 9998.0) IsFound = .true. ! do iDir = 1, 3 ! if (iMagErrorFlag(n,j,iDir) > iMaxErrorFlag) IsFound = .true. ! enddo ! ! if (.not. IsFound) then if (mglat(j)*ihsoln > 0) then MagPredicted(n,j,1) = hgm(n,j) + diff_pot(nE) MagPredicted(n,j,2) = dgm(n,j) + diff_pot(nE+1) MagError(1) = EField_wt(nE) * diff_pot(nE)**2 / & (1.0 - dcd_pot(nE)*EField_wt(nE)) MagError(2) = EField_wt(nE+1) * diff_pot(nE+1)**2 / & (1.0 - dcd_pot(nE+1)*EField_wt(nE+1)) if (use_zed) then MagPredicted(n,j,3) = zgm(n,j) + diff_pot(nE+2) MagError(3) = EField_wt(nE+2) * diff_pot(nE+2)**2 / & (1.0 - dcd_pot(nE+2)*EField_wt(nE+2)) nE = nE + 3 else MagError(3) = 0 nE = nE + 2 endif NewError = NewError + sum(MagError) if (MagError(1) > 0.0) then do iDir = 1, 3 iMagErrorFlag(n,j,iDir) = floor(sqrt(MagError(iDir))/2) if (iMagErrorFlag(n,j,iDir) > 10) iMagErrorFlag(n,j,iDir) = 10 enddo ! We are only going to consider bad stations to be stations ! that have either a large difference or more than 1 component ! not matching. Hence the > 1, instead of > 0 (which would ! mean any error causes the solution to be redone). IsFound = .false. do iDir = 1, 3 if (iMagErrorFlag(n,j,iDir) > iMaxErrorFlag) then if (DebugLevel > 1) & write(*,"(a,i5,i4,i2,5f6.1,i2)") & "==> Bad Magnetometer found : ", & n,j,iDir, MagPredicted(n,j,iDir), & hgm(n,j), dgm(n,j), zgm(n,j),& MagError(iDir), iMagErrorFlag(n,j,iDir) IsFound = .true. endif enddo if (IsFound) then nBad = nBad + 1 else nGood = nGood + 1 endif endif endif enddo if (DebugLevel > 2) & write(*,*) "===> NewError, OldError ", NewError, OldError if (NewError > OldError) then IsWorse = .true. if (DebugLevel > 1) then write(*,*) "==> Yikes! Fit is worse!!!" write(*,*) "==> Raising the bar... ", iMaxErrorFlag-1 endif do j=1,nmagnetometers do iDir = 1, 3 if (iMagErrorFlag(n,j,iDir) <= iMaxErrorFlag) & iMagErrorFlag(n,j,iDir) = 0 enddo enddo iMaxErrorFlag = iMaxErrorFlag-1 else if (nBad > 5) then if (DebugLevel > 1) then write(*,*) "==> Too many bad stations found : ",nBad write(*,*) "==> Lowering the bar... ", iMaxErrorFlag+1 endif iMaxErrorFlag = iMaxErrorFlag+1 endif endif if (nBad > 0 .and. DebugLevel > 0) & write(*,*) "----------- REPEATING INVERSION !! -----------" if (Indices(n,ae_) > 200.0) then ! Have to add back Dst, since we are going to recompute it and ! do the subtraction again (in the calc_dst routine. do j = 1, nmagnetometers if (wth(n,j) /= 0.0) & hgm(n,j) = hgm(n,j) + & Indices(n,dst_) * sin((90. - mglat(j)) * dtor) enddo endif call calc_dst(n,iMaxErrorFlag+1) nInversions = nInversions + 1 OldError = NewError if (nBad == 0) IsDone = .true. if (nInversions >= 5) IsDone = .true. if (IsDone .and. DebugLevel > 2) then write(*,*) "===> I am Done!! ", nBad, nInversions, IsWorse endif enddo call calc_ae(n) if (OutputData) then nGood = 0 do j=1,nmagnetometers if (mglat(j)*ihsoln > 0) nGood = nGood + 1 enddo write(lun_dataoutput, "(a)") "" write(lun_dataoutput, "(a14)") "#MAGNETOMETERS" write(lun_dataoutput, *) nGood do j=1,nmagnetometers if (mglat(j)*ihsoln > 0) then if (use_zed) then ztemp = zgm(n,J) ztempwt = MagPredicted(n,j,3) else ztemp = 1.0e5 ztempwt = 0.0 endif write(lun_dataoutput,"(2f6.2,6f8.1,3i3,a)") & mglat(j), amltmg(j), & hgm(n,j), MagPredicted(n,j,1), & dgm(n,j), MagPredicted(n,j,2), & ztemp, ztempwt, iMagErrorFlag(n,j,:), " "//mgname(j) endif enddo nBad = 0 do j=1,nmagnetometers if (sum(iMagErrorFlag(n,j,1:3)) > 0) nBad = nBad + 1 enddo write(lun_dataoutput, "(a)") "" write(lun_dataoutput, "(a)") "#MAGFLAGS" write(lun_dataoutput, *) nBad write(lun_dataoutput, *) iMaxErrorFlag do j=1,nmagnetometers if (sum(iMagErrorFlag(n,j,1:3)) > 0) then write(lun_dataoutput, "(a,3i3)") mgname(j), iMagErrorFlag(n,j,:) endif enddo endif Magnetic_Coef = 0.0 DO K=1,IMX Magnetic_Coef(1:imx,current_function_) = & Magnetic_Coef(1:imx,current_function_) + psfma(1:imx,k)*a(k) Magnetic_Coef(1:imx,current_potential_) = & Magnetic_Coef(1:imx,current_potential_) + tfma(1:imx,k)*a(k) Magnetic_Coef(1:imx,vertical_magnetic_) = & Magnetic_Coef(1:imx,vertical_magnetic_) + zfma(1:imx,k)*a(k) Magnetic_Coef(1:imx,magnetic_potential_) = & Magnetic_Coef(1:imx,magnetic_potential_) + vfma(1:imx,k)*a(k) enddo ! COMPUTE ELECTRIC POTENTIAL (EFPOT), ELECTRIC FIELD (EE,EN), CURRENT ! FUNCTION (PSI), MAGNETIC POTENTIAL (VPOT), AND current potential (CPOT). EField_Solution = 0.0 DO M=-MMX,MMX EField_Coef = 0.0 MM = IABS(M) IXX = NS(MM) + 1 DO I=IBM(M),IEM(M) IX = I - IBM(M) EField_Coef(0:ithmx,potential_) = & EField_Coef(0:ithmx,potential_) + A(I)*Q(2*IX+NS(MM),0:ITHMX) EField_Coef(0:ithmx,efield_north_) = & EField_Coef(0:ithmx,efield_north_) + & A(I)*DQ(2*IX+NS(MM),0:ITHMX)/RI EField_Coef(0:ithmx,efield_east_) = & EField_Coef(0:ithmx,efield_east_) + & A(I)*Q(2*IX+NS(MM),0:ITHMX)*M/RI/ST(0:ITHMX) EField_Coef(0:ithmx,current_potential_) = & EField_Coef(0:ithmx,current_potential_) + & Magnetic_Coef(i,current_potential_)*Q(2*IX+NS(MM),0:ITHMX) EField_Coef(0:ithmx,current_function_) = & EField_Coef(0:ithmx,current_function_) + & Magnetic_Coef(i,current_function_)*Q(2*IX+IXX,0:ITHMX) EField_Coef(0:ithmx,vertical_magnetic_) = & EField_Coef(0:ithmx,vertical_magnetic_) + & Magnetic_Coef(i,vertical_magnetic_)*Q(2*IX+IXX,0:ITHMX) EField_Coef(0:ithmx,magnetic_potential_) = & EField_Coef(0:ithmx,magnetic_potential_) + & Magnetic_Coef(i,magnetic_potential_)*Q(2*IX+IXX,0:ITHMX) enddo IF (M < 0) then EField_Coef(0,efield_east_) = EField_Coef(0,efield_north_) else EField_Coef(0,efield_east_) = -EField_Coef(0,efield_north_) endif ! SUBTRACT MEAN VALUE OF POTENTIAL AT ITHPLT (SAME AS AT ITHMX) ! SUBTRACT MEAN VALUE OF POTENTIAL AT ITHTRNS (4/94) IF (M == 0) THEN EField_Coef(0:ithmx,potential_) = & EField_Coef(0:ithmx,potential_) - & EField_Coef(ithtrns,potential_) ENDIF ! I am not sure what to do about this.... ! ! IF (MM == 1) then ! IX = (M + 3)/2 ! pot1(ix,0:ithmx) = POTM(0:ithmx) ! endif DO ITH=0,ITHMX Efield_Solution(1:lonmx,ith,potential_) = & Efield_Solution(1:lonmx,ith,potential_) + & EField_Coef(ith,potential_)*amie_f(M,1:lonmx) Efield_Solution(1:lonmx,ith,magnetic_potential_) = & Efield_Solution(1:lonmx,ith,magnetic_potential_) + & EField_Coef(ith,magnetic_potential_)*amie_f(M,1:lonmx) Efield_Solution(1:lonmx,ith,vertical_magnetic_) = & Efield_Solution(1:lonmx,ith,vertical_magnetic_) + & EField_Coef(ith,vertical_magnetic_)*amie_f(M,1:lonmx) Efield_Solution(1:lonmx,ith,current_potential_) = & Efield_Solution(1:lonmx,ith,current_potential_) + & EField_Coef(ith,current_potential_)*amie_f(M,1:lonmx) Efield_Solution(1:lonmx,ith,current_function_) = & Efield_Solution(1:lonmx,ith,current_function_) + & EField_Coef(ith,current_function_)*amie_f(M,1:lonmx) Efield_Solution(1:lonmx,ith,efield_north_) = & Efield_Solution(1:lonmx,ith,efield_north_) + & EField_Coef(ith,efield_north_)*amie_f(M,1:lonmx) Efield_Solution(0,ith,efield_north_) = & Efield_Solution(lonmx,ith,efield_north_) Efield_Solution(1:lonmx,ith,efield_east_) = & Efield_Solution(1:lonmx,ith,efield_east_) + & EField_Coef(ith,efield_east_)*amie_f(-M,1:lonmx) Efield_Solution(0,ith,efield_east_) = & Efield_Solution(lonmx,ith,efield_east_) enddo enddo Efield_Solution(0,:,potential_) = Efield_Solution(lonmx,:,potential_) Efield_Solution(0,:,current_potential_) = & Efield_Solution(lonmx,:,current_potential_) Efield_Solution(0,:,current_function_) = & Efield_Solution(lonmx,:,current_function_) Efield_Solution(0,:,vertical_magnetic_) = & Efield_Solution(lonmx,:,vertical_magnetic_) Efield_Solution(0,:,magnetic_potential_) = & Efield_Solution(lonmx,:,magnetic_potential_) Efield_Solution(:,:,current_east_) = & Calculated_Total_Pedersen(:,:)*Efield_Solution(:,:,efield_east_) + & Calculated_Total_Hall(:,:) *Efield_Solution(:,:,efield_north_) Efield_Solution(:,:,current_north_) = & Calculated_Total_Pedersen(:,:)*Efield_Solution(:,:,efield_north_) - & Calculated_Total_Hall(:,:) *Efield_Solution(:,:,efield_east_) Efield_Solution(:,:,current_horizontal_) = & sqrt(Efield_Solution(:,:,current_east_)**2 + & Efield_Solution(:,:,current_north_)**2) JouleHeating(:,:) = Calculated_Total_Pedersen(:,:) * & (Efield_Solution(:,:,current_east_)**2 + & Efield_Solution(:,:,current_north_)**2) end subroutine Solve_for_EField subroutine EField_Cu_set use ModQset use ModSize use ModAMIE use ModEField use ModIndices, only : Indices, Kp_ use ModConstants implicit none ! ANISO = factor affecting auroral zone EN/EE anisotropy ! XK1 = relative importance of (E/(div)**(EXNNP1/2)E)**2 suppression ! XK2 = factor affecting relative polar-cap electric-field strength ! (with respect to Heppner-Maynard LT boundary) ! XK5 = factor affecting relative polar-cap electric-field strength ! (with respect to Heppner-Maynard midnight boundary) ! XK3 = relative importance of (grad.divE/(div)**EXNNP1/2)E)**2 suppress ! XK4 = minimum value of WT ! F1 = fraction of Heppner-Maynard midnight low-latitude radius relating ! to anisotropy maximum ! F2 = fraction of Heppner-Maynard LT low-latitude radius relating to ! maximum of electric field strength ! EXPO = exponent on colatitude for auroral-zone weighting functions ! EXNNP1 = exponent on n(n+1) and m*m real, parameter :: aniso = 30. real, parameter :: xk1 = 100. real, parameter :: xk3 = 0. real, parameter :: xk2 = .2 ! real, parameter :: xk2 = .02 real, parameter :: xk5 = .0 real, parameter :: xk4 = .004 real, parameter :: f1 = 1.1 real, parameter :: f2 = .75 real, parameter :: expoe = 14. real, parameter :: expoan = 11. real, parameter :: exnnp1 = 1.25 real :: xko2, gam2, ri_norm, area, kpfac, kp, Midnight_Boundary real :: th1, th12, th1n, th2, th2n, wt1, wt2, f0, sp, cp, sth real :: fint, dfint, x, xm2, xx, dpotdp, dpotdt, xmult real, dimension(ithmx) :: wt real, dimension(24) :: shield_lat real, dimension(ithtrns) :: thne, thnan real, dimension(nq) :: xnnp1, xn2x real, dimension(mcoef) :: d2ee, divee, xlappot, xlapee, xlapes, z ! real, dimension(0:mmx) :: m2x integer :: i, j, k, n, m, mm, ix, ith, lon, imlt integer, dimension(mcoef) :: ipvt ! DIMENSION XNNP1(NQ),THNE(ITHTRNS),XLAPPOT(MCOEF), ! | LAP2POT(MCOEF),DIVEE(MCOEF),D3EE(MCOEF), ! | XLAPEE(MCOEF),XLAPES(MCOEF),D2EE(MCOEF), ! | THNAN(ITHTRNS),XN2X(NQ),M2X(0:MMX) ! DIMENSION IPVT(MCOEF),Z(MCOEF) real*8 :: cu_double(MCOEF,MCOEF) logical :: IsPhysical ! SHIELD_LAT is the low latitude convection boundary for Kp=3+ from 1-24 ! MLT from Heppner and Maynard (1987, JGR, Fig 10, p 4482) DATA shield_lat/58.4,58.4,58.5,58.6,58.7,59.3,60.7,63.8,67.2,69.1,& 70.,69.6,68.7,67.3,64.5,62.,60.,59.4,59.3,59.2,58.7,& 58.5,58.4,58.4/ kp = Indices(iteration_number, kp_) if (kp < 0.52) kp = 0.52 ! FACTORS FOR CU FORMULA xko2 = ri/3.e5 xko2 = xko2*xko2 gam2 = 10. xnnp1(1:nq) = qnorm(1:nq)*qnorm(1:nq)/(1.-ct(ithtrns)) xn2x(1:nq) = xnnp1(1:nq)**exnnp1 ! do m=1,mmx ! m2x(m) = (m*m)**exnnp1 ! enddo ! m2x(0) = 0. do ith=1,ithtrns thne(ith) = (ith*dth)**expoe thnan(ith) = (ith*dth)**expoan enddo ri_norm = 1. cu = 0.0 cu_double = 0.0 do j=1,mcoef ax(:,j+imx) = 0. enddo area = 0. kpfac = kp*2.5 - 7.5 !\ ! This insures some solution !/ if (kp <= 0.0) kpfac = 3.0*2.5 - 8.3 if (kp > 2.and.kp < 3.) kpfac = (kp-2.)*1.7 -2.5 if (kp >= 3.) kpfac = kp*2.5 - 8.3 Midnight_Boundary = (90. - shield_lat(24) + kpfac)/rad th1 = f1*Midnight_Boundary th12 = th1*th1 th1n = th1**expoe th2 = f2*Midnight_Boundary th2n = th2**expoe wt1 = 2.*th2n*th1n/(th2n*th2n+th1n*th1n)+ & 2.*xk2*th2n/(th2n+th1n)+ & 2.*xk5*th1n/(th1n+th1n) wt1 = 1./wt1 + xk4 f0 = aniso*wt1 do lon=1,lonmx sp = sin(twopi*float(lon)/float(lonmx)) cp = cos(twopi*float(lon)/float(lonmx)) call fcmp(mmx,1,cp,sp,ff) imlt = 24.*float(lon)/float(lonmx) + .5 if (imlt == 0) imlt = 24 th2 = f2 * (90. - shield_lat(imlt) + kpfac)/rad th2n = th2**expoe ! th = (85.5 - shield_lat(imlt) + kpfac)/rad ! create a belt starting 3 degrees poleward of the convection boundary ! and going 31.5 degrees where the electric field is forced to be ! smaller and smaller. (belt need to extend to the basis function ! transition latitude.) ! do 360 lats=1,21 ! th = th + 1.5/rad do ith=1,ithtrns ! th = ith*dth wt1 = 2.*th2n*thne(ith)/(th2n*th2n+thne(ith)*thne(ith)) + & 2.*xk2*th2n/(th2n+thne(ith)) + & 2.*xk5*th1n/(th1n+thne(ith)) ! latitudes. changed in 9/91. ! if (lats.le.6) wtlat = exp((float(lats)-7.0)/2.0) wt1 = 1./wt1 + xk4 wt2 = 2.*th1n*thnan(ith)/(th1n*th1n+thnan(ith)*thnan(ith))*f0 sth = st(ith) do m=-mmx,mmx fint = ff(m)/ri_norm dfint = m*ff(-m)/sth/ri_norm mm = iabs(m) x = m/sth/ri_norm xm2 = 0. if (m.ne.0) xm2 = (x*x)**exnnp1*ff(m) ! assume that there are equal numbers of symmetric and antisymmetric ! qs for a given value of m. ns(mm) gives the first index value ! for the first subscript of q for wavenumber mm. do i=ibm(m),iem(m) ix = i - ibm(m) xx = xn2x(2*ix+ns(mm))/ri_norm/ri_norm dpotdp = q(2*ix+ns(mm),ith)*dfint dpotdt =dq(2*ix+ns(mm),ith)*fint xlappot(i) = xx*q(2*ix+ns(mm),ith)*fint*ri_norm ! lap2pot(i) = xx*xlappot(i) divee(i) = q(2*ix+ns(mm),ith)*xm2 ! d3ee(i) = divee(i)*x*x d2ee(i) = divee(i)*x qs(i) = dpotdp dqs(i) = dpotdt xlapee(i) = xx*qs(i) xlapes(i) = xx*dqs(i) enddo enddo if (ith.eq.ithtrns) sth=sth*.5 do j=1,imx cu_double(1:j,j) = cu_double(1:j,j) + & sth*(wt1*(xk3*(xlapee(1:j)*xlapee(j) + & xlapes(1:j)*xlapes(j)) + xlappot(1:j)*xlappot(j)+xk1* & (qs(1:j)*qs(j)+dqs(1:j)*dqs(j)))+& wt2*(xk3*d2ee(1:j)*d2ee(j)+divee(1:j)*divee(j)+& xk1*qs(1:j)*qs(j)) ) enddo area = area + sth enddo enddo xmult = 1./area ! fill upper-right part of cu with cui do j=1,imx cu_double(1:j,j) = cu_double(1:j,j)*xmult enddo do m=-mmx,mmx mm = iabs(m) do k=ibm(m),iem(m) cuid(k) = cu_double(k,k) do i=1,k ax(i,k) = cu_double(i,k) ax(k,i) = ax(i,k) enddo ax(k,k+imx) = 1. enddo enddo call solver(ax,imx,2*imx,mcoef,ipvt,z) cu = cu_double do i=1,imx cu(i,1:i) = .5*(ax(i,1+imx:i+imx) + ax(1:i,i+imx)) enddo ! analysis of idm oct 1981 data showed a power spectrum of 1/k**2 ! test whether cu is physical if (DebugLevel > 2) then do i=1,imx do j=1,i-1 if (cu(i,j)*cu(i,j).gt.cu(i,i)*cu(j,j)) then write(6,3025) i,j,cu(i,j),cu(i,i),cu(j,j) endif enddo enddo else IsPhysical = .true. do i=1,imx do j=1,i-1 if (cu(i,j)*cu(i,j).gt.cu(i,i)*cu(j,j)) IsPhysical = .false. enddo enddo if (.not.IsPhysical) & write(*,*) ' Warning : E-Field CU matrix has elements '// & 'which are non-physical' endif 3025 format(1x,'i=',i4,'j=',i4,'cu(i,j)=',e10.3,'cu(i,i)=',e10.3,& 'cu(j,j)=',e10.3) return end