program mercon C 11/93: ~emery/geomag/mercon.f C This is a combination of the old merge.f and convert.f codes, with the C additional ability to add or replace .rem(m,c)a files to a gmag file. C Merge found average standard deviations of each component in bins of C 10 degrees of latitude. These were supposed to be used as the standard C deviations for hourly extrapolated files (.rem(m,c)x), but a standard C number of 20 nT was used instead. This will not be done here, but will C be dealt with in the codes for the hourly files. Convert converted from C geographic (XYZ) or local magnetic (DHZ) coordinates to magnetic (usually C apex) coordinates. It also removed any of the quiet day Dst from Xm at C the very end (10/91). C 11/94: Add another file output of the calculation of AE, AU, AL, C AL(21-01MLT), Dst for plotting later in an IDL program. C NOTE: Should smooth out removal of quiet day Dst, and add addition C of active day Dst here instead of in AMIE. C ntims >= number of times per day * number of days. C parameter (ntims = 1440*1+1) parameter (ntims = 1728*1) parameter (nstatp = 280) REAL GGLAT(nstatp),GGLONG(nstatp),CGLAT(nstatp), | CGLONG(nstatp),DEC(nstatp),psi(nstatp) dimension cosp(nstatp),sinp(nstatp),cosd(nstatp),sind(nstatp) real dstq(24,9),qdst(nstatp),qdstc(24,nstatp),dstqav(9) integer jqd(nstatp),iskp5(nstatp),iadrv(nstatp), | num(3,nstatp,ntims) dimension dstlon(nstatp),dstval(nstatp),dst1h(24,15), | ifhxmn(nstatp) dimension dsthrly(360),tmhrly(360) integer mn(ntims),ih(ntims),iy(ntims),im(ntims),id(ntims) C Must be dimensioned by nsta and ntims, the number of stations and times integer d(nstatp,ntims),h(nstatp,ntims),z(nstatp,ntims) integer dstd(nstatp,ntims),hstd(nstatp,ntims), | zstd(nstatp,ntims) integer iqdst(nstatp,ntims),orig(3,nstatp,ntims), | mean(3,nstatp,ntims) integer ixmn(nstatp,ntims),iymn(nstatp,ntims), | izmn(nstatp,ntims) integer oldxd(nstatp,ntims),oldyh(nstatp,ntims), | mean3(3,nstatp),num3(3,nstatp) character*80 line character*3 stat(nstatp),comp1(nstatp),comp2(nstatp), | dhz3,ihxid(nstatp) character*3 statp(nstatp),stan(12) character*40 fname(nstatp),filein character*80 latfile, magfile, aefile ! INTEGER D,H,Z,dstd,hstd,zstd,spval data dhz3/'DHZ'/ data spval /99999/ DATA D2R/ 0.017453293 / C iskp5=1 means skip over 5 lines of info on daily averages (true of C newer .remma or .remca files) data iskp5 /nstatp*1/ C Order in apex latitude C **************************** may change *********************************** C iadrv = -1,0,1,2 if ignore, use old, add or revise this station in gmag. C If all equal 1, then there is no gmag to replace. A value of -1 means C the station is considered for a future gmag, but not on any here. The C total number of 1's and of -1's in the first nsta values is the additional C number of stations in the psi file read and in the fname list compared C to the older gmag read. c data iadrv/nstatp*1/ C data iadrv/ 46*1,104*-1/ data iadrv/nstatp*1/ c data (stan(i),i=1,12)/'TRO','DIX','CCS','TIK','WEL','BRW', c | 'CMO','YKC','FCC','PDB','NAQ','LRV'/ C for May 96 data (stan(i),i=1,12)/'TRO','DIK','CCS','TIX','CWE','BRW', | 'CMO','YKC','FCC','PDB','NAQ','LRV'/ C ifhxmn = 1,0 if the X or H component uses a mean (ic1,2npav=288 in C remcurve. If so, remove the mean of the quiet Dst, not the curve. data ifhxmn/ nstatp*0/ C May give a list of nsta input file names as fname(i),i=1,nsta C Can use up to 9 quiet days. nqd is the number of days used, and jqd C gives the number (1-nqd) of the particular day used for each station. data nqd/1/ data jqd /nstatp*1/ nsta = 103 ndays = 2 nda1 = 2 ntmpd = 288 write(6,*) 'Enter number of stations :' read(5,*) nsta write(6,*) 'Enter number of days :' read(5,*) ndays write(6,*) 'Enter beginning day :' read(5,*) nda1 write(6,*) 'Enter number of time periods per day :' read(5,*) ntmpd ntime = ndays*ntmpd do i=1,nsta write(6,*) 'Enter file for station number ',i,' :' read(5,*) fname(i) enddo do i=1,24 write(6,*) 'Enter dst for quiet day, hour ',i,' :' read(5,*) dstq(i,1) enddo do j=1,ndays do i=1,24 write(6,*) 'Enter dst for day',j+nda1-1,' hour ',i,' :' read(5,*) dst1h(i,j) enddo enddo write(6,*) 'Enter file which contains mlats for stations :' read(5,'(A80)') latfile write(6,*) 'Enter file to output magnetometer data :' read(5,'(A80)') magfile write(6,*) 'Enter file to output dst and ae :' read(5,'(A80)') aefile C nda1 = day of month of first of the ndays (assume for now in 1 month) nhr0 = 0 nhrpd = 24 open(1,file=latfile,status='old' ) open(22,file=magfile, status='new') open (31,file=aefile,status='new') isumign = 0 isumold = 0 isumadd = 0 isumrev = 0 do 300 i=1,nsta go to (291,292,293,294), iadrv(i)+2 291 isumign = isumign + 1 go to 300 292 isumold = isumold + 1 go to 300 293 isumadd = isumadd + 1 go to 300 294 isumrev = isumrev + 1 300 continue C write(6,"('isumign,isumold,isumadd,isumrev',4i4)") C | isumign,isumold,isumadd,isumrev C write (6,"(1x,'# stations on old gmag = # old + # revised =', C | i3,'+',i3,'=',i3/1x,' + # add =',i3,'=',i3,'= # on new gmag', C | 1x,'# stations ignored (on psi, but not on old or new gmag) =', C | i3,' Sum of all =',i3,'= nsta =',i3)") isumold,isumrev, C | isumold+isumrev,isumadd,isumold+isumrev+isumadd,isumign, C | isumold+isumrev+isumadd+isumign,nsta C if (isumold+isumrev .ne. 0) C | open (21,file='gmag.jan92C33b',status='old') C **************************** end change *********************************** RAD = 180./3.1415926535898 ntmph = ntmpd / nhrpd ntm = ndays * ntmpd if (ntm .gt. ntims) then write (6,"(1x,'mercon: increase ntims to >= ntm: ntm ntims=', | 2i6)") ntm,ntims stop endif C Put ndays Dst into dsthrly and compute the time j=0 nhrs = ndays*24 do 40 n=1,ndays do 40 i=1,24 j = j + 1 dsthrly(j) = dst1h(i,n) 40 tmhrly(j) = float(n+nda1-1) + (float(i)-0.5)/24. C Compute average quiet day Dst for stations which only remove a mean do 55 k=1,nqd dstqav(k) = 0. do 50 i=1,24 50 dstqav(k) = dstqav(k) + dstq(i,k) 55 dstqav(k) = dstqav(k) / 24. write (6,"(1x,'Average quiet day Dst =',9f7.1)") (dstqav(k), | k=1,nqd) C Read number of stations and number of times from psi file read (1,"(i3,6x,i5)") ns,ntime2 read (1,"(a80)") line nae = 0 ndst = 0 DO 100 I = 1,NS READ (1,"(i3,1x,a3,6x,a3,4x,2f7.2,1x,2f7.2,f9.2,f7.2)") | ndx,statp(i),ihxid(i),gglat(i),gglong(i),cglat(i), + cglong(i),dec(i),psi(i) if (abs(cglat(i)).ge.55. .and. abs(cglat(i)).le.76.) nae=nae+1 if (abs(cglat(i)) .le. 40.) ndst = ndst + 1 if (abs(cglat(i)).gt.55) ifhxmn(i) = 1 filein = fname(i) cosp(i) = cos(psi(i)*d2r) sinp(i) = sin(psi(i)*d2r) cosd(i) = cos(dec(i)*d2r) sind(i) = sin(dec(i)*d2r) C if (ihxid(i) .eq. 'HDZ') ihxid(i) = dhz3 cosa = cos(abs(cglat(i))*d2r) nq = jqd(i) qdst(i) = dstqav(nq)*cosa do 944 k=1,24 944 qdstc(k,i) = dstq(k,nq)*cosa C write (6,"(1x,4i3,2f6.1,1x,a4,3(2f7.2,1x) )") i,ndx,ifhxmn(i), C | iadrv(i),qdst(i),qdstc(1,i),statp(i),gglat(i),gglong(i), C | cglat(i),cglong(i),dec(i),psi(i) 100 CONTINUE close(1) write (6,"(1x,'read ns =',i3,' stations in PSI file')") ns if (ns .ne. nsta) then write (6,"(1x,'ns .ne. nsta, ns nsta =',2i4,' stop')") | ns,nsta stop endif C Write out header for AE, Dst file write (31,"('AE is calculated from',i3,' stations between 55-76' | ' mlat'/'Dst is calculated from',i3,' stations below 40 mlat'/ | 'Provisional Dst comes from WDC-C2 from SJG, HON, KAK, HER')") | nae,ndst c write (31,"(' yr mo da hr mn nae AE(nT) AU(nT) AL(nT)', c | ' # AL(21-1) ndst Dst(nT) prov Dst')") c write (31,"(' yr mo da hr mn nae AE(nT) AU(nT) AL(nT)', c | ' # AL(21-1) ndst Dst(nT) prov Dst Stand AE')") write (31,"(' yr mo da hr mn nae AE(nT) AU(nT) AL(nT)', | ' # AL(21-1) ndst Dst(nT) prov Dst Stand AE', | ' ASDup ASDdw')") C Read .rem,ca files first do 400 i = 1,nsta if (iadrv(i) .le. 0) go to 400 open(1,file=fname(i),status='old' ) filein = fname(i) C Read over 5 lines at beginning on daily quiet aves for the latest files if (iskp5(i) .eq. 1) then read (1,"(a75)") char75 read (1,"(a75)") char75 do 58 k=1,3 read (1,"(3x,i4,i6)") num3(k,i),mean3(k,i) do 58 j=1,ntime num(k,i,j) = num3(k,i) 58 mean(k,i,j) = mean3(k,i) endif n = 1 icloc = 1 do while (filein(icloc:icloc).ne.'.') icloc = icloc + 1 enddo icloc = icloc + 4 C do 60 n=1,ndays do while (n.le.ndays) do 60 kh=nhr0,nhr0+nhrpd-1 do 59 m = 1,ntmph j = (n-1)*ntmpd + kh*ntmph + m if (filein(icloc:icloc) .eq. 'c') then read(1,"(5i2,a3,1x,a3,15i6)") iy(j),im(j), | id(j),ih(j),mn(j),stat(i),comp1(i), | h(i,j),d(i,j),z(i,j),dstd(i,j),hstd(i,j), | zstd(i,j),(orig(k1,i,j),k1=1,3), | (mean(k2,i,j),k2=1,3),(num(k3,i,j),k3=1,3) C***temp*** change signs for TIX,CHD and KTN on Oct/95 C if (iy(j) .eq. 95 .and. im(j) .eq. 10) then C if (stat(i) .eq. 'TIX' .or. C | stat(i) .eq. 'CHD' .or. C | stat(i) .eq. 'KTN') then C print *,'change sign for ',stat(i) C if (h(i,j) .ne. 99999) h(i,j) = -h(i,j) C if (d(i,j) .ne. 99999) d(i,j) = -d(i,j) C if (z(i,j) .ne. 99999) z(i,j) = -z(i,j) C endif C endif else read(1,"(5i2,a3,2x,a3,9i6)") iy(j),im(j), | id(j),ih(j),mn(j), stat(i),comp1(i),h(i,j), | d(i,j),z(i,j),dstd(i,j),hstd(i,j), | zstd(i,j),(orig(k1,i,j),k1=1,3) C***temp*** change signs for TIX,CHD and KTN on Oct/95 C if (iy(j) .eq. 95 .and. im(j) .eq. 10) then C if (stat(i) .eq. 'TIX' .or. C | stat(i) .eq. 'CHD' .or. C | stat(i) .eq. 'KTN') then C if (h(i,j) .ne. 99999) h(i,j) = -h(i,j) C if (d(i,j) .ne. 99999) d(i,j) = -d(i,j) C if (z(i,j) .ne. 99999) z(i,j) = -z(i,j) C endif C endif endif 59 continue 60 continue if (id(j).ge.nda1 .and. id(j).lt.nda1+ndays) n = n+1 enddo close(1) 400 continue C Now read old gmag file if appropriate, and create new one C xyz data is held as: x y z C dhz data is held as: d h z write(6,*) "ntime : ",ntime DO 900 j=1,ntime C Get interpolated Dst from hourly provisional data tmhr = float(id(j)) + float(ih(j))/24. + float(mn(j))/1440. call iterp (tmhrly,dsthrly,tmhr,nhrs,dstint,nn,dx1,dx2) C Calculate Ae from available magnetometer stations using d(X), if C X is in the horizontal magnetic north direction. (Is very close C to d(H), which is the official definition of Ae). C Only use stations between 55 and 75 degrees magnetic latitude. nae = 0 AUCALC = -9999. ALCALC = 9999. austan = -9999. alstan = 9999. asdup = -9999. asddw = 9999. istan = 0 C Also find Al for 21-01 MLT region to see evidence of substorms ALMIDN = 9999. nalm = 0 C ROT = magnetic longitude of midnight ROT = 70. - 15.*ih(j) - mn(j)/4. C CALCULATE Dst NDST = 0 dstcalc = 0. DO 600 i=1,nsta DSTLON(i) = 0. 600 DSTVAL(i) = 0. DO 200 I = 1,nsta ndx = i if (iadrv(i) .eq. -1) go to 200 if (iadrv(i) .eq. 0) then go to 190 endif if (stat(i) .ne. statp(i)) then write (6,"(1x,'i stat statp =',i4,1x,a3,1x,a3,' warning!!!')") | i,stat(i),statp(i) c stop endif comp2(i) = 'XYZ' ipr = 0 if (j .eq. 1 .or. j .eq. 289 .or. j .eq. 577 .or. j .eq. 864) | ipr = 1 if (ipr .eq. 1) write (6,"(1x,'orig:',3i2,i3,i2,1x,a3,1x,a3, | 6i6)") iy(j),im(j),id(j),ih(j),mn(j),stat(i),comp1(i),h(i,j), | d(i,j),z(i,j),dstd(i,j),hstd(i,j),zstd(i,j) C Find deviations from mean (0 if use mean, non-zero if use curve) IF ( COMP1(i).EQ.'XYZ' ) THEN xgmn = mean(1,i,j) - mean3(1,i) ygmn = mean(2,i,j) - mean3(2,i) ELSE hmn = mean(1,i,j) - mean3(1,i) dmn = mean(2,i,j) - mean3(2,i) xgmn = hmn*cosd(i) - dmn*sind(i) ygmn = hmn*sind(i) + dmn*cosd(i) ENDIF izmn(i,j) = mean(3,i,j) - mean3(3,i) xmn = xgmn*cosp(i) + ygmn*sinp(i) ymn = ygmn*cosp(i) - xgmn*sinp(i) ixmn(i,j) = nint(xmn) iymn(i,j) = nint(ymn) IF ( (abs(D(i,j)).LT.99999).AND.(abs(H(i,j)).LT.99999) ) THEN IF ( COMP1(i).EQ.'XYZ' ) THEN xg = h(i,j) yg = d(i,j) xgstd = hstd(i,j) ygstd = dstd(i,j) ELSE xg = h(i,j)*cosd(i) - d(i,j)*sind(i) yg = h(i,j)*sind(i) + d(i,j)*cosd(i) xgstd = sqrt( (hstd(i,j)*cosd(i))**2 + (dstd(i,j)* | sind(i))**2 ) ygstd = sqrt( (hstd(i,j)*sind(i))**2 + (dstd(i,j)* | cosd(i))**2 ) ENDIF xm = xg*cosp(i) + yg*sinp(i) ym = yg*cosp(i) - xg*sinp(i) oldxd(i,j) = h(i,j) oldyh(i,j) = d(i,j) h(i,j) =nint(xm) d(i,j) =nint(ym) xmsd = sqrt( (xgstd*cosp(i))**2 + (ygstd*sinp(i))**2 ) ymsd = sqrt( (xgstd*sinp(i))**2 + (ygstd*cosp(i))**2 ) hstd(i,j) =nint(xmsd) dstd(i,j) =nint(ymsd) ELSE oldxd(i,j) = spval oldyh(i,j) = spval h(i,j) = spval d(i,j) = spval hstd(i,j) = spval dstd(i,j) = spval ENDIF C Add back any quiet day Dst rqdst = qdst(i) C Do different Dst correction if remove a curve if (ifhxmn(i) .eq. 0) then C Interpolate quiet hourly Dst to avoid discontinuities, especially C at the start UT of the quiet day. rqdst = qdstc(ih(j)+1,i) ihm = ih(j) if (ihm .eq. 0) ihm = 24 ihp = ih(j)+2 if (ihp .eq. 25) ihp = 1 rqdstm = qdstc(ihm,i) rqdstp = qdstc(ihp,i) if (mn(j) .gt. 30) then f1 = (90.-mn(j))/60. f2 = (mn(j)-30.)/60. rqdst = f1*rqdst + f2*rqdstp else f1 = (30.-mn(j))/60. f2 = (mn(j)+30.)/60. rqdst = f1*rqdstm + f2*rqdst endif endif iqdst(i,j) = nint(rqdst) if (d(i,j) .ne. spval) h(i,j) = nint(xm + rqdst) if (ipr .eq. 1) write (6,"(1x,'final:',12x,a3,1x,a3, | 7i6)") stat(i),comp1(i),h(i,j), | d(i,j),z(i,j),hstd(i,j),dstd(i,j),zstd(i,j),iqdst(i,j) 190 continue write (22,"(3i2,i3,i2,1x,a3,1x,a3,6i6,i7,1x,3i6,1x,a3,2i6,1x, | 3i6,1x,6i6)") iy(j),im(j),id(j),ih(j),mn(j),stat(i),comp2(i), | h(i,j),d(i,j),z(i,j),hstd(i,j),dstd(i,j),zstd(i,j),iqdst(i,j), | ixmn(i,j),iymn(i,j),izmn(i,j), comp1(i),oldxd(i,j), | oldyh(i,j),(orig(k1,i,j),k1=1,3),(mean(k2,i,j),k2=1,3), | (num(k3,i,j),k3=1,3) C Eliminate if missing data IF (h(i,j) .gt. spval-1.) go to 200 C Calculate AE IF (ABS(cglat(i)) .LT. 55. .OR. ABS(cglat(i)) .GT. 76.) GO TO 195 dx = h(i,j) nae = nae + 1 AUCALC = AMAX1(DX,AUCALC) ALCALC = AMIN1(DX,ALCALC) C Calculate standard AE do kk=1,12 if(stat(i) .eq. stan(kk)) then istan = istan + 1 austan = amax1(dx,austan) alstan = amin1(dx,alstan) endif enddo C Find MLT RMLT = AMOD(cglong(i)-ROT+360.,360.) / 15. IF (RMLT .GE. 1. .AND. RMLT .LE. 21.) GO TO 195 nalm = nalm + 1 ALMIDN = AMIN1(DX,ALMIDN) 575 CONTINUE 195 CONTINUE C Calculate Dst C IF (ABS(cglat(i)) .GT. 40.) go to 200 IF (stat(i) .eq. 'AMS') go to 197 IF (stat(i) .eq. 'HER') go to 197 IF (ABS(cglat(i)) .GT. 40.) go to 200 197 continue CLATMGJ = (90. - cglat(i)) / RAD NDST = NDST + 1 DSTLON(NDST) = cglong(i) DSTVAL(NDST) = h(i,j)/SIN(CLATMGJ) IF (NDST.GT.1) THEN DO 650 K=1,NDST-1 IF (DSTLON(NDST).GE.DSTLON(K)) GO TO 650 DSTLSV = DSTLON(NDST) DSTVSV = DSTVAL(NDST) DO 630 LL=K,NDST-1 L = NDST + K - LL DSTLON(L) = DSTLON(L-1) 630 DSTVAL(L) = DSTVAL(L-1) DSTLON(K) = DSTLSV DSTVAL(K) = DSTVSV GO TO 660 650 CONTINUE 660 continue ENDIF 200 CONTINUE AECALC = AUCALC - ALCALC aestan = austan - alstan 701 IF (NDST.EQ.0) GO TO 751 DSTLON(NDST+1) = DSTLON(1) + 360. DSTLON(0) = DSTLON(NDST) - 360. DO 750 i=1,NDST 750 dstcalc = dstcalc + DSTVAL(i)*(DSTLON(i+1) - DSTLON(i-1)) dstcalc = dstcalc / 720. 751 CONTINUE C Calculate ASD (the asymmetric ring current index) asdup = -9999. asddw = 9999. do I = 1,nsta if (abs(cglat(i)) .le. 40. .and. h(i,j) .lt. 9999.) then CLATMGJ = (90. - cglat(i)) / RAD dx = h(i,j) - dstcalc*sin(CLATMGJ) asdup = AMAX1(DX,asdup) asddw = AMIN1(DX,asddw) endif enddo C Write out AE, Dst info for later IDL plotting c write (31,"(5i3,2x,i3,3f8.1,i3,f8.1,2x,i3,f8.1,2x,f8.1)") iy(j), c | im(j),id(j),ih(j),mn(j),nae,aecalc,aucalc,alcalc,nalm,almidn, c | ndst,dstcalc,dstint c write (31,"(5i3,2x,i3,3f8.1,i3,f8.1,2x,i3,f8.1,2x,f8.1,2x,f8.1)") c | iy(j),im(j),id(j),ih(j),mn(j),nae,aecalc,aucalc,alcalc,nalm, c | almidn,ndst,dstcalc,dstint,aestan write(31,"(5i3,2x,i3,3f8.1,i3,f8.1,2x,i3,f8.1,2x,f8.1,2x,3f8.1)") | iy(j),im(j),id(j),ih(j),mn(j),nae,aecalc,aucalc,alcalc,nalm, | almidn,ndst,dstcalc,dstint,aestan,asdup,asddw 900 CONTINUE END CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC SUBROUTINE ITERP(X,Y,XI,N,TERP,NN,DX1,DX2) SAVE C **** C **** PERFORMS LINEAR INTERPOLATION IN TABLE OF N (X,Y) POINTS C **** X(N) IS ARRAY OF MONOTONICALLY INCREASING ABSCISSAE C **** Y(N) IS ARRAY OF CORRESPONDING ORDINATES C **** XI IS ABSCISSA AT WHICH INTERPOLATION IS REQUIRED C **** N IS NUMBER OF POINTS IN TABLE C **** INTERPOLATED VALUE IS RETURNED IN TERP C Interpolation is performed between points NN and NN+1, where C DX1=XI-X(NN) and DX2=X(NN+1)-XI C **** C **** DIMENSION X(N),Y(N) C **** C **** SEARCH FOR INTERVAL IN WHICH XI FALLS C **** NN=1 DX1=0. DX2=0. IF(N.EQ.1)THEN TERP=Y(1) ELSE DO 1 I=1,N-1 IF(XI.GT.X(I))NN=I 1 CONTINUE C **** C **** PERFORM LINEAR INTERPOLATION C **** DX1 = XI-X(NN) DX2 = X(NN+1)-XI TERP=((X(NN+1)-XI)*Y(NN)+(XI-X(NN))*Y(NN+1))/(X(NN+1)-X(NN)) ENDIF RETURN END CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC