C&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
C                                                                          &
C          SUPPLEMENTARY CODES:   LOCATE  and   CROSSING                   &
C                                                                          &
C&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&



      SUBROUTINE LOCATE 
     *(XN_PD,VEL,XGSM,YGSM,ZGSM,XMGNP,YMGNP,ZMGNP,DIST,ID)
C
C     THIS SUBROUTINE DEFINES THE POSITION OF A POINT (XMGNP,YMGNP,ZMGNP)
C            AT THE MODEL MAGNETOPAUSE, CLOSEST TO A GIVEN POINT OF SPACE
C               (XGSM,YGSM,ZGSM),   AND THE DISTANCE BETWEEN THEM (DIST)
C
C INPUT:  XN_PD - EITHER SOLAR WIND PROTON NUMBER DENSITY (PER C.C.) (IF VEL>0)
C                    OR THE SOLAR WIND RAM PRESSURE IN NANOPASCALS   (IF VEL<0)
C         VEL - EITHER SOLAR WIND VELOCITY (KM/SEC) 
C                  OR ANY NEGATIVE NUMBER, WHICH INDICATES THAT XN_PD STANDS
C                     FOR THE SOLAR WIND PRESSURE, RATHER THAN FOR THE DENSITY
C
C         XGSM,YGSM,ZGSM - COORDINATES OF THE OBSERVATION POINT IN EARTH RADII
C
C  OUTPUT: XMGNP,YMGNP,ZMGNP - COORDINATES OF A POINT AT THE MAGNETOPAUSE,
C                                    CLOSEST TO THE POINT  XGSM,YGSM,ZGSM
C          DIST -  THE DISTANCE BETWEEN THE ABOVE TWO POINTS, IN RE,
C          ID -    INDICATOR; ID=+1 AND ID=-1 MEAN THAT THE POINT 
C                       (XGSM,YGSM,ZGSM)  LIES INSIDE OR OUTSIDE 
C                        THE MODEL MAGNETOPAUSE, RESPECTIVELY
C
C   THE PRESSURE-DEPENDENT MAGNETOPAUSE IS THAT USED IN THE T96_01 MODEL
C
c            CODED BY:  N.A. TSYGANENKO, AUG.1, 1995;  REVISED  JUNE 22, 1996
C
      REAL XN_PD ,VEL   ,XGSM  ,
     *     YGSM  ,ZGSM  ,XMGNP ,
     *     YMGNP ,ZMGNP ,DIST  
      INTEGER ID 

      IF (VEL.LT.0.) THEN 
       PD=XN_PD
                  ELSE
       PD=1.94E-6*XN_PD*VEL**2  ! PD IS THE SOLAR WIND DYNAMIC PRESSURE 
C                                      (IN NANOPASCALS)
      ENDIF
C
      RAT=PD/2.0  ! RATIO OF PD TO THE AVERAGE PRESSURE, ASSUMED AS 2 nPa
C
      RAT16=RAT**0.14 ! THE POWER IN THE SCALING FACTOR IS THE BEST-FIT VALUE
C                         OBTAINED FROM DATA IN THE T96_01 VERSION OF THE MODEL
C
      A0=70.
      S00=1.08
      X00=5.48    !  VALUES OF THE MAGNETOPAUSE PARAMETERS FOR  PD = 2 nPa
C
      A=A0/RAT16
      S0=S00
      X0=X00/RAT16   !  VALUES OF THE MAGNETOPAUSE PARAMETERS, SCALED TO THE 
C                         ACTUAL PRESSURE
C
       XM=X0-A    !  THIS IS THE X-COORDINATE OF THE "SEAM" BETWEEN THE
C                             ELLIPSOID AND THE CYLINDER
C
C        (FOR DETAILS ON THE ELLIPSOIDAL COORDINATES, SEE THE PAPER:  
C            N.A.TSYGANENKO, SOLUTION OF CHAPMAN-FERRARO PROBLEM FOR AN
C             ELLIPSOIDAL MAGNETOPAUSE, PLANET.SPACE SCI., V.37, P.1037, 1989).
C
          IF (YGSM.NE.0..OR.ZGSM.NE.0.) THEN 
             PHI=ATAN2(YGSM,ZGSM)
              ELSE
             PHI=0.
          ENDIF
C           
          RHO=SQRT(YGSM**2+ZGSM**2)
C
         IF (XGSM.LT.XM) THEN
           XMGNP=XGSM
           RHOMGNP=A*SQRT(S0**2-1)
           YMGNP=RHOMGNP*SIN(PHI)
           ZMGNP=RHOMGNP*COS(PHI)
           DIST=SQRT((XGSM-XMGNP)**2+(YGSM-YMGNP)**2+(ZGSM-ZMGNP)**2)
           IF (RHOMGNP.GT.RHO) ID=+1
           IF (RHOMGNP.LT.RHO) ID=-1
           RETURN
         ENDIF
C
          XKSI=(XGSM-X0)/A+1.
          XDZT=RHO/A
          SQ1=SQRT((1.+XKSI)**2+XDZT**2)
          SQ2=SQRT((1.-XKSI)**2+XDZT**2)
          SIGMA=0.5*(SQ1+SQ2)
          TAU=0.5*(SQ1-SQ2)
C
C  NOW CALCULATE (X,Y,Z) FOR THE CLOSEST POINT AT THE MAGNETOPAUSE
C
          XMGNP=X0-A*(1.-S0*TAU)
          RHOMGNP=A*SQRT((S0**2-1.)*(1.-TAU**2))
          YMGNP=RHOMGNP*SIN(PHI)
          ZMGNP=RHOMGNP*COS(PHI)
C
C  NOW CALCULATE THE SHORTEST DISTANCE BETWEEN THE POINT XGSM,YGSM,ZGSM AND THE
C            MAGNETOPAUSE 
C
      DIST=SQRT((XGSM-XMGNP)**2+(YGSM-YMGNP)**2+(ZGSM-ZMGNP)**2)
C
      IF (SIGMA.GT.S0) ID=-1   !  ID=-1 MEANS THAT THE POINT LIES OUTSIDE
      IF (SIGMA.LT.S0) ID=+1   !  ID=+1 MEANS THAT THE POINT LIES INSIDE
C                                           THE MAGNETOSPHERE
      RETURN
      END
C
C
C----------------------------------------------------------------------
C
      SUBROUTINE CROSSING(XGSM,YGSM,ZGSM,PDYN)
C
C  THIS SUBROUTINE ACCEPTS THE GSM COORDINATES OF AN OBSERVED MAGNETOPAUSE
C    CROSSING AND RETURNS AN APPROXIMATE ESTIMATE OF THE SOLAR WIND DYNAMIC
C    PRESSURE, BASED ON THE ELLIPSOIDAL APPROXIMATION FOR THE BOUNDARY AND
C    ASSUMING A SELF-SIMILAR SCALING OF THE MAGNETOSPHERIC DIMENSIONS IN RES-
C    PONSE TO CHANGES IN THE RAM PRESSURE.
C
C    INPUT:   XGSM,YGSM,ZGSM - COORDINATES OF THE CROSSING POINT IN EARTH RADII
C    OUTPUT:  PDYN - THE ESTIMATE OF THE SOLAR WIND PRESSURE IN NANOPASCALS
C
C       THE PRESSURE-DEPENDENT MAGNETOPAUSE IS THAT USED IN THE T95_01 MODEL
C          (TSYGANENKO, 1995)
C
C     WARNING:  THIS SUBROUTINE SHOULD NOT BE USED TAILWARD OF XGSM=-60 RE !!!
C
c            CODED BY:  N.A. TSYGANENKO, AUG.2, 1995.
C
      PD0=2.
      A=70.
      S0=1.08
      X0=5.48    !  VALUES OF THE MAGNETOPAUSE PARAMETERS FOR  PD = 2 nPa
C
C
C        (FOR DETAILS ON THE ELLIPSOIDAL COORDINATES, SEE THE PAPER:  
C            N.A.TSYGANENKO, SOLUTION OF CHAPMAN-FERRARO PROBLEM FOR AN
C             ELLIPSOIDAL MAGNETOPAUSE, PLANET.SPACE SCI., V.37, P.1037, 1989).
C
          AMX0=A-X0
          RHO2=YGSM**2+ZGSM**2
          S02=S0**2
C
          EPS=(XGSM*AMX0+S0*SQRT(A**2*(XGSM**2+RHO2*S02/(S02-1.))
     *           -AMX0**2*RHO2/(S02-1.)))/((A*S0-AMX0)*(A*S0+AMX0))
          PDYN=PD0/EPS**(1./0.14)  ! THE POWER IN THE SCALING FACTOR CORRES-
