! !------------------------------------------------------------------------ ! ! Subroutine for - reading in the E-field data from a file ! !------------------------------------------------------------------------ subroutine read_efields use ModAMIE use ModFiles use ModData use ModConductance use ModIndices implicit none integer :: iLunEfield, iyr, imon, iday, ihr, minute, isec, ierr integer :: nlines_to_read, i, j, n, k, ll, l, ifile, npt, iNC real*8 :: time_start, time, dt_magfile, time_before_all, ctime real :: alon, alat, hp character (len=100) :: line integer, dimension(7) :: itime iLunEfield = lun_efield npts_efield = 0 do ifile = 1, nefield_files if (DebugLevel > 2) write(*,*) "efield file --->", & efield_file(ifile),"<---" if (index(efield_file(ifile),efield_dir(1:2)) /= 1) then call merge_str(efield_dir, efield_file(ifile)) endif iLunEfield = iLunEfield + ifile - 1 open(iLunEfield, file = efield_file(ifile), status='old', iostat = ierr) if (ierr /= 0) then write(6,*) 'Error opening E-field file : -->',efield_file(ifile),'<--' endif ! We want to determine whether this file is a SuperDARN file or ! an amiein file. One difference between the files is that the ! SuperDARN data just starts right off with data, while the ! amiein file has a header. So, lets try reading some data to ! so if we get an error: read(iLunEfield,*, iostat = ierr) itime(1:6) if (ierr == 0) then ! file is a SuperDARN file rewind(iLunEfield) write(*,*) "I am opening SuperDarn file!" call read_superdarn(ifile, iLunEfield) write(*,*) "Over calling read_superdarn subrouting." else ! file is an amiein file rewind(iLunEfield) ! Now lets test to see if it is a DMSP data file from Marc read(iLunEfield,*) line if (index(line,"f") == 1) then if (DebugLevel > 5) write(*,*) "======> DMSP File!!" rewind(iLunEfield) call read_dmsp(ifile, iLunEfield) else rewind(iLunEfield) read(iLunEfield,*) line read(iLunEfield,*) line read(iLunEfield,*) line ctime = currenttime itime(7) = 0 if (DebugLevel > 2) & write(*,*) " ===> read_efields, n = ", ifile, ntimes time = 0.0 ctime = ctime - window_dt/2.0 if (DebugLevel > 2) & write(*,*) "ctime, window_dt = ", ctime, window_dt do while (time < ctime) read(iLunEfield,"(i4,5i3)", iostat = ierr) itime(1:6) if (ierr /= 0) then time = 1.0e32 else call time_int_to_real(itime, time) endif enddo if (time < 1.0e32) backspace(iLunEfield) if (DebugLevel > 2) & write(*,*) "time after while loop : ",time npts_efield(ifile) = 1 if (timelist(ntimes) > -1.0e32) then ctime = timelist(ntimes) else ctime = currenttime + dt * (ntimes-1) endif ctime = ctime + window_dt/2.0 do while (time < ctime) itime = 0 iNC = npts_efield(ifile) read(iLunEfield,"(i4,5i3,f6.2,f7.2,f6.2,i5,5f7.2)", & iostat = ierr) itime(1:6), & EfieldLocations(ifile, iNC, mag_lat_), & alon, EfieldLocations(ifile, iNC, mlt_), & npt, EfieldAngle(ifile,iNC), & EfieldData(ifile, iNC, 1), EfieldDataError(ifile, iNC, 1), & EfieldData(ifile, iNC, 2), EfieldDataError(ifile, iNC, 2) !\ ! There are 2 possible reasons why ierr /= 0: ! (1) EOF was reached. ! (2) a data value was missing from the row. ! Therefore, we have to check to see which is true. !/ !\ ! This is option 2 above !/ if (ierr /= 0 .and. itime(1) /= 0) ierr = 0 if (ierr /= 0) then !\ ! This is option 1 above !/ time = ctime + 1.0 else alat = EfieldLocations(ifile, iNC, mag_lat_) if (efield_flag(ifile) == 'S') then EfieldData(ifile, iNC, 2) = -999. EfieldDataError(ifile, iNC, 2) = -999. endif if (ihsoln*alat > 0) then call time_int_to_real(itime, time) EfieldTimes(ifile, iNC) = time ! ! this implies iNC = iNC + 1 ! npts_efield(ifile) = npts_efield(ifile) + 1 endif endif enddo endif endif close(iLunEfield) npts_efield(ifile) = npts_efield(ifile) - 1 if (DebugLevel > 1) write(*,*) " ==> ",npts_efield(ifile), & " points read from file ",efield_file(ifile) enddo end subroutine read_efields !------------------------------------------------------------------------ ! ! Subroutine for - reading in the E-field data from a file ! !------------------------------------------------------------------------ subroutine read_superdarn(ifile, iLunEfield) use ModAMIE use ModFiles use ModData use ModConductance use ModIndices implicit none integer, external :: jday integer, intent(in) :: ifile, iLunEfield integer :: iYear, iJulianDay, ihr, minute, isec, ierr integer :: nlines_to_read, i, j, n, k, ll, l, npt, iNC real*8 :: time_start, time, dt_magfile, time_before_all, ctime real :: alon, alat, amlt, rtime, glat,glon real :: v, zmag, ymag, xmag, alontest, alattest, a, bmag real, dimension(20) :: data integer :: nRadars, nParams, nPts character (len=100) :: line integer, dimension(7) :: itime logical :: IsVerbose if (timelist(1) > -1.0e32) then ctime = timelist(1) else ctime = currenttime endif call time_real_to_int(ctime, itime) iYear = itime(1) iJulianDay = jday(iYear, itime(2), itime(3)) itime(7) = 0 time = 0.0 ctime = ctime - window_dt/2.0 do while (time < ctime) read(iLunEfield,*, iostat = ierr) itime(1:6) if (ierr /= 0) then time = 1.0e32 else call time_int_to_real(itime, time) if (time < ctime) then ! This line is a "2" read(iLunEfield, *) line read(iLunEfield, *) nRadars, nParams ! Three header lines read(iLunEfield, *) line read(iLunEfield, *) line read(iLunEfield, *) line do i=1,nRadars read(iLunEfield, *) line enddo ! Three header lines read(iLunEfield, *) nPts, nParams ! Three MORE header lines read(iLunEfield, *) line read(iLunEfield, *) line read(iLunEfield, *) line do i=1,nPts read(iLunEfield,*) data(1:nParams) enddo endif endif enddo if (time < 1.0e32) backspace(iLunEfield) npts_efield(ifile) = 1 ! if (timelist(ntimes) > -1.0e32) then ! ctime = timelist(ntimes) ! else ! ctime = currenttime + dt * (ntimes-1) ! endif call time_real_to_int(ctime, itime) ctime = ctime + window_dt do while (time < ctime) itime = 0 IsVerbose = .false. iNC = npts_efield(ifile) read(iLunEfield,*, iostat = ierr) itime(1:6) write(*,*) "=> Reading SuperDARN ",itime call time_int_to_real(itime, time) ! if (itime(4) == 11 .and. itime(5) == 20) IsVerbose=.true. if (time < ctime .and. ierr == 0) then if (DebugLevel > 1) & write(*,*) " ==> Reading from SuperDARN file, time = ", itime call find_subsolar_info(iYear,iJulianDay,itime(4),itime(5),itime(6)) ! This line is a "2" read(iLunEfield, *) line read(iLunEfield, *) nRadars, nParams ! Three header lines read(iLunEfield, *) line read(iLunEfield, *) line read(iLunEfield, *) line do i=1,nRadars read(iLunEfield, *) line enddo read(iLunEfield, *) nPts, nParams ! Three MORE header lines read(iLunEfield, *) line read(iLunEfield, *) line read(iLunEfield, *) line i = 1 do while (i <= nPts) read(iLunEfield,*) data(1:nParams) alat = data(2) if (ihsoln*alat > 0) then alon = data(1) call magloctm(alon, subsolar_lat, subsolar_lon, & pole_colat, pole_east_lon, amlt) rtime = float(iyear) + float(ijulianday)/jday(iyear,12,31) call find_geographic_point(rtime, alat, alon, glat, glon, IsVerbose) call apex(rtime, glat, glon, 110.0, a, & alattest, alontest, bmag, xmag, ymag, zmag, v) EfieldTimes(ifile, iNC) = time EfieldLocations(ifile, iNC, mag_lat_) = alat EfieldLocations(ifile, iNC, mlt_) = amlt ! EfieldAngle(ifile,iNC) = data(4) + 90.0 EfieldAngle(ifile,iNC) = data(3) + 90.0 ! EfieldData(ifile, iNC, 1) = data(3)*abs(zmag)*1.0e-9*1.0e3 EfieldData(ifile, iNC, 1) = data(7)*abs(zmag)*1.0e-9*1.0e3 ! EfieldDataError(ifile, iNC, 1) = data(7)*abs(zmag)*1.0e-9*1.0e3 EfieldDataError(ifile, iNC, 1) = data(8)*abs(zmag)*1.0e-9*1.0e3 EfieldData(ifile, iNC, 2) = 0.0 EfieldDataError(ifile, iNC, 2) = -999.0 iNC = iNC + 1 if (iNC >= Max_Data) then if (DebugLevel > 1) & write(*,*) " ==> Too many data points from this file!!" ! call stop_mpi("Tell Aaron to fix this problem!!") iNC = Max_Data + 2 i = nPts time = ctime + 1.0 endif endif i = i + 1 enddo npts_efield(ifile) = iNC else time = ctime + 1.0 endif enddo if (npts_efield(ifile) > Max_Data) npts_efield(ifile) = Max_Data write(*,*) "finish reading SuperDARN data!!" end subroutine read_superdarn !------------------------------------------------------------------------ ! ! Subroutine for - reading in the E-field data from a DMSP data file ! !------------------------------------------------------------------------ subroutine read_dmsp(ifile, iLunEfield) use ModAMIE use ModFiles use ModData implicit none integer, external :: jday integer, intent(in) :: ifile, iLunEfield integer :: iYear, iJulianDay, iError character (len=100) :: line integer, dimension(7) :: itime real :: rDate, rTime real*8 :: cTime, Time integer :: iQuality1, iQuality2, iNC real :: altitude, glat, glon, mlat, amlt, vx, vy, vz, rmsx, rmsy, rmsz real :: alattest, alontest, xmag, ymag, zmag, v, bmag, a real :: old_mlt, old_mlat, dn, de, x,y logical :: flag1,flag2 iNC = 1 if (timelist(1) > -1.0e32) then ctime = timelist(1) else ctime = currenttime endif call time_real_to_int(ctime, itime) if (DebugLevel > 2) write(*,*) "===> itime : ",itime itime(7) = 0 time = 0.0 ctime = ctime - window_dt/2.0 read(iLunEfield, *) line read(iLunEfield, *) line read(iLunEfield, *) line do while (time < ctime) read(iLunEfield,101,iostat = iError) & rDate, rTime itime(1) = rDate/1000 if (itime(1) >= 100) iTime(1) = iTime(1) - 100 itime(2) = 1 itime(3) = mod(int(rDate),1000) itime(4) = 0 itime(5) = 0 itime(6) = rTime call time_int_to_real(itime, time) call time_real_to_int(time, itime) if (DebugLevel > 2) write(*,*)'===>itime from DMSP : ', itime(1:6) if (iError /= 0) then ! This may be another orbit stuck together. backspace(iLunEfield) ! Now lets test to see if it is a DMSP data file from Marc read(iLunEfield,*) line if (index(line,"f") == 1) then write(*,*) "Another DMSP Orbit was found!!" read(iLunEfield, *) line read(iLunEfield, *) line else time = 1.0e32 endif endif enddo if (time < 1.0e32) then ! backspace(iLunEfield) backspace(iLunEfield) ! We need to do this to make sure that we can calculate an ! angle, which takes at least two points... if (DebugLevel > 2) write(*,*)'===> Begin to read DMSP data: ' read(iLunEfield,101, iostat = iError) rDate, rTime, & iQuality1, iQuality2, Altitude, & glat, glon, old_mlat, old_mlt, & vx, vy, vz, rmsx, rmsy, rmsz 101 FORMAT(f10.0,f8.1,2i2,f7.1,f8.2,f8.2,f8.2,f8.2,3f8.1,f8.2,& 2f8.1,1e15.4,3f9.2,2f7.0,i7) endif npts_efield(ifile) = 1 ctime = ctime + window_dt do while (time < ctime) read(iLunEfield,101, iostat = iError) rDate, rTime, & iQuality1, iQuality2, Altitude, & glat, glon, mlat, amlt, & vx, vy, vz, rmsx, rmsy, rmsz itime(1) = rDate/1000 if (itime(1) >= 100) itime(1) = itime(1) - 100 itime(2) = 1 itime(3) = mod(int(rDate),1000) itime(4) = 0 itime(5) = 0 itime(6) = rtime if (DebugLevel > 4) write(*,*)'=====> time from DMSP:', itime(1:6) call time_int_to_real(itime, time) flag1=(vx > -8000.0) .and. (vy > -8000.0) .and. (vz > -8000.0) flag2=(iQuality1 < 3.0) .and. (iQuality2 <3.0) if (mlat * ihsoln > 45.0 .and. flag1 .and. flag2 ) then rtime = float(itime(1)) + float(itime(3))/jday(iyear,12,31) if (rtime < 65) rtime = rtime+2000 if (rtime < 1900) rtime = rtime+1900 call apex(rtime, glat, glon, 110.0, a, & alattest, alontest, bmag, xmag, ymag, zmag, v) EfieldTimes(ifile, iNC) = time EfieldLocations(ifile, iNC, mag_lat_) = mlat EfieldLocations(ifile, iNC, mlt_) = amlt de = (amlt-old_mlt)*sin(mlat*3.14159/180.0)*3.14159/24.0 dn = (mlat-old_mlat)*3.14159/180.0 x=min(1.0,dn/sqrt(de*de+dn*dn)) y = asin(min(1.0,x))*180.0/3.14159 EfieldAngle(ifile,iNC) = y EfieldData(ifile, iNC, 1) = vy*abs(zmag)*1.0e-9*1.0e3 EfieldDataError(ifile, iNC, 1) = rmsy*abs(zmag)*1.0e-9*1.0e3 EfieldData(ifile, iNC, 2) = 0.0 EfieldDataError(ifile, iNC, 2) = -999.0 iNC = iNC + 1 if (iNC >= Max_Data) then if (DebugLevel > 1) & write(*,*) " ==> Too many data points from this file!!" call stop_mpi("Tell Aaron to fix this problem!!") endif endif old_mlat = mlat old_mlt = amlt ! We don't need to check for more orbit files, since the starting ! latitude is going to be 0 degrees, which is completely outside ! of the AMIE area. If we do past this, then we are done anyways... if (iError /= 0) time = 1.0e32 enddo npts_efield(ifile) = iNC if (DebugLevel > 2) write(*,*) "===> Finish reading DMSP data : ", iNC end subroutine read_dmsp