!------------------------------------------------------------------------ ! ! Subroutine for reading in the master_psi file for magnetometers ! !------------------------------------------------------------------------ subroutine read_psi_file use ModAMIE use ModFiles use ModData implicit none integer :: iunit, ierr, n character (len=100) :: line real :: glat, glon, decst, psi iunit = lun_magnetometer open(iunit, file = psi_file, status='old', iostat = ierr) if (ierr /= 0) then write(6,*) 'Error opening file :',psi_file stop endif read(iunit,*) line read(iunit,*) line ierr = 0 nmagnetometers = 1 do while (ierr == 0) read(iunit,'(4X,A3,13X,2F7.2,F8.2,F7.2,2X,3F7.2)', & iostat=ierr) mgname(nmagnetometers),glat,glon, & mglat(nmagnetometers),mglon(nmagnetometers),decst,psi, & weights(nmagnetometers) decmg(nmagnetometers) = 0. if (ierr == 0) nmagnetometers = nmagnetometers + 1 enddo close(iunit) end subroutine read_psi_file !------------------------------------------------------------------------ ! ! Subroutine for - reading in the magnetometer data from a file ! - calculating AE ! - removing Dst from magnetometers ! !------------------------------------------------------------------------ subroutine read_magnetometers use ModAMIE use ModFiles use ModData use ModIndices use ModConstants implicit none integer :: iunit, iyr, imon, iday, ihr, minute, isec, ierr integer :: nlines_to_read, i, j, n, k, ll, l real*8 :: time_start, time, dt_magfile, time_before_all, ctime real, dimension(MAX_MAGNETOMETERS) :: xmg,ymg,zmg,xer,yer,zer character (len=3), dimension(MAX_MAGNETOMETERS) :: stname real, dimension(MAX_MAGNETOMETERS) :: dstval real, dimension(0:MAX_MAGNETOMETERS) :: dstlon real, dimension(3) :: MagMean, MagStd integer, dimension(3) :: MagNPts logical :: IsBad integer, dimension(7) :: itime character (len=3) :: name, xyz character (len=1) :: dummy real :: snrsq, dx, erzgm, erhgm, erdgm, x, clat iunit = lun_magnetometer open(iunit, file = magnetometer_file, status='old', iostat = ierr) if (ierr /= 0) then write(6,*) 'Error opening file :',magnetometer_file endif !\ ! Figure out how many magnetometers there are by reading in all of the ! lines which have the same time !/ nmagnetometers = 0 read (iunit,"(3i2,1x,2i2)") iyr,imon,iday,ihr,minute itime(1) = iyr itime(2) = imon itime(3) = iday itime(4) = ihr itime(5) = minute itime(6) = 0 itime(7) = 0 call time_int_to_real(itime, time_start) time = time_start do while (time == time_start) read (iunit,"(3i2,1x,2i2)") iyr,imon,iday,ihr,minute itime(1) = iyr itime(2) = imon itime(3) = iday itime(4) = ihr itime(5) = minute call time_int_to_real(itime, time) nmagnetometers = nmagnetometers + 1 enddo dt_magfile = time - time_start nlines_to_read = (currenttime - time_start) / dt_magfile * nmagnetometers rewind(iunit) do i = 1, nlines_to_read read (iunit,'(a)') dummy enddo Indices(:,au_) = 0.0 Indices(:,al_) = 0.0 Indices(:,ae_) = 0.0 Indices(:,dst_) = 0.0 nae = 0 ndst = 0 if (timelist(1) > -1.0e32) then ctime = timelist(1) else ctime = currenttime endif iMagErrorFlag = 0 do n = 1, ntimes if (DebugLevel > 2) write(*,*) " ===> read_magnetometers, n = ", n, ntimes do j = 1, nmagnetometers read (iunit,"(3i2,1x,2i2,1x,a3,4x,6f6.0)") & iyr,imon,iday,ihr,minute,stname(j),xmg(j), & ymg(j),zmg(j),xer(j),yer(j),zer(j) itime(1) = iyr itime(2) = imon itime(3) = iday itime(4) = ihr itime(5) = minute call time_int_to_real(itime, time) if (time /= ctime) then write(*,*) "Error in reading in magnetometer data." write(*,*) "Expected time : ", currenttime write(*,*) "Recieved time : ", time stop endif if (stname(j) .ne. mgname(j)) THEN write(6,*) 'station mismatch!!!' write(6,*) 'read : ',stname(j) write(6,*) 'expected : ',mgname(j) stop endif ierr = 0 if (abs(xmg(j)).gt.5000) ierr = 1 if (abs(ymg(j)).gt.5000) ierr = 1 if (abs(zmg(j)).gt.5000) ierr = 1 if (abs(xer(j)).gt.5000) ierr = 1 if (abs(yer(j)).gt.5000) ierr = 1 if (abs(zer(j)).gt.5000) ierr = 1 if (weights(j).eq.0) ierr = 1 if (ierr.ne.0) then xmg(j) = 99999.0 ymg(j) = 99999.0 zmg(j) = 99999.0 xer(j) = 0.0 yer(j) = 0.0 zer(j) = 0.0 endif ! Change sign of D (or Y) and Z in the Southern Hemisphere (since is ! solved as if it were in the North). (9/91) if (mglat(j) .lt. 0.) then if (ymg(j) .lt. 9998.) ymg(j) = -ymg(j) if (zmg(j) .lt. 9998.) zmg(j) = -zmg(j) endif zgm(n,j) = zmg(j) snrsq=zer(j)**2 erzgm = sqrt(snrsq) hgm(n,j) = xmg(j) snrsq=xer(j)**2 erhgm = sqrt(snrsq) dgm(n,j) = ymg(j) snrsq=yer(j)**2 erdgm = sqrt(snrsq) enddo MagMean = 0.0 MagStd = 0.0 MagNPts = 0 do j = 1, nmagnetometers if (hgm(n,j) < 9998.0 .and. abs(mglat(j)) > 55.0) then MagMean(1) = MagMean(1) + abs(hgm(n,j)) MagNpts(1) = MagNpts(1) + 1 endif if (dgm(n,j) < 9998.0 .and. abs(mglat(j)) > 55.0) then MagMean(2) = MagMean(2) + abs(dgm(n,j)) MagNpts(2) = MagNpts(2) + 1 endif if (zgm(n,j) < 9998.0 .and. abs(mglat(j)) > 55.0) then MagMean(3) = MagMean(3) + abs(zgm(n,j)) MagNpts(3) = MagNpts(3) + 1 endif enddo if (MagNpts(1) > 0) then MagMean = MagMean/MagNpts do j = 1, nmagnetometers if (hgm(n,j) < 9998.0 .and. abs(mglat(j)) > 55.0) then MagStd(1) = MagStd(1) + abs(MagMean(1)-abs(hgm(n,j)))/MagNpts(1) endif if (hgm(n,j) < 9998.0 .and. abs(mglat(j)) > 55.0) then MagStd(2) = MagStd(2) + abs(MagMean(2)-abs(hgm(n,j)))/MagNpts(2) endif if (zgm(n,j) < 9998.0 .and. abs(mglat(j)) > 55.0) then MagStd(3) = MagStd(3) + abs(MagMean(3)-abs(zgm(n,j)))/MagNpts(3) endif enddo do j = 1, nmagnetometers IsBad = .false. if (abs(hgm(n,j)) > MagMean(1) + 20.0*MagStd(1) .and. & hgm(n,j)<1.0e4) then IsBad = .true. iMagErrorFlag(n,j,1) = 10 endif if (abs(dgm(n,j)) > MagMean(2) + 20.0*MagStd(2) .and. & dgm(n,j)<1.0e4) then IsBad = .true. iMagErrorFlag(n,j,2) = 10 endif if (abs(zgm(n,j)) > MagMean(3) + 20.0*MagStd(3) .and. & zgm(n,j)<1.0e4) then IsBad = .true. iMagErrorFlag(n,j,3) = 10 endif if (IsBad) then write(*,*) "Found Bad Magnetometer Station: ", stname(j) write(*,*) " at time :",itime write(*,*) " H, mean(H), stddev(H) : ", & hgm(n,j), MagMean(1), MagStd(1) write(*,*) " E, mean(E), stddev(E) : ", & dgm(n,j), MagMean(2), MagStd(2) write(*,*) " Z, mean(Z), stddev(Z) : ", & zgm(n,j), MagMean(3), MagStd(3) endif enddo endif call calc_ae(n) !\ ! Calculate weights on different magnetometers !/ do j = 1, nmagnetometers x = 0.5 * (abs(mglat(j)) - 30.0) if (x > 15.0) x = 15.0 if (x < 15.0) x = 5.0 x = x*x*(1.0 + ae(n)*ae(n)/40000.0) if (hgm(n,j) < 9998.0 .and. iMagErrorFlag(n,j,1) < 10 .and. & iMagErrorFlag(n,j,2) < 10) then wth(n,j) = weight_on_x / (erhgm*erhgm + x) wtd(n,j) = weight_on_y / (erdgm*erdgm + x) else wth(n,j) = 0.0 wtd(n,j) = 0.0 endif ! wth(n,j) = weight_on_x ! wtd(n,j) = weight_on_y if (use_zed) then if (abs(mglat(j)) >= 70.0) & x = x + (1.5*(abs(mglat(j)) - 70.0))**2 if (hgm(n,j) < 9998.0 .and. iMagErrorFlag(n,j,2) < 10) then wtzed(n,j) = weight_on_z / (erzgm*erzgm + x) else wtzed(n,j) = 0.0 endif ! wtzed(n,j) = weight_on_z else zgm(n,j) = 99999.0 wtzed(n,j) = 0.0 endif enddo call calc_dst(n,5,1) if (n < ntimes) then if (timelist(n+1) > -1.0e32) then ctime = timelist(n+1) else ctime = ctime + dt endif nlines_to_read = (ctime-time-dt_magfile)/dt_magfile*nmagnetometers if (nlines_to_read > 1) then do i = 1, nlines_to_read read (iunit,'(a)') dummy enddo endif endif enddo close(iunit) end subroutine read_magnetometers subroutine ahn use ModAMIE use ModData use ModConductance use ModConstants implicit none integer :: n, j real :: latwt, Hal, Ped, MaxDh, LatMaxDh n = iteration_number !\ ! Calculate Ahn conductance from magnetometer observations !/ if (use_ahn) then NAhn = 0 MaxDh = 0.0 LatMaxDh = -999.0 do j=1,nmagnetometers if (mglat(j)*ihsoln >= 50.0 .and. & mglat(j)*ihsoln <= 78.0 .and. & wth(n,j) > 0.0 .and. & abs(hgm(n,j)) > MaxDh) then MaxDh = abs(hgm(n,j)) LatMaxDh = mglat(j)*ihsoln endif enddo do j=1,nmagnetometers if (mglat(j)*ihsoln >= 50.0 .and. & mglat(j)*ihsoln <= 78.0 .and. & wth(n,j) > 0.0 .and. & wtd(n,j) > 0.0) then if (index(ahn_model,'83') > 0) & call ahn83(hgm(n,j),Hal,Ped) if (index(ahn_model,'91') > 0) & call ahn91(hgm(n,j),zgm(n,j),Hal,Ped) if (index(ahn_model,'98') > 0) & call ahn98(amltmg(j),hgm(n,j),zgm(n,j),Hal,Ped) if (hal > 0.0 .and. ped > 0.0) then NAhn = NAhn + 1 AhnLat(NAhn) = mglat(j) AhnMLT(NAhn) = amltmg(j) AhnPed(NAhn) = ped AhnHal(NAhn) = hal if (AhnPed(NAhn) < ConductanceBackground(pedersen_)) & AhnPed(NAhn) = ConductanceBackground(pedersen_) if (AhnHal(NAhn) < ConductanceBackground(hall_)) & AhnHal(NAhn) = ConductanceBackground(hall_) ! if (UseExperimentalCode) then ! latwt = abs(mglat(j)-LatMaxDh)/10.0 ! AhnWeight(NAhn) = 4.0 / (1.+latwt*latwt) ! else latwt = ((90.-abs(mglat(j)))*dtor-0.426)/0.12 AhnWeight(NAhn) = 4.0 / (1.+0.25*latwt*latwt) ! endif ! Reduce weights for stations between 50 and 60 MLAT latwt = 0.5 + (1.0 - exp(-(abs(mglat(j))-50.0)**2.0/50.0))*0.5 AhnWeight(NAhn) = AhnWeight(NAhn) * latwt ! if (abs(mglat(j)) > 50. .and. abs(mglat(j)) < 60.) & ! AhnWeight(NAhn) = AhnWeight(NAhn) * 0.5 AhnWeight_EFlux(NAhn) = AhnWeight(NAhn) * 0.3 AhnWeight_AveE(NAhn) = AhnWeight(NAhn) * 2.0 if (AhnWeight(NAhn) > 10.0) AhnWeight(NAhn) = 10.0 if (AhnWeight_EFlux(NAhn) > 20.0) AhnWeight_EFlux(NAhn) = 20.0 if (AhnWeight_AveE(NAhn) > 20.0) AhnWeight_AveE(NAhn) = 20.0 endif endif enddo else AhnHal = -1.0 AhnPed = -1.0 AhnWeight = -1.0 NAhn = 0 endif end subroutine ahn subroutine ahn91(h,z,Hal,Ped) real, intent(in) :: h,z real, intent(out) :: Hal, Ped if (h < 0.0) then if (z > 9998.0) then Ped = 2.432*(-h)**.235 Hal = 4.355*(-h)**.340 Hal = 0.5 * Hal else if (h > -30.0) then if (z >= 0.0) then Hal = 1.22 * (-h)**0.43 else Hal = 1.98 * (-h)**0.46 endif else if (z >= 0.) then Hal = 1.53 * (-h)**0.42 else Hal = 1.73 * (-h)**0.51 endif endif if (z >= 0.0) then Ped = exp(alog(Hal/0.86)/1.23) else Ped = exp(alog(Hal/1.63)/1.04) endif endif endif if (h > 0.0) then if (z > 9998.0) then Ped = 0.699*(h)**0.432 Hal = 2.096*(h)**0.357 Hal = 0.5 * Hal else if (h < 30.0) then if (z >= 0.0) then Hal = 1.23 * (h)**0.34 else Hal = 0.85 * (h)**0.48 endif else if (z >= 0.) then Hal = 0.40 * (h)**0.63 else Hal = 0.19 * (h)**0.85 endif endif Ped = exp(alog(Hal/0.60)/1.33) endif endif if (h == 0.0) then Ped = -1.0 Hal = -1.0 endif end subroutine ahn91 !----------------------------------------------------------------------- subroutine ahn98 (rmlt,h,z,Hal,Ped) ! Use Ahn et al (1998, JGR, 103, 14,769-14,780) model ! If z=missing, or if h=0 (h=missing checked with wth=0 before this ! call), or if mlt region not there, skip this point. ! If all z missing, use ahn83, and change call in putcsrc. ! If abs(z)<15 nT(=zs), average solutions for z + and -. ! Have no model for 9.5>rmlt>12.5, so use ahn83 model then use ModAMIE, only : UseExperimentalCode implicit none real, intent(in) :: rmlt, h, z real, intent(out) :: Hal, Ped real, dimension(29) :: ah, bh, ap, bp integer :: is real, parameter :: zs = 15.0 ! Hal = ah*abs(h)**bh, Ped = ap*abs(h)**bp, where conductances ! are in Siemens (mhos) and delta(H) and delta(Z) are in nT. ! 31 categories: ! 1-4 for h>0,z<-zs (16-18, 19-21, 22-23, 00-02 MLT) ! 5-9 for h>0,z>zs (13-15, 16-18, 19-21, 22-23, 00-02 MLT) ! 10-14 for h<0,z>zs (19-21, 22-23, 00-02, 03-04, 05-06 MLT) ! 15-20 for h<0,z<-zs (19-21, 22-23, 00-02, 03-04, 05-06, 07-09 MLT) ! For -zs0,z small (16-18, 19-21, 22-23, 00-02 MLT) (5 for 13-15 MLT) ! 25-29 for h<0,z small(19-21,22-23,00-02,03-04,05-06 MLT)(20 for 07-09 MLT) data ah/ 1.190,5.393,8.519,5.247, 3.385,0.921,0.527,0.940,0.830, & 0.465,0.211,1.570,0.111,4.407, 4.820,2.536,2.966,2.490,1.501, & 4.541, 1.0555,2.9600,4.7295,3.0385, 2.6425,1.3735,2.2680, & 1.3005,2.9540/ data bh/ 0.373,0.174,0.125,0.136, 0.129,0.424,0.600,0.618,0.665, & 0.693,0.793,0.419,0.870,0.260, 0.316,0.465,0.426,0.429,0.509, & 0.286, 0.3985,0.3870,0.3715,0.4005, 0.5045,0.6290,0.4225, & 0.6495,0.3845/ data ap/ 2.446,3.813,6.509,2.146, 2.382,0.310,1.193,1.073,1.796, & 1.953,0.727,2.876,0.504,5.119, 12.015,4.699,5.017,3.048,1.244, & 1.574, 1.3780,2.5030,3.7910,1.971, 6.9840,2.7130,3.9465, & 1.7760,3.1815/ data bp/ 0.260,0.223,0.113,0.358, 0.201,0.657,0.392,0.524,0.376, & 0.338,0.466,0.223,0.517,0.116, 0.013,0.197,0.185,0.240,0.364, & 0.284, 0.4585,0.3075,0.3185,0.3670, 0.1755,0.3315,0.2040, & 0.3785,0.2400/ ! Check for 9>MLT>13 ! if (UseExperimentalCode) then ! if (rmlt > 9.5 .and. rmlt < 16.0) then ! call ahn91(h,z,Hal,Ped) ! return ! endif ! else if (rmlt > 9.5 .and. rmlt < 12.5) then call ahn83(h,Hal,Ped) return endif ! endif is = 0 ! Check for h=0 or for missing z if (abs(h) > 0. .and. abs(z) < 9998.) then ! h positive (Eastward electrojet): if (h > 0.0) then ! Check for MLT out of bounds if ((rmlt > 2.5 .and. rmlt < 12.5) .or. & (rmlt.gt.2.5 .and. rmlt.lt.15.5 .and. z.lt.-zs)) then is = 0 else ! 13-15 MLT, only have equatorward (z>0) eastward electrojet if (rmlt .ge. 12.5 .and. rmlt .lt. 15.5) is = 5 ! Small z if (abs(z) .le. zs) then if (rmlt .ge. 15.5 .and. rmlt .lt. 18.5) is=21 if (rmlt .ge. 18.5 .and. rmlt .lt. 21.5) is=22 if (rmlt .ge. 21.5 .and. rmlt .lt. 23.5) is=23 if (rmlt .ge. 23.5 .or. rmlt .le. 2.5) is=24 else ! Larger z if (z .lt. -zs) then ! Poleward half of eastward electrojet (z negative and < -zs) if (rmlt .ge. 15.5 .and. rmlt .lt. 18.5) is=1 if (rmlt .ge. 18.5 .and. rmlt .lt. 21.5) is=2 if (rmlt .ge. 21.5 .and. rmlt .lt. 23.5) is=3 if (rmlt .ge. 23.5 .or. rmlt .le. 2.5) is=4 else ! Equatorward half of eastward electrojet ! (z positive and > zs) if (rmlt .ge. 15.5 .and. rmlt .lt. 18.5) is=6 if (rmlt .ge. 18.5 .and. rmlt .lt. 21.5) is=7 if (rmlt .ge. 21.5 .and. rmlt .lt. 23.5) is=8 if (rmlt .ge. 23.5 .or. rmlt .le. 2.5) is=9 endif endif endif else ! h negative (Westward electrojet): ! Check for MLT out of bounds if ((rmlt .gt. 9.5 .and. rmlt .lt. 18.5) .or. & (rmlt.gt.6.5 .and. rmlt.lt.18.5 .and. z.gt.zs)) then is = 0 else ! 07-09 MLT, only have equatorward (z<0) westward electrojet if (rmlt .ge. 6.5 .and. rmlt .lt. 15.5) is = 20 ! Small z if (abs(z) .le. zs) then if (rmlt .ge. 18.5 .and. rmlt .lt. 21.5) is=25 if (rmlt .ge. 21.5 .and. rmlt .lt. 23.5) is=26 if (rmlt .ge. 23.5 .or. rmlt .lt. 2.5) is=27 if (rmlt .ge. 2.5 .and. rmlt .lt. 4.5) is=28 if (rmlt .ge. 4.5 .and. rmlt .lt. 6.5) is=29 else ! Larger z ! Poleward half of westward electrojet (z positive and > zs) if (z .lt. -zs) then if (rmlt .ge. 18.5 .and. rmlt .lt. 21.5) is=10 if (rmlt .ge. 21.5 .and. rmlt .lt. 23.5) is=11 if (rmlt .ge. 23.5 .or. rmlt .lt. 2.5) is=12 if (rmlt .ge. 2.5 .and. rmlt .lt. 4.5) is=13 if (rmlt .ge. 4.5 .and. rmlt .lt. 6.5) is=14 else ! Equatorward half of westward electrojet (z negative and < -zs) if (rmlt .ge. 18.5 .and. rmlt .lt. 21.5) is=15 if (rmlt .ge. 21.5 .and. rmlt .lt. 23.5) is=16 if (rmlt .ge. 23.5 .or. rmlt .lt. 2.5) is=17 if (rmlt .ge. 2.5 .and. rmlt .lt. 4.5) is=18 if (rmlt .ge. 4.5 .and. rmlt .lt. 6.5) is=19 endif endif endif endif endif if (is > 0) then Hal = ah(is)*abs(h)**bh(is) Ped = ap(is)*abs(h)**bp(is) else Hal = -1.0 Ped = -1.0 endif return end subroutine ahn98 !----------------------------------------------------------------------- subroutine ahn83 (h,Hal,Ped) real, intent(in) :: h real, intent(out) :: Hal, Ped ! Use Ahn's old formulas based also on delta H alone. ! Ahn et al, PSS, 641-653, 1996. if (h < 0.) then Ped = 2.432*(-H)**.235 Hal = 4.355*(-H)**.340 else ! Positive del(H) (and zero!) Ped = 0.699*(H)**.432 Hal = 2.096*(H)**.357 endif return end subroutine ahn83 !----------------------------------------------------------------------- subroutine calc_ae(n) use ModIndices use ModData use ModConstants implicit none integer, intent(in) :: n integer :: j, iMin, iMax real :: dx, mini, maxi !\ ! AE Calculation !/ if (.not. UseMinMaxAE) then do j = 1, nmagnetometers if (abs(mglat(j)) < 76.0 .and. & abs(mglat(j)) > 55.0 .and. & hgm(n,j) < 9998.0 .and. & iMagErrorFlag(n,j,1) < 2) then if (hgm(n,j) > maxi) then maxi = hgm(n,j) iMax = j endif if (hgm(n,j) < mini) then mini = hgm(n,j) iMin = j endif endif enddo else mini = -1.0e32 maxi = 1.0e32 iMin = 1 iMax = 1 endif do j = 1, nmagnetometers dx = 0.0 if (abs(mglat(j)) < 76.0 .and. & abs(mglat(j)) > 55.0 .and. & hgm(n,j) < 9998.0 .and. & iMagErrorFlag(n,j,1) < 2) then nae(n) = nae(n) + 1 dx = hgm(n,j) ! Here we take the SECOND largest and SECOND smallest values. ! This may eliminate spikes. ! (only really true if UseMinMaxAE is .false.) if (dx > au(n) .and. dx < maxi) au(n) = dx if (dx < al(n) .and. dx > mini) al(n) = dx endif enddo ae(n) = au(n) - al(n) Indices(n,au_) = au(n) Indices(n,al_) = al(n) Indices(n,ae_) = ae(n) if (maxi > au(n) + 50.0 .and. .not.UseMinMaxAE) then iMagErrorFlag(n,iMax,:) = 11 if (DebugLevel > 0) & write(*,*) "=> Bad Station in AE : ",n,mgname(iMax), maxi, au(n) endif if (mini < al(n) - 100.0 .and. .not.UseMinMaxAE) then iMagErrorFlag(n,iMin,:) = 11 if (DebugLevel > 0) & write(*,*) "=> Bad Station in AE : ",n,mgname(iMin), mini, al(n) endif end subroutine calc_ae !----------------------------------------------------------------------- subroutine calc_dst(n,iMaxErrorFlag, iMethod) use ModIndices use ModData use ModConstants implicit none integer, intent(in) :: n,iMaxErrorFlag, iMethod integer :: j, k, l, ll, iMini, iMaxi, nStats real, dimension(MAX_MAGNETOMETERS) :: dstval real, dimension(0:MAX_MAGNETOMETERS) :: dstlon real :: dstlsv, dstvsv, clat, maxi, mini, maxi2, mini2 !\ ! Calculate Dst !/ maxi = -1000.0 mini = 1000.0 maxi2 = -1000.0 mini2 = 1000.0 nDst(n) = 0 nStats = 0 do j = 1, nmagnetometers if (abs(mglat(j)) < 40.0 .and. wth(n,j) /= 0.0 .and. & hgm(n,j) < 9998.0 .and. & iMagErrorFlag(n,j,1) <= iMaxErrorFlag) then nStats = nStats + 1 endif if (abs(mglat(j)) < 40.0 .and. wth(n,j) /= 0.0 .and. & hgm(n,j) < 9998.0 .and. & iMagErrorFlag(n,j,1) <= iMaxErrorFlag .and. & hgm(n,j) > maxi) then maxi = hgm(n,j) iMaxi = j endif if (abs(mglat(j)) < 40.0 .and. wth(n,j) /= 0.0 .and. & hgm(n,j) < 9998.0 .and. & iMagErrorFlag(n,j,1) <= iMaxErrorFlag .and. & hgm(n,j) < mini) then mini = hgm(n,j) iMini = j endif enddo !write(*,*) "nStats : ", nStats, iMaxErrorFlag, mini, maxi ! If we don't have enough stations, then we need to loosen our ! restrictive criteria if (nStats <= 2 .or. UseMinMaxDst) then mini = -1.0e32 maxi = 1.0e32 endif do j = 1, nmagnetometers ! Notice that the < and > don't include =. if (abs(mglat(j)) < 40.0 .and. wth(n,j) /= 0.0 .and. & hgm(n,j) < 9998.0 .and. & iMagErrorFlag(n,j,1) <= iMaxErrorFlag .and. & hgm(n,j) > mini .and. hgm(n,j) < maxi) then if (hgm(n,j) < mini2) mini2 = hgm(n,j) if (hgm(n,j) > maxi2) maxi2 = hgm(n,j) clat = (90. - mglat(j)) * dtor ndst(n) = ndst(n) + 1 dstlon(ndst(n)) = mglon(j) dstval(ndst(n)) = hgm(n,j)/sin(clat) if (ndst(n).gt.1) then k = 1 do while (k <= ndst(n)-1) if (dstlon(ndst(n)) < dstlon(k)) then dstlsv = dstlon(ndst(n)) dstvsv = dstval(ndst(n)) do ll=k,ndst(n)-1 l = ndst(n) + k - ll dstlon(l) = dstlon(l-1) dstval(l) = dstval(l-1) enddo dstlon(k) = dstlsv dstval(k) = dstvsv k = ndst(n) else k = k + 1 endif enddo endif endif enddo Indices(n,dst_) = 0.0 if (ndst(n) > 0) then dstlon(ndst(n)+1) = dstlon(1) + 360. dstlon(0) = dstlon(ndst(n)) - 360. do j=1,ndst(n) Indices(n,dst_) = Indices(n,dst_) + & dstval(j)*(dstlon(j+1) - dstlon(j-1)) enddo if (iMethod == 2) & Indices(n,dst_) = Indices(n,dst_)/720. if (iMethod == 1 .or. iMethod == 3) & Indices(n,dst_) = sum(dstval(1:ndst(n)))/ndst(n) endif if (nStats > 2 .and. .not.UseMinMaxDst) then if (maxi > maxi2 + 30.0) then iMagErrorFlag(n,iMaxi,:) = 11 ! if (DebugLevel > 0) & write(*,*) "=> Bad Station in Dst : ",mgname(iMaxi), maxi, maxi2 endif if (mini < mini2 - 30.0) then iMagErrorFlag(n,iMini,:) = 11 ! if (DebugLevel > 0) & write(*,*) "=> Bad Station in Dst : ",mgname(iMini), mini, mini2 endif endif !write(*,*) "Dst : ",n,ndst(n), Indices(n,dst_) end subroutine calc_dst