C*
      SUBROUTINE ERPFBOXW(ERM,ANT,GRD,REFF,YSAT,SUN,KAPPA,MONTH,
     1                    BLKNUM,SVN,MJD,ACCEL)
CC
CC NAME       :  ERPFBOXW
CC
CC PURPOSE    :  COMPUTATION OF EARTH RADIATION PRESSURE ACTING ON A
CC               BOW-WING SATELLITE
CC
CC PARAMETERS :
CC         IN :  ERM       : EARTH RADIATION MODEL
CC                           0 = NONE
CC                           1 = EARTH RADIATION PRESSURE (ANALYTICAL)
CC                           2 = EARTH RADIATION PRESSURE (CERES DATA)
CC               CERES DAT : CERES DATA (EXTERNAL ASCII FILES) -> SET PATH IN DATPATH
CC               ANT       : 0 = NO ANTENNA THRUST
CC                         : 1 = WITH ANTENNA THRUST
CC               GRD       : 1 = 2.5 degree grid
CC                         : 2 = 5.0 degree grid
CC                         : 3 = 10  degree grid
CC               REFF      : ACCELERATION IN REFERENCE FRAME:
CC                           0 = INERTIAL
CC                           1 = BODY FIXED (Z,Y,X)
CC                           2 = SUN FIXED (D,Y,B)
CC                           3 = ORBITAL (RADIAL, ALONG- AND CROSS-TRACK)
CC               YSAT      : SATELLITE POSITION [m] AND VELOCITY [m/s] (INERTIAL), 
CC                           (POSITION,VELOCITY) = (RX,RY,RZ,VX,VY,VZ)
CC               SUN       : SUN POSITION VECTOR [m] (INERTIAL)
CC               KAPPA     : ROTATION MATRIX FROM INERTIAL TO EARTH-FIXED
CC                           REFERENCE FRAME (FOR ERM=2)
CC               MONTH     : MONTH NUMBER (1 ... 12), FOR ERM=2
CC               BLKNUM    : BLOCK NUMBER
CC                             1 = GPS-I
CC                             2 = GPS-II
CC                             3 = GPS-IIA
CC                             4 = GPS-IIR
CC                             5 = GPS-IIR-A
CC                             6 = GPS-IIR-B
CC                             7 = GPS-IIR-M
CC                             8 = GPS-IIF
CC                           101 = GLONASS
CC                           102 = GLONASS-M
CC                           103 = GLONASS-K (not yet available)
CC               SVN       : SPACE VEHICLE NUMBER          
CC               MJD       : MODIFIED JULIAN DAY
CC
CC        OUT : ACCEL      : ACCELERATION VECTOR [m/s^2]
CC
CC AUTHOR     : C.J. RODRIGUEZ-SOLANO
CC              rodriguez@bv.tum.de
CC
CC VERSION    : 1.0 (OCT 2010)
CC
CC CREATED    : 2010/10/18
CC
C*
C
      IMPLICIT NONE
C
      INTEGER*4 BLKNUM,ERM,ANT,GRD,REFF,SBLK,INDB,MONTH,SVN
      INTEGER*4 IFIRST,II,JJ,K,LFNLOC,IOSTAT,LENPATH
      INTEGER*4 LATK,LONK,LATKMX,LONKMX,GRDCER
C
      REAL*8 PI,C,AU,S0,TOA,ALB,MASS
      REAL*8 YSAT(6),SUN(3),FORCE(3),ACCEL(3)
C
C     EARTH AND SATELLITE PROPERTIES
      REAL*8 CERES_R(72,144),CERES_E(72,144)
      REAL*8 CERGRE(72,144),CERGEM(72,144)
      REAL*8 D_AREA_ALL(72,144),V_NS_ALL(72,144,3)
      REAL*8 AREA(4,2,15),REFL(4,2,15),DIFU(4,2,15),ABSP(4,2,15)
      REAL*8 AREA2(4,2,15),REFL2(4,2,15),DIFU2(4,2,15),ABSP2(4,2,15)
      REAL*8 REFLIR(4,2,15),DIFUIR(4,2,15),ABSPIR(4,2,15)
      REAL*8 AREAS(4,2),REFLS(4,2),DIFUS(4,2),ABSPS(4,2)
      REAL*8 AREA2S(4,2),REFL2S(4,2),DIFU2S(4,2),ABSP2S(4,2)
      REAL*8 REFLIRS(4,2),DIFUIRS(4,2),ABSPIRS(4,2)
C
C     ATTITUDE
      REAL*8 RADVEC(3),ALGVEC(3),CRSVEC(3),ESUN(3)
      REAL*8 RSUN,ABSPOS,ABSCRS,ABSY0V,Y0SAT(3)
      REAL*8 Z_SAT(3),D_SUN(3),Y_SAT(3),B_SUN(3),X_SAT(3)
      REAL*8 ATTSURF(3,4)
C
      REAL*8 ABSNCFVI,ABSNCFIR,ALBFAC,PHASEVI,PHASEIR
      REAL*8 NCFVEC(3)
      REAL*8 D_ANG,GRDANG,GRDNUM,LATIND,LONIND
      REAL*8 D_AREA,LAT_IN,LON_IN,DIST2,ABSDIST
      REAL*8 COSLAT,PSI
      REAL*8 V_NS(3),V_INS(3),V_DIST(3),V_SAT(3)
      REAL*8 COS_IN,COS_RE,PSIDOT,ABSSUN
      REAL*8 REFL_CF,EMIT_CF,E_REFL,E_EMIT
      REAL*8 FORCE_LL(3),FREF(3)
      REAL*8 ANTFORCE,ANTPOW,MJD
      REAL*8 KAPPA(3,3)
C
      CHARACTER*2 F_MONTH(12)
      CHARACTER*100 DATPATH,FILREFL,FILEMIT
      CHARACTER*5000 LINE
      CHARACTER*25 ITEM
C
      DATA IFIRST/1/
      DATA (F_MONTH(K),K=1,12)/
     . '01', '02', '03', '04', '05', '06', 
     . '07', '08', '09', '10', '11', '12'/

C DATA FILES PATH AND NAME (CHANGE TO LOCAL PATH)
      DATPATH = './CERES/'

C COMPLETE FILE NAMES
      II = LEN(DATPATH)
      DO WHILE (DATPATH(II:II).EQ.' ')
         II = II-1
      ENDDO
      LENPATH = II

      FILREFL = DATPATH(1:LENPATH) // 'REFLMO' // F_MONTH(MONTH)
      FILEMIT = DATPATH(1:LENPATH) // 'EMITMO' // F_MONTH(MONTH)

C Constants needed
      PI = 4D0*DATAN(1D0)
      C  = 299792458D0
      AU = 149597870691D0
C Solar Constant(W/m2)
      S0 = 1367D0
C Top of Atmosphere for CERES
      TOA = 6371000D0 + 30000D0
C Albedo of the Earth
      ALB = 0.3D0
C Initialization of force vector
      FORCE(1) = 0D0
      FORCE(2) = 0D0
      FORCE(3) = 0D0

      ACCEL(1) = 0D0
      ACCEL(2) = 0D0
      ACCEL(3) = 0D0

      LFNLOC = 999

C ----------------------------------------
C LOAD SATELLITE PROPERTIES AND CERES DATA
C ----------------------------------------
      IF ((IFIRST==1).AND.(ERM.GT.0)) THEN

C     PROPERTIES FOR ALL SATELLITES BLOCKS
         DO SBLK = 1,105
            IF((SBLK.LE.10).OR.(SBLK.GT.100))THEN 

            CALL PROPBOXW(SBLK,AREAS,REFLS,DIFUS,ABSPS,AREA2S,REFL2S,
     1                    DIFU2S,ABSP2S,REFLIRS,DIFUIRS,ABSPIRS)
            
            IF(SBLK.LE.10)THEN
               INDB = SBLK
            ELSEIF(SBLK.GT.100)THEN
               INDB = SBLK - 90
            ENDIF

            DO II = 1,4
               DO JJ = 1,2
                  AREA(II,JJ,INDB) = AREAS(II,JJ)
                  REFL(II,JJ,INDB) = REFLS(II,JJ)
                  DIFU(II,JJ,INDB) = DIFUS(II,JJ)
                  ABSP(II,JJ,INDB) = ABSPS(II,JJ)

                  AREA2(II,JJ,INDB) = AREA2S(II,JJ)
                  REFL2(II,JJ,INDB) = REFL2S(II,JJ)
                  DIFU2(II,JJ,INDB) = DIFU2S(II,JJ)
                  ABSP2(II,JJ,INDB) = ABSP2S(II,JJ)

                  REFLIR(II,JJ,INDB) = REFLIRS(II,JJ)
                  DIFUIR(II,JJ,INDB) = DIFUIRS(II,JJ)
                  ABSPIR(II,JJ,INDB) = ABSPIRS(II,JJ)
               ENDDO
            ENDDO
            ENDIF
         ENDDO


         IF(ERM.EQ.2)THEN




C     REFLECTIVITY 
            CERES_R = 0
            OPEN(UNIT=LFNLOC,FILE=FILREFL,STATUS='UNKNOWN',
     1           FORM='FORMATTED',IOSTAT=IOSTAT)
            DO II=1,72
               READ(LFNLOC,"(A)")LINE
               DO JJ=1,144
                  ITEM = LINE((JJ-1)*25+1:JJ*25)
                  IF (INDEX(ITEM,'NaN') /= 0) THEN
                     CERES_R(II,JJ) = 0D0
                  ELSE
                     READ(ITEM,*) CERES_R(II,JJ)
                  ENDIF
               ENDDO
            ENDDO
            CLOSE(LFNLOC)
C     EMISSIVITY 
            CERES_E = 0
            OPEN(UNIT=LFNLOC,FILE=FILEMIT,STATUS='UNKNOWN',
     1           FORM='FORMATTED',IOSTAT=IOSTAT)
            DO II=1,72
               READ(LFNLOC,"(A)")LINE
               DO JJ=1,144
                  ITEM = LINE((JJ-1)*25+1:JJ*25)
                  IF (INDEX(ITEM,'NaN') /= 0) THEN
                     CERES_E(II,JJ) = 0D0
                  ELSE
                     READ(ITEM,*) CERES_E(II,JJ)
                  ENDIF
               ENDDO
            ENDDO
            CLOSE(LFNLOC)

C     PRE-INTEGRATION, INDEPENDENT OF SATELLITE POSITION
            GRDANG = 2.5D0
            LATIND = 36.5D0
            LONIND = 72.5D0
            IF(GRD.EQ.1)THEN
               GRDNUM = 1D0
               LATKMX = 72
               LONKMX = 144
            ELSEIF(GRD.EQ.2)THEN
               GRDNUM = 2D0
               GRDCER = 2
               LATKMX = 36
               LONKMX = 72
            ELSEIF(GRD.EQ.3)THEN
               GRDNUM = 4D0
               GRDCER = 4
               LATKMX = 18
               LONKMX = 36
            ENDIF
 
            GRDANG = GRDANG*GRDNUM
            LATIND = (LATIND-0.5D0)/GRDNUM + 0.5D0
            LONIND = (LONIND-0.5D0)/GRDNUM + 0.5D0
              
            D_ANG = (PI*GRDANG/180D0)**2

            DO LATK = 1,LATKMX
               DO LONK = 1,LONKMX

                  LAT_IN = (LATK-LATIND)*GRDANG*(PI/180D0)
                  LON_IN = (LONK-LONIND)*GRDANG*(PI/180D0)

C                 Sphere normal vector and differential of area
                  COSLAT = DCOS(LAT_IN)
                  D_AREA_ALL(LATK,LONK) = (TOA**2)*COSLAT*D_ANG
                  V_NS_ALL(LATK,LONK,1) = COSLAT*DCOS(LON_IN)
                  V_NS_ALL(LATK,LONK,2) = COSLAT*DSIN(LON_IN)
                  V_NS_ALL(LATK,LONK,3) = DSIN(LAT_IN)

C                 New matrix of Reflectivity and Emissivity
                  CERGRE(LATK,LONK) = 0D0
                  CERGEM(LATK,LONK) = 0D0
                  IF(GRD.EQ.1)THEN
                     CERGRE(LATK,LONK) = CERES_R(LATK,LONK)
                     CERGEM(LATK,LONK) = CERES_E(LATK,LONK)
                  ELSEIF((GRD.EQ.2).OR.(GRD.EQ.3))THEN
                     DO II = 0,(GRDCER-1)
                        DO JJ = 0,(GRDCER-1)
                           CERGRE(LATK,LONK) = CERGRE(LATK,LONK) 
     1                     + CERES_R(GRDCER*LATK-II,GRDCER*LONK-JJ)
                           CERGEM(LATK,LONK) = CERGEM(LATK,LONK)
     1                     + CERES_E(GRDCER*LATK-II,GRDCER*LONK-JJ)
                        ENDDO
                     ENDDO
                     CERGRE(LATK,LONK) = CERGRE(LATK,LONK)/(GRDNUM**2)
                     CERGEM(LATK,LONK) = CERGEM(LATK,LONK)/(GRDNUM**2)
                  ENDIF

               ENDDO
            ENDDO
C
          ENDIF
         IFIRST = 0
      ENDIF


C --------------------------
C NOMINAL SATELLITE ATTITUDE
C --------------------------

      ABSPOS = DSQRT(YSAT(1)**2+YSAT(2)**2+YSAT(3)**2)
      DO K=1,3
         RADVEC(K) = YSAT(K)/ABSPOS
      ENDDO
      
      CRSVEC(1) = YSAT(2)*YSAT(6)-YSAT(3)*YSAT(5)
      CRSVEC(2) = YSAT(3)*YSAT(4)-YSAT(1)*YSAT(6)
      CRSVEC(3) = YSAT(1)*YSAT(5)-YSAT(2)*YSAT(4)
      ABSCRS = DSQRT(CRSVEC(1)**2+CRSVEC(2)**2+CRSVEC(3)**2)
      DO K=1,3
         CRSVEC(K) = CRSVEC(K)/ABSCRS
      ENDDO

      ALGVEC(1) = CRSVEC(2)*RADVEC(3)-CRSVEC(3)*RADVEC(2)
      ALGVEC(2) = CRSVEC(3)*RADVEC(1)-CRSVEC(1)*RADVEC(3)
      ALGVEC(3) = CRSVEC(1)*RADVEC(2)-CRSVEC(2)*RADVEC(1) 

      IF(ERM.GT.0)THEN

C        DISTANCE FROM SATELLITE TO SUN
         RSUN = DSQRT((YSAT(1)-SUN(1))**2+(YSAT(2)-SUN(2))**2+
     1           (YSAT(3)-SUN(3))**2)

C        D VECTOR AND Z VECTOR
         DO K=1,3
            D_SUN(K) = (SUN(K)-YSAT(K))/RSUN
            Z_SAT(K) = -YSAT(K)/ABSPOS
         ENDDO

C        Y VECTOR
         Y0SAT(1) = Z_SAT(2)*D_SUN(3)-Z_SAT(3)*D_SUN(2)
         Y0SAT(2) = Z_SAT(3)*D_SUN(1)-Z_SAT(1)*D_SUN(3)
         Y0SAT(3) = Z_SAT(1)*D_SUN(2)-Z_SAT(2)*D_SUN(1)
         ABSY0V = DSQRT(Y0SAT(1)**2 + Y0SAT(2)**2 + Y0SAT(3)**2)
         DO K=1,3
            Y_SAT(K) = Y0SAT(K)/ABSY0V
         ENDDO

C        B VECTOR
         B_SUN(1) = Y_SAT(2)*D_SUN(3) - Y_SAT(3)*D_SUN(2)
         B_SUN(2) = Y_SAT(3)*D_SUN(1) - Y_SAT(1)*D_SUN(3)
         B_SUN(3) = Y_SAT(1)*D_SUN(2) - Y_SAT(2)*D_SUN(1)

C        X VECTOR
         X_SAT(1) = Y_SAT(2)*Z_SAT(3) - Y_SAT(3)*Z_SAT(2)
         X_SAT(2) = Y_SAT(3)*Z_SAT(1) - Y_SAT(1)*Z_SAT(3)
         X_SAT(3) = Y_SAT(1)*Z_SAT(2) - Y_SAT(2)*Z_SAT(1)

         DO K=1,3
            ATTSURF(K,1) = Z_SAT(K)
            ATTSURF(K,2) = Y_SAT(K)
            ATTSURF(K,3) = X_SAT(K)
            ATTSURF(K,4) = D_SUN(K)
         ENDDO
      ENDIF

C ---------------------------- 
C OPTICAL PROPERTIES PER BLOCK
C ----------------------------

      IF(ERM.GT.0)THEN
         IF(BLKNUM.LE.10)THEN
            INDB = BLKNUM
         ELSEIF(BLKNUM.GT.100)THEN
            INDB = BLKNUM - 90
         ENDIF

         DO II = 1,4
            DO JJ = 1,2
               AREAS(II,JJ) = AREA(II,JJ,INDB)
               REFLS(II,JJ) = REFL(II,JJ,INDB)
               DIFUS(II,JJ) = DIFU(II,JJ,INDB)
               ABSPS(II,JJ) = ABSP(II,JJ,INDB)

               AREA2S(II,JJ) = AREA2(II,JJ,INDB)
               REFL2S(II,JJ) = REFL2(II,JJ,INDB)
               DIFU2S(II,JJ) = DIFU2(II,JJ,INDB)
               ABSP2S(II,JJ) = ABSP2(II,JJ,INDB)

               REFLIRS(II,JJ) = REFLIR(II,JJ,INDB)
               DIFUIRS(II,JJ) = DIFUIR(II,JJ,INDB)
               ABSPIRS(II,JJ) = ABSPIR(II,JJ,INDB)
            ENDDO
         ENDDO
      ENDIF

C ----------------------
C EARTH RADIATION MODELS
C ----------------------

      IF(ERM.GT.0)THEN
         ABSSUN = DSQRT(SUN(1)**2 + SUN(2)**2 + SUN(3)**2)
         DO K=1,3
            ESUN(K) = SUN(K)/ABSSUN
         ENDDO

         PSIDOT = ESUN(1)*RADVEC(1)+ESUN(2)*RADVEC(2)+ESUN(3)*RADVEC(3)
         IF(DABS(PSIDOT).GT.(1D0-1D-6))THEN
            PSI = 0D0
         ELSE
            PSI = DACOS(PSIDOT)
         ENDIF
         S0 = S0*(AU/ABSSUN)**2
      ENDIF

C     ANALYTICAL MODEL
      IF(ERM.EQ.1)THEN

         NCFVEC(1) = RADVEC(1)
         NCFVEC(2) = RADVEC(2)
         NCFVEC(3) = RADVEC(3)
         ALBFAC = (PI*TOA**2)*(S0/C)/(ABSPOS**2)
         PHASEVI = (2*ALB/(3*PI**2))*((PI-PSI)*DCOS(PSI)+DSIN(PSI))
         PHASEIR = (1-ALB)/(4*PI)
         ABSNCFVI = ALBFAC*PHASEVI
         ABSNCFIR = ALBFAC*PHASEIR

         CALL SURFBOXW(AREAS,REFLS,DIFUS,ABSPS,
     1                    AREA2S,REFL2S,DIFU2S,ABSP2S,
     2                    REFLIRS,DIFUIRS,ABSPIRS,
     3                    ABSNCFVI,ABSNCFIR,NCFVEC,ATTSURF,FORCE)



C     NUMERICAL MODEL (CERES DATA)
      ELSEIF(ERM.EQ.2)THEN

         ABSSUN = DSQRT(SUN(1)**2 + SUN(2)**2 + SUN(3)**2)
         DO K=1,3
            ESUN(K) = SUN(K)/ABSSUN
         ENDDO

         DO LATK = 1,LATKMX
            DO LONK = 1,LONKMX

               D_AREA = D_AREA_ALL(LATK,LONK)
               V_NS(1) = V_NS_ALL(LATK,LONK,1)
               V_NS(2) = V_NS_ALL(LATK,LONK,2)
               V_NS(3) = V_NS_ALL(LATK,LONK,3)

               DO II=1,3
                  V_INS(II)=0D0
                  DO JJ=1,3
                     V_INS(II) = V_INS(II) + KAPPA(JJ,II)*V_NS(JJ)
                  ENDDO
               ENDDO

C              Distance and direction from point in the Earth to satellite
               V_DIST(1) = YSAT(1)-TOA*V_INS(1) 
               V_DIST(2) = YSAT(2)-TOA*V_INS(2)
               V_DIST(3) = YSAT(3)-TOA*V_INS(3)
               DIST2 = V_DIST(1)**2 +V_DIST(2)**2 +V_DIST(3)**2
               ABSDIST = DSQRT(DIST2)
               V_SAT(1) = V_DIST(1)/ABSDIST
               V_SAT(2) = V_DIST(2)/ABSDIST
               V_SAT(3) = V_DIST(3)/ABSDIST

C              Cosine of angles of incident and reflected radiation
               COS_IN = ESUN(1)*V_INS(1) + ESUN(2)*V_INS(2) 
     1                + ESUN(3)*V_INS(3)
               COS_RE =V_SAT(1)*V_INS(1) +V_SAT(2)*V_INS(2)
     1                +V_SAT(3)*V_INS(3)

               IF(COS_RE.GE.0)THEN

C              Reflectivity and emissivity coefficients
                  REFL_CF = CERGRE(LATK,LONK)
                  EMIT_CF = CERGEM(LATK,LONK)

C                 Reflected Irradiance
                  IF(COS_IN.GE.0)THEN
                     E_REFL=(REFL_CF/(PI*DIST2))*COS_RE*COS_IN*S0*D_AREA
                  ELSE
                     E_REFL=0D0
                  ENDIF

C                 Emitted Irradiance
                  E_EMIT = (EMIT_CF/(4*PI*DIST2))*COS_RE*S0*D_AREA

C                 Non-conservative force
                  ABSNCFVI = E_REFL/C
                  ABSNCFIR = E_EMIT/C
                  NCFVEC(1) = V_SAT(1)
                  NCFVEC(2) = V_SAT(2)
                  NCFVEC(3) = V_SAT(3)

                  CALL SURFBOXW(AREAS,REFLS,DIFUS,ABSPS,
     1                    AREA2S,REFL2S,DIFU2S,ABSP2S,
     2                    REFLIRS,DIFUIRS,ABSPIRS,
     3                    ABSNCFVI,ABSNCFIR,NCFVEC,ATTSURF,FORCE_LL)

                  FORCE(1) = FORCE(1) + FORCE_LL(1)
                  FORCE(2) = FORCE(2) + FORCE_LL(2)
                  FORCE(3) = FORCE(3) + FORCE_LL(3)
               ENDIF

            ENDDO
         ENDDO

      ENDIF


C     ANTENNA POWER OF GPS SATELLITES (IN WATTS)
C     IGS MODEL (JIM RAY, 2011)

C     GPS BLOCK IIA (ASSUMED THE SAME FOR BLOCK I AND II) 
      IF(BLKNUM.LE.3)THEN
         ANTPOW = 76D0

C     GPS BLOCK IIR
      ELSEIF((BLKNUM.GE.4).AND.(BLKNUM.LE.6))THEN
         ANTPOW = 85D0

C     GPS BLOCK IIR-M
      ELSEIF(BLKNUM.EQ.7)THEN
         ANTPOW = 198D0

C     GPS BLOCK IIF
      ELSEIF(BLKNUM.EQ.8)THEN
         ANTPOW = 249D0
         IF((SVN.EQ.62).AND.(MJD.GE.55656D0))THEN
C           NO M-CODE FOR SVN62/PRN25 STARTING 05APR2011; NANU 2011026
            ANTPOW = 154D0
         ENDIF
      ENDIF

C     NO ANTENNA POWER INFORMATION FOR GLONASS SATELLITES


C     NAVIGATION ANTENNA THRUST (SIMPLE MODEL)
      IF((ANT.EQ.1).AND.(BLKNUM.LT.100))THEN
         ANTFORCE = ANTPOW/C
         FORCE(1) = FORCE(1) + ANTFORCE*RADVEC(1)
         FORCE(2) = FORCE(2) + ANTFORCE*RADVEC(2)
         FORCE(3) = FORCE(3) + ANTFORCE*RADVEC(3)
      ENDIF

      IF((REFF.GT.0).AND.((ERM.GT.0).OR.(ANT.EQ.1)))THEN
         DO K=1,3
            FREF(K) = 0D0
         ENDDO

C        FORCE IN BODY-FIXED REFERENCE FRAME
         IF(REFF.EQ.1)THEN
            DO K=1,3
               FREF(1) = FREF(1) + FORCE(K)*Z_SAT(K)
               FREF(2) = FREF(2) + FORCE(K)*Y_SAT(K)
               FREF(3) = FREF(3) + FORCE(K)*X_SAT(K)
            ENDDO

C        FORCE IN SUN-FIXED REFERENCE FRAME
         ELSEIF(REFF.EQ.2)THEN
            DO K=1,3
               FREF(1) = FREF(1) + FORCE(K)*D_SUN(K)
               FREF(2) = FREF(2) + FORCE(K)*Y_SAT(K)
               FREF(3) = FREF(3) + FORCE(K)*B_SUN(K)
            ENDDO

C        FORCE IN ORBITAL REFERENCE FRAME
         ELSEIF(REFF.EQ.3)THEN
            DO K=1,3
               FREF(1) = FREF(1) + FORCE(K)*RADVEC(K)
               FREF(2) = FREF(2) + FORCE(K)*ALGVEC(K)
               FREF(3) = FREF(3) + FORCE(K)*CRSVEC(K)
            ENDDO
         ENDIF

         DO K=1,3
            FORCE(K) = FREF(K)
         ENDDO
      ENDIF


C     MASS OF SATELLITES
      IF(BLKNUM.EQ.1)THEN
         MASS = 500D0
      ELSEIF(BLKNUM.EQ.2)THEN 
         MASS = 880D0
      ELSEIF(BLKNUM.EQ.3)THEN
         MASS = 975D0
      ELSEIF((BLKNUM.GE.4).AND.(BLKNUM.LE.7))THEN
         MASS = 1100D0
      ELSEIF(BLKNUM.EQ.8)THEN
         MASS = 1555D0
      ELSEIF((BLKNUM.EQ.101).OR.(BLKNUM.EQ.102))THEN
         MASS = 1415D0
      ELSEIF(BLKNUM.EQ.103)THEN
         MASS = 935D0
      ENDIF

C     CONVERSION TO ACCELERATION
      DO K=1,3
         ACCEL(K) = FORCE(K)/MASS
      ENDDO

      END SUBROUTINE
