      DOUBLE PRECISION   TRISE,      TCA,        TSET,                          
     $                   TNODE,      TSTART,     TSTOP,                         
     $                   TNEXT,      EQXYDA,     DMOD                           
      DIMENSION          ISAT(10),   PER(10),    WMTN(10),                      
     $                   EQXYDA(10), EQXMIN(10), EQXLON(10),                    
     $                   TRISE(10),  TCA(10),    TSET(10),                      
     $                   TNODE(10),  IWEEK(2,7), MONTH(2,12),                   
     $                   IECA(10),   ICADIR(10), IDIR(10)                       
      DIMENSION          ITCLOS(5),  ITSET(5),   IP(5),                         
     $                   ISIDE(10),  ITIME(5),   ITRISE(5),                     
     $                   JYEAR(10)                                              
      EQUIVALENCE        (IP(1),LUCON),                                         
     $                   (IP(2),LULIST),                                        
     $                   (IP(3),IONLIN),                                        
     $                   (IP(4),NSAT)                                           
      DATA                                                                      
     $     IWEEK   /2HMO,2HN ,2HTU,2HE ,2HWE,2HD ,2HTH,2HU ,                    
     $              2HFR,2HI ,2HSA,2HT ,2HSU,2HN /,                             
     $     MONTH   /2HJA,2HN ,2HFE,2HB ,2HMA,2HR ,2HAP,2HR ,                    
     $              2HMA,2HY ,2HJU,2HN ,2HJU,2HL ,2HAU,2HG ,                    
     $              2HSE,2HP ,2HOC,2HT ,2HNO,2HV ,2HDE,2HC /                    
C$                                                                              
C NAME        ALERT                                                             
C                                                                               
C TYPE        PROGRAM                                                           
C                                                                               
C PURPOSE     COMPUTE TRANSIT SATELLITE ALERTS FOR UP TO 10 SATELLITES          
C                                                                               
C AUTHOR/DATE D.WELLS 800130   (DEFAULT DATA 800124)                            
C                                                                               
C EXTERNALS   RMPAR,EXEC,FMOD,FLOAT,UPDAT,NEXT,MOD                              
C                                                                               
C CALLING     RU,ALERT,LUCON,LULIST,IONLIN,NSAT                                 
C                                                                               
C PARAMETERS  LUCON  = CONSOLE LU                                               
C             LULIST = LISTING LU (DEFAULT = LUCON)                             
C             IONLIN = 0 FOR OFFLINE ALERTS (DEFAULT)                           
C                      1 FOR ONLINE ALERTS                                      
C             NSAT   = NUMBER OF SATELLITES (DEFAULT = 6)                       
C                                                                               
C LANGUAGE    FORTRAN                                                           
C                                                                               
C COMMENTS    PROGRAM FIRST CHECKS TO SEE IF FILE "NODAL"                       
C             EXISTS, AND IF SO, WHETHER IT CONTAINS UP-TO-DATE                 
C             NODAL CROSSING DATA (PLACED THERE BY PROGRAM SATFX).              
C             IF NOT, DEFAULT DATA FROM DATA STATEMENT IS USED.                 
C                                                                               
C             ONLINE MODE COMPUTES NEXT 10 ALERTS FOR CURRENT POSITION          
C                                                                               
C             OFFLINE MODE LETS OPERATOR SELECT POSITION, START AND             
C             STOP TIMES.                                                       
C                                                                               
C             PROGRAM "ACTIV" SCHEDULES ALERT IN ONLINE MODE ONCE               
C             PER HOUR.                                                         
C                                                                               
C                                                                               
C**********PROGRAM ALERT INPUT**********                                        
C                                                                               
C        COLUMNS   FORMAT    NAME      DESCRIPTION                              
C        _________________________________________                              
C                                                                               
C  CARD 1                                                                       
C                                                                               
C        1-4       I4        LUCON     INPUT LOGICAL UNIT NO.(DEFAULT=5)        
C        5-8       I4        LULIST    OUTPUT LOGICAL UNIT NO.(DEFAULT=6)       
C        9-12      I4        IONLIN    =0 FOR OFFLINE ALERTS                    
C                                      =1 FOR ONLINE ALERTS                     
C        13-16     I4        NSAT      NUMBER OF SATELLITES                     
C                                                                               
C                                                                               
C  CARD 2                                                                       
C                                                                               
C       1-10      I10       IYEAR     YEAR FOR WHICH ALERTS REQUIRED            
C       11-20     F10.4     TSTAX     FIRST DAY FOR WHICH ALERTS REQUIRED       
C       21-30     F10.4     TENDX     LAST DAY FOR WHICH ALERTS REQUIRED        
C                                     (DEFAULT=TSTAX)                           
C       31-40     F10.4     OBSLAX    LATITUDE OF OBSERVING STATION             
C       41-50     F10.4     OBSLOX    LONGITUDE OF OBSERVING STATION            
C                                                                               
C                                                                               
C  CARD 3 +     NOTE:  ONE CARD OF THIS TYPE IS NEEDED FOR EACH SATELLITE       
C                      FOR WHICH ALERTS ARE TO BE COMPUTED                      
C                                                                               
C        1-10      F10.4     SATX      SATELLITE NUMBER                         
C        11-20     F10.4     PERX      SATELLITE ORBITAL PERIOD(MINUTES)        
C        21-30     F10.4     WMNTX     SATELLITE WESTWARD MOTION(DEGREES)       
C        31-40     F10.4     EQXDAX    DATE OF OBSERVATION (YEAR . DAY)         
C        41-50     F10.4     EQXMIX    TIME OF EQUATOR CROSSING(MINUTES)        
C        51-60     F10.4     EQXLOX    LONGITUDE OF EQUATOR CROSSING            
C                                      (DEGREES (EAST))                         
C                                                                               
C                                                                               
C***********************************************************************        
C                                                                               
C$                                                                              
C GET RUN PARAMETERS                                                            
      READ(5,*)(IP(I),I=1,4)                                                    
      WRITE(LULIST,9000)                                                        
 9000 FORMAT(1H1,/,/,10X,'ALERT INPUT',                                         
     &             /,10X,'===========')                                         
      WRITE(LULIST,9001)(IP(I),I=1,4)                                           
 9001  FORMAT(3X,'RUN PARAMETERS : ',5I4)                                       
C1005 FORMAT(5I4)                                                               
      IF(LUCON.EQ.0)LUCON=5                                                     
      IF(LULIST.EQ.0)LULIST=6                                                   
C READ FILE NODAL FOR INPUT                                                     
      CALL NODAL(LUCON,LULIST,IONLIN,NSAT,                                      
     $           ISAT,PER,WMTN,EQXYDA,EQXMIN,EQXLON,                            
     $           OBSLAT,OBSLON,IYEAR,TSTART,TSTOP)                              
C FOR OFFLINE MODE CONFIRM INPUT WITH OPERATOR                                  
      IF(IONLIN .NE. 0) GO TO 5                                                 
      CALL OPERA(LUCON,LULIST,NSAT,                                             
     $           ISAT,PER,WMTN,EQXYDA,EQXMIN,EQXLON,                            
     $           OBSLAT,OBSLON,IYEAR,TSTART,TSTOP)                              
    5 CONTINUE                                                                  
      OBSLON = FMOD(OBSLON,360.,N)                                              
      NALERT = 0                                                                
      LBLDAY = TSTART - 1D0                                                     
C COMPUTE TNODE IN DAYS                                                         
      DO 10 I = 1,NSAT                                                          
        JYEAR(I) = EQXYDA(I)                                                    
        EQXDAY = DMOD(EQXYDA(I),1D0) * 1000D0                                   
        CALL XYEAR(IYEAR,JYEAR(I),EQXDAY)                                       
   10   TNODE(I) = EQXDAY + EQXMIN(I) / 1440D0                                  
C FILL ALERT BUFFER                                                             
      DO 20 I = 1,NSAT                                                          
        NORB = (TSTART - TNODE(I)) * 1440D0 / PER(I) - 1                        
        CALL UPDAT(LULIST,PER(I),WMTN(I),TNODE(I),EQXLON(I),NORB)               
   15   CALL NEXT(LULIST,OBSLAT,OBSLON,                                         
     $            PER(I),WMTN(I),TNODE(I),EQXLON(I),                            
     $            TRISE(I),TCA(I),TSET(I),                                      
     $            IECA(I),ICADIR(I),IDIR(I),ISIDE(I))                           
        IF(TCA(I) .LT. TSTART) GO TO 15                                         
   20   CONTINUE                                                                
C FIND EARLIEST ALERT IN BUFFER                                                 
   30 TNEXT = TCA(1)                                                            
      J = 1                                                                     
      DO 40 I = 2,NSAT                                                          
        IF(TCA(I) .GT. TNEXT) GO TO 40                                          
        TNEXT = TCA(I)                                                          
        J = I                                                                   
   40   CONTINUE                                                                
C QUIT AFTER TSTOP                                                              
      IF(TCA(J) .GT. TSTOP) GO TO 99                                            
C LABEL EACH DAYS ALERTS                                                        
      IF(IDINT(TCA(J)) .LE. LBLDAY) GO TO 45                                    
      LBLDAY = LBLDAY + 1                                                       
      ND = LBLDAY                                                               
      IY = IYEAR                                                                
      CALL DATE(IY,ND,MNUM,MDAY,IW)                                             
      OUTLON = OBSLON                                                           
      IF(OUTLON .GT. 180.) OUTLON = OUTLON - 360.                               
      WRITE(LULIST,1001) IY,ND,                                                 
     $       IWEEK(1,IW),IWEEK(2,IW),MDAY,MONTH(1,MNUM),MONTH(2,MNUM),          
     $                   OBSLAT,OUTLON                                          
C COMPUTE RISE/CA/SET   HR/MIN                                                  
   45 CALL TIME(TRISE(J),ITRISE)                                                
      CALL TIME(TCA(J),ITCLOS)                                                  
      CALL TIME(TSET(J),ITSET)                                                  
C WRITE LATEST ALERT                                                            
      WRITE(LULIST,1002) ISAT(J),IECA(J),(ITRISE(7-I),I=3,4),                   
     $                   (ITCLOS(7-I),I=3,4),(ITSET(7-I),I=3,4),                
     $                   ISIDE(J),IDIR(J),ICADIR(J)                             
C REFILL EMPTY BUFFER SLOT                                                      
      CALL NEXT(LULIST,OBSLAT,OBSLON,                                           
     $          PER(J),WMTN(J),TNODE(J),EQXLON(J),                              
     $          TRISE(J),TCA(J),TSET(J),                                        
     $          IECA(J),ICADIR(J),IDIR(J),ISIDE(J))                             
C CHECK FOR MAXIMUM OF 10 ALERTS ONLINE                                         
      NALERT = NALERT + 1                                                       
      IF(IONLIN .NE. 0 .AND. NALERT .GE. 10) GO TO 999                          
      GO TO 30                                                                  
C ALERT DONE MESSAGE                                                            
   99 WRITE(LULIST,1003)                                                        
  999 CONTINUE                                                                  
 1001 FORMAT(1H1,/,/,'    ALERTS FOR ',I5,' DAY ', I4,                          
     $          '    (',2A2,I3,1X,2A2,')  ',                                    
     $          ' LAT',F6.1,' LON',F6.1,/,/,                                    
     $  '    SAT  ELEV       RISE         CA       SET      SAT DIR')           
 1002 FORMAT(1X,2I6,3(I7,I4),5X,2A2,I7)                                         
 1003 FORMAT(/,/,/,'  ALERTS DONE|')                                            
      STOP                                                                      
      END                                                                       
      SUBROUTINE NODAL(LUCON,LULIST,IONLIN,NSAT,                                
     $                 ISAT,PER,WMTN,EQXYDA,EQXMIN,EQXLON,                      
     $                 OBSLAT,OBSLON,IYEAR,TSTART,TSTOP)                        
      DOUBLE PRECISION TSTART,TSTOP,DAY,TSTAX,EQXYDA,EQXDAX                     
      DIMENSION        ISAT(10),   PER(10),    WMTN(10),                        
     $                 EQXYDA(10), EQXMIN(10), EQXLON(10),                      
     $                 NPROG(3),   NFILE(3),   IDCB(144),                       
     $                 IBUF(30),   ITIME(6),   IBLANK(30)                       
      DATA    NPROG    /2HAL,2HER,2HT /,                                        
     $        NFILE    /2HNO,2HDA,2HL /,                                        
     $        IBLANK   /30*2H  /                                                
C$                                                                              
C NAME        NODAL                                                             
C                                                                               
C TYPE        SUBROUTINE                                                        
C                                                                               
C PURPOSE     READ FILE "NODAL" IF IT EXISTS                                    
C                                                                               
C CALLING     CALL NODAL(LUCON,LULIST,IONLIN,NSAT,                              
C                        ISAT,PER,WMTN,EQXYDA,EQXMIN,EQXLON,                    
C                        OBSLAT,OBSLON,TSTART,TSTOP)                            
C$                                                                              
C CHECK WHETHER FILE NODAL EXISTS                                               
      WRITE(LULIST,9100)                                                        
 9100 FORMAT(/,3X,'YEAR       TSAT      TENDX     OBSLAX     OBSLOX')           
      READ(5,*) IYEAR,TSTAX,TENDX,OBSLAX,OBSLOX                                 
      WRITE(LULIST,9300) IYEAR,TSTAX,TENDX,OBSLAX,OBSLOX                        
 9300 FORMAT(3X,I4,1X,F10.4,1X,F10.4,1X,F10.4,1X,F10.4,1X,F10.4)                
      CALL TIME(TSTAX,ITIME)                                                    
      OBSLAT = OBSLAX                                                           
      OBSLON = OBSLOX                                                           
C READ REST OF FILE (LATEST NODAL PARAMETERS)                                   
      WRITE(LULIST,9350)                                                        
 9350 FORMAT(/,1X,'  SATX      PERX       WMTNS',                               
     &'      EQXDAX      EQXMIN      EQXLON')                                   
      DO 30  I=1,NSAT                                                           
      READ(5,*) SATX,PERX,WMTNX,EQXDAX,EQXMIX,EQXLOX                            
      WRITE(LULIST,9400)SATX,PERX,WMTNX,EQXDAX,EQXMIX,EQXLOX                    
 9400 FORMAT(1X,F6.1,6(F10.4,2X))                                               
      J=I                                                                       
      ISAT(J)=SATX                                                              
        PER(J)  = PERX                                                          
        WMTN(J) = WMTNX                                                         
        EQXYDA(J) = EQXDAX                                                      
        EQXMIN(J) = EQXMIX                                                      
        EQXLON(J) = EQXLOX                                                      
   30   CONTINUE                                                                
C OFFLINE- COMPUTE 5 DAYS  ALERTS FROM LAST FIX TIME                            
      ITIME(1) = 0                                                              
      ITIME(2) = 0                                                              
      ITIME(3) = 0                                                              
      TSTART = DAY(ITIME)                                                       
      TSTOP = TSTART + 7D0 / 24D0                                               
      IF(IONLIN.EQ.0) TSTOP = TENDX + 1.D0                                      
      IF(TENDX.EQ.0.) TSTOP = TSTART + 1.D0                                     
      RETURN                                                                    
 1001 FORMAT(I10,4F10.4)                                                        
 1002 FORMAT(6F10.4)                                                            
      END                                                                       
      SUBROUTINE OPERA(LUCON,LULIST,NSAT,                                       
     $                 ISAT,PER,WMTN,EQXYDA,EQXMIN,EQXLON,                      
     $                 OBSLAT,OBSLON,IYEAR,TSTART,TSTOP)                        
      DOUBLE PRECISION TSTART,  TSTOP,  EQXYDA                                  
      DIMENSION        ISAT(10),PER(10), WMTN(10),                              
     $                 EQXYDA(10),EQXMIN(10),EQXLON(10)                         
C$                                                                              
C NAME        OPERA                                                             
C                                                                               
C TYPE        SUBROUTINE                                                        
C                                                                               
C PURPOSE     INPUT-CHECKING DIALOGUE WITH OPERATOR                             
C                                                                               
C CALLING     CALL OPERA(LUCON,LULIST,NSAT,                                     
C                        OBSLAT,OBSLON,TSTART,TSTOP)                            
C$                                                                              
C ECHO INPUT ON LULIST                                                          
      IDAY = TSTART                                                             
      IHR = (TSTART - IDAY) * 24D0                                              
      JDAY = TSTOP                                                              
      JHR = (TSTOP - JDAY) * 24D0                                               
      WRITE(LULIST,1006)OBSLAT,OBSLON,IYEAR,IDAY,IHR,JDAY,JHR,                  
     $                   (ISAT(I),PER(I),WMTN(I),EQXYDA(I),                     
     $                    EQXMIN(I),EQXLON(I),I=1,NSAT)                         
C FORMAT STATEMENTS                                                             
 1006 FORMAT(1H1,/,'   ALERT NORTH LATITUDE =    ',F10.4,/,                     
     $       '         EAST LONGITUDE =    ',F10.4,/,                           
     $       '         START YR/DA/HR =    ',3I4,/,                             
     $       '         STOP  DAY/HOUR =    ',2I4,//,                            
     $      '   SAT   PERIOD    WMOTN   YEAR.DAY     EQXMIN    EQXLON',/,       
     $       10(3X,I3,5F10.4/))                                                 
      RETURN                                                                    
      END                                                                       
      SUBROUTINE NEXT(LU,OBSLAT,OBSLON,                                         
     $                PER,WMTN,TNODE,EQXLON,                                    
     $                TRISE,TCA,TSET,                                           
     $                IECA,ICADIR,IDIR,ISIDE)                                   
      DOUBLE PRECISION TNODE,TRISE,TCA,TSET                                     
      DIMENSION INS(2),IEW(2)                                                   
      DATA PI/3.14159/,                                                         
     $     REARTH /6365./,                                                      
     $     GRAV   /398600.5/,                                                   
     $     INS    /2HNS,2HSN/,                                                  
     $     IEW    /2HW ,2HE /                                                   
C$                                                                              
C NAME        NEXT                                                              
C                                                                               
C TYPE        SUBROUTINE                                                        
C                                                                               
C PURPOSE     COMPUTE ALERT PARAMETERS FOR NEXT PASS OF ONE SATELLITE           
C                                                                               
C AUTHOR/DATE D.WELLS 790517                                                    
C                                                                               
C EXTERNALS   COS,SIN,ATAN,ATAN2                                                
C                                                                               
C CALLING     CALL NEXT(LU,OBSLAT,OBSLON,                                       
C                       PER,WMTN,TNODE,EQXLON,                                  
C                       TRISE,TCA,TSET,IECA,ICADIR,IDIR,ISIDE)                  
C                                                                               
C PARAMETERS  LU = LISTING LU                                                   
C            OBSLAT = OBSERVER'S LATITUDE (DEGREES NORTH)                       
C            OBSLON = OBSERVER'S LONGITUDE (DEGREES EAST)                       
C            PER = SATELLITE ORBITAL PERIOD (MINUTES)                           
C            WMTN = SATELLITE WESTWARD MOTION (DEGREES)                         
C            TNODE = NODAL CROSSING TIME (DAYS)                                 
C            EQXLON = NODAL CROSSING LONGITUDE (DEGREES EAST)                   
C            TRISE = RISE TIME (DAYS)                                           
C            TCA   = CLOSEST APPROACH TIME (DAYS)                               
C            TSET  = SET TIME (DAYS)                                            
C            IECA = ELEVATION AT TCA (DEGREES)                                  
C            ICADIR = SATELLITE DIRECTION AT TCA (DEGREES)                      
C            IDIR = SATELLITE DIRECTION (1 = NS, 2 = SN)                        
C            ISIDE = SATELLITE SIDE (1 = WEST, 2 = EAST)                        
C                                                                               
C LANGUAGE    FORTRAN                                                           
C$                                                                              
C USE KEPLER'S THIRD LAW  TO GET ORBIT RADIUS FROM PERIOD                       
      RORBIT = ((PER * 60. / (2. * PI))**2 * GRAV) ** (1./3.)                   
      COSAMX = REARTH / RORBIT                                                  
      SOBLAT = SIN(OBSLAT*PI/180.)                                              
      COBLAT = COS(OBSLAT*PI/180.)                                              
      TOBLAT = SOBLAT / COBLAT                                                  
C COMPUTE SATELLITE POSITION AT TCA (SUBLAT,SUBLON)                             
   10 IVIS = 0                                                                  
      DLONG = FMOD(EQXLON-OBSLON,360.,N)                                        
      SDLON = SIN(DLONG*PI/180.)                                                
      CDLON = COS(DLONG*PI/180.)                                                
      TDLON = SDLON / CDLON                                                     
      SUBLAT = OBSLAT                                                           
      IF(ABS(DLONG) .GE. 90. .AND. ABS(DLONG) .LE. 270.)                        
     $         SUBLAT = 180. - OBSLAT                                           
      DO 20 ITER = 1,10                                                         
        SUB2 = SUBLAT                                                           
        SUBLON = FMOD(DLONG - WMTN * SUBLAT / 360., 360.,N)                     
        SLON = SIN(SUBLON*PI/180.)                                              
        CLON = COS(SUBLON*PI/180.)                                              
        SUBLAT = ATAN2(TOBLAT + WMTN * SLON / 360.,CLON)*180. / PI              
        IF(ABS(SUBLAT - SUB2) .LT. 0.1) GO TO 25                                
   20   CONTINUE                                                                
      GO TO 40                                                                  
C IS PASS VISIBLE AT TCA                                                        
   25 SLAT = SIN(SUBLAT*PI/180.)                                                
      CLAT = COS(SUBLAT*PI/180.)                                                
      CALF = CLAT * CLON * COBLAT + SLAT * SOBLAT                               
      IF(CALF .LT. COSAMX) GO TO 40                                             
C VISIBLE PASS FOUND                                                            
      IVIS = 1                                                                  
      TCA = TNODE + (PER / 1440D0) * (SUBLAT / 360D0)                           
      DURAT = (PER / 1440D0) * (ARCCOS(COSAMX / CALF) / 360.)                   
      TRISE = TCA - DURAT                                                       
      TSET  = TCA + DURAT                                                       
      SALF = SQRT(1. - CALF **2)                                                
      IECA=90                                                                   
      AZCA= 0.                                                                  
      IF( ABS(SALF) .LT. 1.E-5) GO TO 35                                        
      IECA = ATAN((CALF - COSAMX) / SALF) * 180. / PI                           
      AZCA = ARCCOS((SLAT - SOBLAT * CALF) / (COBLAT * SALF))                   
   35 AZCA = (AZCA+ 45.*(1.-SIGN(1.,OBSLAT)))*SIGN(1.,SDLON/CDLON)              
      ICADIR = AZCA - 90. * SIGN(1.,SDLON)                                      
      IF(ICADIR .LT. 0) ICADIR = ICADIR + 360                                   
      IDIR = (IFIX(CDLON/ABS(CDLON)) + 3) / 2                                   
      ISIDE = (IFIX(TDLON/ABS(TDLON)) + 3) / 2                                  
      IDIR = INS(IDIR)                                                          
      ISIDE = IEW(ISIDE)                                                        
C UPDATE TNODE,EQXLON                                                           
   40 CALL UPDAT(LU,PER,WMTN,TNODE,EQXLON,1)                                    
C IF NO ALERT TRY NEXT PASS                                                     
      IF(IVIS .EQ. 0) GO TO 10                                                  
      RETURN                                                                    
      END                                                                       
      SUBROUTINE UPDAT(LU,PER,WMTN,TNODE,EQXLON,NORB)                           
      DOUBLE PRECISION TNODE                                                    
C$                                                                              
C NAME        UPDATE                                                            
C                                                                               
C TYPE        SUBROUTINE                                                        
C                                                                               
C PURPOSE     UPDATE TNODE,EQXLON BY NORB ORBITS                                
C                                                                               
C AUTHOR/DATE D.WELLS 790517                                                    
C                                                                               
C CALLING     CALL UPDAT(LU,PER,WMTN,TNODE,EQXLON,NORB)                         
C                                                                               
C LANGUAGE    FORTRAN                                                           
C$                                                                              
      TNODE = TNODE + NORB * PER / 1440D0                                       
      EQXLON = FMOD(EQXLON - NORB * WMTN,360.,N)                                
      RETURN                                                                    
      END                                                                       
      FUNCTION ARCCOS(COSA)                                                     
      DATA PI/3.14159/                                                          
C$                                                                              
C NAME        ARCCOS                                                            
C                                                                               
C TYPE        FUNCTION                                                          
C                                                                               
C PURPOSE     RETURNS ARC COSINE IN DEGREES                                     
C                                                                               
C CALLING     A = ARCCOS(COSA)                                                  
C$                                                                              
      ARCCOS = 90.                                                              
      IF(COSA .EQ. 0.) RETURN                                                   
      ARCCOS = 0.                                                               
      IF(ABS(COSA) .GE. 1.) RETURN                                              
      ARCCOS = ATAN(SQRT(1. - COSA**2) / COSA) * 180. / PI                      
      RETURN                                                                    
      END                                                                       
      FUNCTION FMOD(VALUE,CYCLE,N)                                              
C$                                                                              
C NAME        FMOD                                                              
C                                                                               
C TYPE        FUNCTION                                                          
C                                                                               
C PURPOSE     GENERALIZED MOD FUNCTION                                          
C                                                                               
C CALLING     X = FMOD(VALUE,CYCLE,N)                                           
C                                                                               
C PARAMETERS  VALUE = INPUT VALUE                                               
C             CYCLE = INPUT MODULUS RANGE                                       
C             N = OUTPUT NUMBER OF CYCLES TO NORMALIZE VALUE                    
C             FMOD = OUTPUT NORMALIZED VALUE,                                   
C                    IN THE RANGE (0,CYCLE) IF CYCLE POSITIVE                   
C                    IN THE RANGE (CYCLE,0) IF CYCLE NEGATIVE                   
C$                                                                              
      N = VALUE / CYCLE                                                         
      FMOD = VALUE - N * CYCLE                                                  
      IF(FMOD * CYCLE .GE. 0.) RETURN                                           
      FMOD = FMOD + CYCLE                                                       
      N = N - 1                                                                 
      RETURN                                                                    
      END                                                                       
      FUNCTION JERR(LUCON,NPROG,NFILE,ID,IER)                                   
      DIMENSION NPROG(3),NFILE(3)                                               
C$                                                                              
C NAME        JERR                                                              
C                                                                               
C TYPE        FUNCTION                                                          
C                                                                               
C PURPOSE     PRINT FMP ERROR MESSAGES                                          
C                                                                               
C CALLING     I = JERR(LUCON,NPROG,NFILE,ID,IER)                                
C                                                                               
C PARAMETERS  LUCON = CONSOLE LU                                                
C             NPROG(3) = NAME OF CALLING PROGRAM                                
C             NFILE(3) = NAME OF DISC FILE WITH FMP ERROR                       
C             ID = FMP OPERATION RESULTING IN ERROR                             
C             IER = ERROR CODE                                                  
C$                                                                              
      JERR = IER                                                                
      IF(IER .GE. 0) RETURN                                                     
      WRITE(LUCON,1001) NPROG,NFILE,IER,ID                                      
      RETURN                                                                    
 1001 FORMAT(1X,3A2,' - FILE ',3A2,' HAS ERROR ',I6,' ON ',A2)                  
      END                                                                       
      SUBROUTINE XYEAR(IYEAR,JYEAR,XDAY)                                        
C$                                                                              
C NAME        XYEAR                                                             
C                                                                               
C TYPE        SUBROUTINE                                                        
C                                                                               
C PURPOSE     GIVEN XDAY REFERRED TO JYEAR, REFER IT TO IYEAR                   
C                                                                               
C CALLING     CALL XYEAR(IYEAR,JYEAR,XDAY)                                      
C$                                                                              
      IF(IYEAR .EQ. JYEAR) RETURN                                               
      NYEAR = IABS(IYEAR - JYEAR)                                               
      KYEAR = JYEAR                                                             
      IF(IYEAR .LT. JYEAR) KYEAR = IYEAR                                        
      IF(IYEAR .LT. JYEAR) XDAY = -XDAY                                         
      DO 10 I = 1,NYEAR                                                         
        XDAY = XDAY + 365.                                                      
        IF(MOD(KYEAR + I - 1,4) .EQ. 0) XDAY = XDAY + 1.                        
   10   CONTINUE                                                                
      RETURN                                                                    
      END                                                                       
      SUBROUTINE DATE(IY,ND,M,ID,IW)                                            
      DIMENSION MONTH(12)                                                       
      DATA      MONTH  /31,28,31,30,31,30,31,31,30,31,30,31/                    
C$                                                                              
C NAME        DATE                                                              
C                                                                               
C TYPE        SUBROUTINE                                                        
C                                                                               
C PURPOSE     GIVEN YEAR AND DAY OF YEAR,                                       
C             FIND MONTH, DAY OF MONTH, DAY OF WEEK                             
C$                                                                              
C CHECK LEAP YEAR                                                               
   10 MONTH(2) = 28                                                             
      IF(MOD(IY,4) .EQ. 0) MONTH(2) = 29                                        
      LEAP = 0                                                                  
      IF(MOD(IY,4) .EQ. 0) LEAP = 1                                             
C MODIFY IY AND ND IF ND PAST END OF YEAR                                       
      IF(ND .LE. 365 + LEAP) GO TO 20                                           
      IY = IY + 1                                                               
      ND = ND - 365 - LEAP                                                      
      GO TO 10                                                                  
C MODIFY IY AND ND IF ND BEFORE START OF YEAR                                   
   20 IF(ND .GT. 0) GO TO 30                                                    
      IY = IY - 1                                                               
      ND = ND + 365                                                             
      IF(MOD(IY,4) .EQ. 0) ND = ND + 1                                          
      GO TO 10                                                                  
C SUBTRACT MONTHS UNTIL ND NEGATIVE.                                            
C   THEN NUMBER OF SUBTRACTIONS = NUMBER OF CURRENT MONTH                       
   30 ND1 = ND                                                                  
      DO 40 I = 1,12                                                            
        IM = I                                                                  
        ND1 = ND1 - MONTH(IM)                                                   
        IF(ND1 .LE. 0) GO TO 50                                                 
   40   CONTINUE                                                                
   50 M = IM                                                                    
C COMPUTE DAY OF WEEK. DAY ZERO FOR 1971 IS THU (IW=4)                          
      ID = ND1 + MONTH(IM)                                                      
      JY = 1971                                                                 
      IW = 4 + (IY - JY)                                                        
      IF(IY .GT. JY) IW = IW + (IY - JY + 2) / 4                                
      IF(IY .LT. JY) IW = IW + (IY - JY - 1) / 4                                
      IW = IW + ND                                                              
      IW = MOD(IW-1,7) + 1                                                      
      RETURN                                                                    
      END                                                                       
      SUBROUTINE TIME(DAYS,ITIME)                                               
      DOUBLE PRECISION DAYS,DBLE,SCALE,XDAY                                     
      DIMENSION ITIME(5),SCALE(4)                                               
      DATA SCALE/100D0,2*60D0,24D0/                                             
C                                                                               
C  NAME       TIME                                                              
C                                                                               
C  TYPE       SUBROUTINE                                                        
C                                                                               
C PURPOSE     CONVERTS TIME IN DAYS TO RTE TIME VECTOR                          
C                                                                               
C AUTHOR/DATE D. WELLS 780605                                                   
C                                                                               
C EXTERNALS   DBLE,FLOAT,IDINT                                                  
C                                                                               
C CALLING     CALL TIME(DAYS,ITIME)                                             
C                                                                               
C PARAMETERS  DAYS = DOUBLE PRECISION TIME IN DAYS                              
C             ITIME(5) = EQUIVALENT RTE TIME VECTOR                             
C                                                                               
C LANGUAGE    FORTRAN                                                           
C                                                                               
      XDAY = DAYS                                                               
      DO 5 I=1,4                                                                
         J=5-I                                                                  
         ITIME(J+1) = IDINT(XDAY)                                               
         XDAY=(XDAY-DBLE(FLOAT(ITIME(J+1)))) * SCALE(J)                         
    5    CONTINUE                                                               
      ITIME(1)=IDINT(XDAY)                                                      
      RETURN                                                                    
      END                                                                       
      DOUBLE PRECISION FUNCTION DAY(ITIME)                                      
      DOUBLE PRECISION DBLE,SCALE,XDAY                                          
      DIMENSION ITIME(5),SCALE(4)                                               
      DATA SCALE/100D0,2*60D0,24D0/                                             
C                                                                               
C NAME        DAY                                                               
C                                                                               
C TYPE        DOUBLE PRECISION FUNCTION                                         
C                                                                               
C PURPOSE     CONVERT RTE TIME VECTOR TO TIME IN DAYS                           
C                                                                               
C AUTHOR/DATE D. WELLS 780605                                                   
C                                                                               
C EXTERNALS DBLE,FLOAT                                                          
C                                                                               
C CALLING     X = DAY(ITIME)                                                    
C                                                                               
C PARAMETERS  ITIME(5) = RTE TIME VECTOR                                        
C             DAY = EQUIVALENT TIME IN DAYS                                     
C                                                                               
C LANGUAGE    FORTRAN                                                           
C                                                                               
C                                                                               
      XDAY = 0.1D0                                                              
      DO 5 I=1,4                                                                
    5   XDAY=(XDAY+DBLE(FLOAT(ITIME(I)))) / SCALE(I)                            
      XDAY=XDAY+DBLE(FLOAT(ITIME(5)))                                           
      DAY = XDAY                                                                
      RETURN                                                                    
      END                                                                       