C                                    PONDS TO THE BEST-FIT VALUE OBTAINED FROM
C                                    DATA IN THE T96_01 VERSION OF THE MODEL.
C
         RETURN
         END

C
C~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
      SUBROUTINE DIPOLE (PS,X,Y,Z,BX,BY,BZ)
C
C     THIS IS A DOUBLE PRECISION ROUTINE, OTHERWISE IDENTICAL TO THE S/R DIP OF GEOPACK
C
C     CALCULATES GSM COMPONENTS OF A GEODIPOLE FIELD WITH THE DIPOLE MOMENT
C     CORRESPONDING TO THE EPOCH OF 2000.
C
C------INPUT PARAMETERS:
C       PS - GEODIPOLE TILT ANGLE IN RADIANS,
C       X,Y,Z - GSM COORDINATES IN RE (1 RE = 6371.2 km)
C
C----OUTPUT PARAMETERS:
C     BX,BY,BZ - FIELD COMPONENTS IN GSM SYSTEM, IN NANOTESLA.
C
C    LAST MODIFICATION: JAN. 5, 2001. THE VALUE OF THE DIPOLE MOMENT WAS UPDATED TO 2000.
C      AND A "SAVE" STATEMENT HAS BEEN ADDED, TO AVOID POTENTIAL PROBLEMS WITH SOME
C      FORTRAN COMPILERS
C
C    WRITTEN BY: N. A. TSYGANENKO
C
      IMPLICIT REAL*8 (A-H,O-Z)
      SAVE M,PSI
      DATA M,PSI/0,5.D0/
      IF(M.EQ.1.AND.DABS(PS-PSI).LT.1.D-5) GOTO 1   !   THIS IS TO AVOID MULTIPLE CALCULATIONS
      SPS=DSIN(PS)                                  !   OF SIN(PS) AND COS(PS), IF THE ANGLE PS
      CPS=DCOS(PS)                                  !   REMAINS UNCHANGED
      PSI=PS
      M=1
  1   P=X**2
      U=Z**2
      V=3.D0*Z*X
      T=Y**2
      Q=30115.D0/DSQRT(P+T+U)**5
      BX=Q*((T+U-2.D0*P)*SPS-V*CPS)
      BY=-3.D0*Y*Q*(X*SPS+Z*CPS)
      BZ=Q*((P+T-2.D0*U)*CPS-V*SPS)
      RETURN
      END

	SUBROUTINE GSMMAG(xgsm,ygsm,zgsm,r,lat,MLT,j)
c	Written by JM, April 1998, corrected Sep 2000
c	lat is in degrees and MLT in hours
c     pridavam promenou rsun  
	IMPLICIT NONE
	REAL xgsm,ygsm,zgsm,r,lat,pi,mlt,x,y,z
	REAL junk(3),teta,phi
	REAL radeg,msun,rsun
	INTEGER j
	DATA pi /3.141592654/

	radeg = 1./ATAN(1.)*45.

c	 poloha osy x (tj slunce) v MAG
	call SMGSM(junk(1),junk(2),junk(3),10.,0.,0.,-1) 
	call MAGSM(x,y,z,junk(1),junk(2),junk(3),-1)
	call SPHCAR(rsun,teta,msun,x,y,z,-1)
	IF (j.GT.0) THEN 
	 call SMGSM(junk(1),junk(2),junk(3),xgsm,ygsm,zgsm,-1)
	 call MAGSM(x,y,z,junk(1),junk(2),junk(3),-1)
	 call SPHCAR(r,teta,phi,x,y,z,-1)
	 lat=(pi/2-teta)*radeg
	 mlt=(phi-msun)*radeg/15.+12.
	 if (mlt.GE.24) mlt=mlt-24
	 if (mlt.LT.0)  mlt=mlt+24
c	call SPHCAR(r,teta,phi,junk(1),junk(2),junk(3),-1)
c	lat=(pi/2-teta)*radeg
c	mlt=phi*radeg/15.-12.
c	IF (mlt.LT.0.) mlt=mlt+24.
	ELSE
	 phi=(mlt-12.)*15./radeg+msun
	 teta=pi/2-lat/radeg
	 call SPHCAR(r,teta,phi,x,y,z,1)
	 call MAGSM(x,y,z,junk(1),junk(2),junk(3),1)
	 call SMGSM(junk(1),junk(2),junk(3),xgsm,ygsm,zgsm,1)
	END IF
	RETURN
      END
