PROGRAM TRACE_GOES c c This program should be linked with the software package GEOPACK (version c 2005) and the model subroutine TS04 C common /data_out/ + xbt1x,xbt1y,xbt1z,xbt2x,xbt2y,xbt2z,xsrcx,xsrcy,xsrcz,xprcx, + xprcy,xprcz,xr11x,xr11y,xr11z,xr12x,xr12y,xr12z,xr21x,xr21y, + xr21z,xr22x,xr22y,xr22z COMMON /GEOPACK1/ ST0,CT0,SL0,CL0,CTCL,STCL,CTSL,STSL,SFI,CFI,SPS, * CPS,SHI,CHI,HI,PSI,XMUT,A11,A21,A31,A12,A22,A32,A13,A23,A33,DS3, * CGST,SGST,BA(6) c COMMON /GEOPACK2/G(105),H(105),REC(105) DIMENSION XC(500),YC(500),ZC(500), PARMOD(10) parameter (NPT0=25000,NPT1=500,NPT2=11000,NPT3=2800) dimension time1(NPT1),ae(NPT1),all(NPT1),dst(NPT1), + time2(NPT2),bby(NPT2),bbz(NPT2),time3(NPT3),den(NPT3), + vel(NPT3),bby_s(NPT2),bbz_s(NPT2),s_param(6,NPT0) dimension idayofmon(12),rr(6),xlamda(6),beta(6),gamma(6) character head*80,filein*50,fileout*25,satid*7 REAL X,Y,Z,BX_T04,BY_T04,BZ_T04,BX,BY,BZ REAL AMPL,R00,AL,THETA,FACTOR,FACTOR2 COMMON /PARAMETERS/ AMPL, R00, AL, THETA COMMON /RUNSCW/ ISCW COMMON /RUNTAIL/ FACTOR,FACTOR2 C FACTOR is the justification factor for the growth phase tail current modification C THETA = 0.0 for Jan 9; THETA=-0.17 for Jan 12 C ISCW = 1, or 0 for run/not run SCW.f DATA ISCW/0/ DATA re/6371.2/ DATA idayofmon/31,29,31,30,31,30,31,31,30,31,30,31/ DATA rr/0.39, 0.7, 0.031, 0.58, 1.15, 0.88/ DATA xlamda/0.39, 0.46, 0.39, 0.42, 0.41, 1.29/ DATA beta/0.80, 0.18, 2.32, 1.25, 1.6, 2.4/ DATA gamma/0.87, 0.67, 1.32, 1.29, 0.69, 0.53/ AMPL = 1.0 !used as the flag for large Pdyn mods R00 = 5.0 AL = 0.50 THETA = 0.35 THETA = 0.17 PRINT *, 'Enter the starting day & UT, and ending day & ut:' c READ *, sday,stime,eday,etime sday = 9 sut = 18. ! eday = 10 ! eut = 12. eday = 9 eut = 19. C PRINT *, ' enter ISCW: 1 or 0:' C READ *, ISCW ISCW = 0 C print *, ' choose IGRF or DIPOLE (0-dipole,>0 for IGRF)' C read *,K k = 0 C if (K .eq. 0) then C print *, ' ENTER THE GEODIPOLE TILT ANGLE IN DEGREES:' C READ *,PS C PSI=PS*0.01745329 C CPS=COS(PSI) C SPS=SIN(PSI) C endif c C print *,' Enter the GOES data file' C read *,filein filein = '/galileo/i/ganglu/d.nov04/GOES12_MAG_nov04.txt' write(satid,"(a4)") filein(27:32) filein = '/galileo/i/ganglu/d.nov04/C1_CP_FGM_SPIN_nov04.txt' write(satid,"(a7)") 'CLUSTER' print *,'satid = ',satid fileout='trace_goes12_nov04.dat' fileout='trace_C1_nov04.dat' fileout='C1_nov04.dat' c OPEN(UNIT=1,FILE=fileout,status='UNKNOWN') c open the AE and DST file open (11,file='/home/ganglu/work/aedst_nov04.dat',status='old') C do i=1,4 do i=1,10084 read(11,"(a80)") head enddo do i=1,NPT1 c read(11,"(5i3,5x,f8.1,8x,f8.1,3x,f8.1,13x,f10.1)") read(11,"(5i3,5x,f8.1,8x,f8.1,3x,f8.1,5x,f8.1)") + iyr,imo,ida,ihour,imn,ae(i),al0,all(i),dst(i) time1(i) = ida*24.+ihour + imn/60. enddo c open IMF data C add 28-imn time delay for ACE data on 9 November 2004 tdelay = 28./60. open (12,file='/galileo/d/ganglu/d.geomag/AC_MFI_nov04.txt', + status='old') do i=1,16480 read(12,"(a80)") head enddo j = 0 do i=1,NPT2 read(12,"(i2,i3,i5,2i3,f7.3,4G17.5)")ida,imo,iyr,ihour, + imn,sec,bbx0,bby0,bbz0,xdis if (bby0 .ne. -1.e+31) then j = j + 1 time2(j) = ida*24.+ihour+imn/60.+sec/3600.+tdelay bby(j) = bby0 bbz(j) = bbz0 endif enddo num_imf = j C call SMOOTH(bby,bby_s,ntp2,10) C call SMOOTH(bbz,bbz_s,ntp2,10) c open solar wind data open (13,file='/galileo/d/ganglu/d.geomag/AC_SWE_nov04.txt', + status='old') do i=1,4160 read(13,"(a80)") head enddo j = 0 do i=1,NPT3 read(13,"(2(i2,1x),i5,2i3,f7.3,8G14.5)")ida,imo,iyr,ihour, + imn,sec,den0,vel0,x_km,y_km,z_km,vx,vy,vz if (den0.ne.-1.e+31 .and. vel0.ne.-1.e+31) then j = j + 1 time3(j) = ida*24.+ihour+imn/60.+sec/3600.+tdelay C as suggest by Tsyganenko August 16, 1996 C press(i) = 1.937E-6*de*vr*vr den(j)=den0 vel(j)=vel0 endif enddo num_sw = j print *,'num_imf, num_sw = ',num_imf, num_sw C midification factors for Feb 11, 2000 C the day,ihr,imn,and delta_imn for which FACTOR changes if (iday .eq. 11) then iday_tail = 11 ihr_tail = 3 imn_tail = 50 delta_imn = 45. !for GOES 8 C duration is the interval (in imns) to keep the thin tail current duration = 90. ut_t0 = iday_tail*24.+ihr_tail+imn_tail/60. ut_t1 = ut_t0+delta_imn/60. ut_t2 = ut_t1+duration/60. endif uts = sday*24.+sut ute = eday*24.+eut ut0 = sday*24.+18. ut0 = sday*24.+22.5 ut0 = sday*24.+20.2 ! ut0 = sday*24.+12. ut0 = sday*24.+18. print *,'filein = ',filein open (9,file=filein,status='old') if (satid .eq. 'CLUSTER') goto 8 if (satid .eq. 'GE') goto 10 print *,'Reading GOES data >>>' c read (9,"(f7.1)") GEOLON c XLON = GEOLON*0.01745329 c COLAT = 90.*0.01745329 c read 1538 head lines for GOES data retrieved from CDAW database C in which B-field are in Bx, By, and Bz in GSM do ii=1,115 read (9,"(a80)") head enddo kount = 0 5 read(9,"(i2,i3,i5,i3,i3,f7.3,6G23.4)", + end=100) ida,imon,iyear,ihr,imn,sec, + hx_gsm,hy_gsm,hz_gsm,X0,Y0,Z0 bt=sqrt(hx_gsm**2+hy_gsm**2+hz_gsm**2) ut = ida*24.+ihr+imn/60.+sec/3600. if (ut .lt. ut0 .or. X0 .eq. -1.e+31) goto 5 if (ut .gt. (eday*24+24.)) go to 100 kount = kount + 1 iymd = (iyear-1900)*1000+ida if(iyear .ge. 2000) iymd = (iyear-2000)*1000+ida XGSE = X0/re YGSE = Y0/re ZGSE = Z0/re goto 15 8 continue print *,'Reading CLUSTER data >>>' kount = 0 do ii=1,125 read (9,"(a80)") head enddo 9 read(9,"(i2,i3,i5,i3,i3,f7.3,4G14.4,3g17.4)", + end=100) ida,imon,iyear,ihr,imn,sec, + bt,bx_gse,by_gse,bz_gse,X0,Y0,Z0 bt=sqrt(bx_gse**2+by_gse**2+bz_gse**2) ut = ida*24.+ihr+imn/60.+sec/3600. if (ut .lt. ut0 .or. X0 .eq. -1.e+31) goto 9 if (ut .gt. (eday*24+eut)) go to 100 kount = kount + 1 iymd = (iyear-1900)*1000+ida if(iyear .ge. 2000) iymd = (iyear-2000)*1000+ida XGSE = X0/re YGSE = Y0/re ZGSE = Z0/re goto 15 10 continue print *,'Reading GEOTAIL data >>>' C Read GEOTAIL data do ii=1,38 read (9,"(a80)") head enddo kount = 0 11 continue do ii=1,20 read(9,"(i2,i3,i5,i3,i3,f7.3,3f15.5,3f14.4)", + end=100) ida,imon,iyear,ihr,imn,sec, + X0,Y0,Z0,hx_gsm,hy_gsm,hz_gsm enddo ht=sqrt(hx_gsm**2+hy_gsm**2+hz_gsm**2) ut = ida*24.+ihr+imn/60.+sec/3600. if (ut .lt. ut0 .or. X0 .eq. -1.e+31) go to 11 if (ut .gt. (eday*24+24.)) go to 100 cc print *, ida,ihr,imn,sec, cc + hx_gse,hy_gse,hz_gse,X0,Y0,Z0 kount = kount + 1 iymd = (iyear-1900)*1000+ida if(iyear .ge. 2000) iymd = (iyear-2000)*1000+ida XGSE = X0 YGSE = Y0 ZGSE = Z0 15 continue jday = 0 do i=1, imon-1 jday = jday + idayofmon(i) enddo jday = jday + ida if (jday .ne. 314) AMPL = 0.0 !used as the flag for large Pdyn mods CALL RECALC(iyear,jday,ihr,imn,0) C if (kount .eq. 1) print *,'iyear,PSI,CHI,SHI,SPS,CPS = ', + iyear,jday,ihr,imn,sec,PSI*57.29578,CHI,SHI,SPS,CPS C print *,'iyear,jday,ihr,imn = ',iyear,jday,ihr,imn call GSMGSE(X,Y,Z,XGSE,YGSE,ZGSE,-1) if (X0 .eq. -1.e+31) X = 99999. if (Y0 .eq. -1.e+31) Y = 99999. if (Z0 .eq. -1.e+31) Z = 99999. if (kount .eq. 1) print *,'X,Y,Z,XGSE,YGSE,ZGSE = ', | x,y,z,xgse,ygse,zgse c C add a time delay for SCW to set up del_t = atan(abs(Y/X)) ut1 = ut - del_t*57.3/120. c if (ut .gt. 294 .and. ut.lt. 296) c + print *,'del_t,ut,ut1 = ',del_t*57.3,ut,ut1 c interpolate AE,DST, Bx, and By values c call ITERP(time3,press,ut,num_sw,PARMOD(1),nn3,x1,x2) call ITERP(time3,den,ut,num_sw,den_fit,nn3,x1,x2) call ITERP(time3,vel,ut,num_sw,vel_fit,nn3,x1,x2) PARMOD(1) = 0.167e-5*den_fit*vel_fit**2 call ITERP(time1,dst,ut,NPT1,PARMOD(2),nn1,x1,x2) call ITERP(time2,bby,ut,num_imf,PARMOD(3),nn2,x1,x2) call ITERP(time2,bbz,ut,num_imf,PARMOD(4),nn2,x1,x2) print *,'time3, ut = ',time3(1),time3(num_sw),time3(nn3), ut bs=0. if (PARMOD(4) .lt. 0.) bs=abs(PARMOD(4)) do ipar=1,6 if (kount .eq. 1 .or. ut.lt.uts .or. ut.gt.ute) then s_param(ipar,kount) = | rr(ipar)/12.*(den_fit/5.)**xlamda(ipar)* | (vel_fit/400.)**beta(ipar)*(bs/5)**gamma(ipar) else s_param(ipar,kount) = s_param(ipar,kount-1) + | rr(ipar)/12.*(den_fit/5.)**xlamda(ipar)* | (vel_fit/400.)**beta(ipar)* | (bs/5)**gamma(ipar)*exp(rr(ipar)*(uts-ut)) endif CCC Note: ut and ut0 are in hour, not minute, as written in TS, JGR, 2005 PARMOD(ipar+4) = (rr(ipar)/12.)*s_param(ipar,kount) c print *,'ipar,den_fit,vel_fit,xlamda,beta,gamma =',ipar, c | den_fit,vel_fit,xlamda(ipar),beta(ipar), c | gamma(ipar),rr(ipar),s_param(ipar,kount),parmod(ipar+4) enddo c print *,'PARMOD(1),bs,PARMOD(5:10) = ',PARMOD(1),bs, c | PARMOD(5),PARMOD(6),PARMOD(7),PARMOD(8),PARMOD(9),PARMOD(10) c if (kount .eq. 10) STOP c if (PARMOD(2) .gt. -10.) PARMOD(2) = -10. C ***T96_01 does not use AE and have take that out in T96_01 code*** IOPT = 0 write(6,"('Main call T04_s: ampl=',e12.4,' r00=',e12.4,' al=', | e12.4,' theta=',e12.4)") ampl,r00,al,theta CALL T04_s(IOPT,PARMOD,PSI,X,Y,Z,BX_T04,BY_T04,BZ_T04) c CALL T04_s(IOPT,PARMOD,PS,X,Y,Z,BX_T04,BY_T04,BZ_T04,xbt1x, c + xbt1y,xbt1z,xbt2x,xbt2y,xbt2z,xsrcx,xsrcy,xsrcz,xprcx, c + xprcy,xprcz,xr11x,xr11y,xr11z,xr12x,xr12y,xr12z,xr21x,xr21y, c + xr21z,xr22x,xr22y,xr22z) c print *,'X,Y,Z,BX,BY,BZ = ',x,y,z,bx_t04,by_t04,bz_t04 c CALL T95_06(IOPT,PARMOD,PSI,X,Y,Z,BX_T95,BY_T95,BZ_T95) c CALL T89M_KP(IOPT,PARMOD,PSI,X,Y,Z,BX_T89,BY_T89,BZ_T89) c CALL T89(IOPT,PSI,X,Y,Z,BX,BY,BZ) BX = BX_T04 BY = BY_T04 BZ = BZ_T04 C calculate SCW current induced B-field BXSCW = 0. BYSCW = 0. BZSCW = 0. if (ISCW .ge. 1) THEN SPS=SIN(PSI) C print *, 'PSI,SPS,X,Y,Z = ',PSI, SPS,X,Y,Z C CALL WARPWEDG(X,Y,Z,SPS,BXSCW,BYSCW,BZSCW) C CALL WARPWEDG(dble(X),dble(Y),dble(Z),dble(SPS),BXD,BYD,BZD) C print *, 'PSI,SPS,X,Y,Z,BXD,BYD,BZD = ', C + PSI, SPS,X,Y,Z,BXD,BYD,BZD BXSCW = real(BXD) BYSCW = real(BYD) BZSCW = real(BZD) endif k=6 IF (K.EQ.0) GOTO 1 CALL IGRF_GSM(X,Y,Z,HX,HY,HZ) GOTO 2 1 CALL DIP(X,Y,Z,HX,HY,HZ) 2 BX=BX+HX BY=BY+HY BZ=BZ+HZ C using only the internal field *** temp C 2 BX=HX C BY=HY C BZ=HZ C using only the SCW *** temp c BX=HX+BXSCW c BY=HY+BYSCW c BZ=HZ+BZSCW print *,'TRACE: XGSE,YGSE,ZGSE,X,Y,Z,HX,HY,HZ = ', + XGSE,YGSE,ZGSE,x,y,z,HX,HY,HZ C convert from GSM to GSE CALL GSMGSE(BX,BY,BZ,BXGSE,BYGSE,BZGSE,1) C convert from GSM to GEO c CALL GEOGSM(BXGEO,BYGEO,BZGEO,BX,BY,BZ,-1) C convert from GEO to cylindrical coordinates c BRR = BXGEO*COS(XLON) + BYGEO*SIN(XLON) c BEE = -BXGEO*SIN(XLON) + BYGEO*COS(XLON) c c write(1,*) iymd,ihr,imn,iday,hp,he,hn,ht, c + BZGEO,-BRR,BEE,X,Y,Z c write(1,*) iymd,ihr,imn,iday,hp,he,hn,ht, if (satid .eq. 'GOES') then write(1,*) iyear,jday,ihr,imn,sec,hx_gsm,hy_gsm,hz_gsm,bt, C + BX,BY,BZ,HX,HY,HZ,X,Y,Z,xgse,ygse,zgse C + BX_T04,BY_T04,BZ_T04,HX,HY,HZ,xbt1x,xbt1y,xbt1z, + BX,BY,BZ,HX,HY,HZ,XGSE,YGSE,ZGSE c + BXGSE,BYGSE,BZGSE,BX_T04,BY_T04,BZ_T04,xbt1x,xbt1y,xbt1z, c + xbt2x,xbt2y,xbt2z,xsrcx,xsrcy,xsrcz,xprcx,xprcy,xprcz, c + xr11x,xr11y,xr11z,xr21x,xr21y,xr21z, c + xbxr12,xbyr12,xbzr12,xr22x,xr22y,xr22z endif if (satid .eq. 'CLUSTER') then write(1,*) iyear,jday,ihr,imn,sec,bx_gse,by_gse,bz_gse,bt, C + BX,BY,BZ,HX,HY,HZ,X,Y,Z,xgse,ygse,zgse C + BX_T04,BY_T04,BZ_T04,HX,HY,HZ,xbt1x,xbt1y,xbt1z, + BXGSE,BYGSE,BZGSE,HX,HY,HZ,XGSE,YGSE,ZGSE c + BXGSE,BYGSE,BZGSE,BX_T04,BY_T04,BZ_T04,xbt1x,xbt1y,xbt1z, c + xbt2x,xbt2y,xbt2z,xsrcx,xsrcy,xsrcz,xprcx,xprcy,xprcz, c + xr11x,xr11y,xr11z,xr21x,xr21y,xr21z, c + xbxr12,xbyr12,xbzr12,xr22x,xr22y,xr22z endif print *,'satid = ',satid if (satid .eq. 'GE') goto 11 if (satid .eq. 'CLUSTER') goto 9 goto 5 C 100 stop END CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC SUBROUTINE ITERP(X,Y,XI,N,TERP,NN,DX1,DX2) SAVE 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 DIMENSION X(N),Y(N) 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 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC SUBROUTINE SMOOTH(X,Y,N,M) SAVE C **** PERFORMS BOX-CAR AVERAGE OF N (X,Y) POINTS C **** N IS NUMBER OF POINTS IN TABLE C **** M is the number of points to be averaged over DIMENSION X(N),Y(N) diff = 0. DO I=1,N YY = 0. JJ = 0 IF (I .LT. M/2) THEN DO J=1,M YY = YY + X(J) JJ = JJ +1 ENDDO ELSE IF (I .GT. N-M/2) THEN DO J=I-M/2,N YY = YY + X(J) JJ = JJ +1 ENDDO ELSE DO J=I-M/2,I+M/2 YY = YY + X(J) JJ = JJ +1 ENDDO ENDIF ENDIF Y(I) = YY/JJ diff = diff + abs(X(I)-Y(I)) ENDDO c print *,'Call SMOOTH: N,M,JJ,diff = ',N,M,JJ,diff RETURN END CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC