text:
      INTEGER FUNCTION GVRINI(IFUNC,IPARMS)

C *** McIDAS Revision History ***
C 1 NVXGVAR.DLM 5-Jul-94,13:01:32,`KENB' Initial release of GVAR navigation
C      module. (4663)
C 2 NVXGVAR.DLM 25-Aug-94,13:30:14,`USER' Released
C 3 NVXGVAR.DLM 29-Aug-94,10:56:58,`RUSSN' Corrected revision history cards
C 4 NVXGVAR.DLM 13-Sep-94,12:09:06,`USER' Released
C 5 NVXGVAR.DLM 28-Oct-94,15:34:28,`KENB' adopted fix from m/f for changes
C 6 NVXGVAR.DLM 31-Oct-94,12:12:46,`USER' Released
C 7 NVXGVAR.DLM 19-Feb-96,16:07:56,`DWS' reglue: modified file
C 8 NVXGVAR.DLM 20-Feb-96,12:02:20,`USER' Released
C 9 NVXGVAR.DLM 11-Jul-96,14:36:20,`JMB' comment out ddests (6616)
C 10 NVXGVAR.DLM 29-Jul-96,12:57:46,`USER' Released
C *** McIDAS Revision History ***
C J. Cowie/COMET	 1/95	Changed names of routines to avoid conflicts
C				with GOES nav
C
C				NVXINI --> GVRINI
C      				NVXSAE --> GVRSAE
C				NVXEAS --> GVREAS
C				NVXOPT --> GVROPT
      IMPLICIT NONE

C     INCLUDE 'ELCONS.INC'
C=========================== DIELCONS =============================

C $   (JR)
C $   THIS INCLUDE FILE IS PART OF THE SUPPLIED GVAR NAV SOFTWARE
C $
C $   DESCRIPTION:
C $
C $
C $$  DIELCONS = INCLUDE

      REAL*8 PI
           PARAMETER (PI=3.141592653589793D0)
      REAL*8 DEG
           PARAMETER (DEG=180.D0/PI)
      REAL*8 RAD
           PARAMETER (RAD=PI/180.D0)
C                    DEGREES TO RADIANS CONVERSION PI/180
      REAL*8 NOMORB
           PARAMETER (NOMORB=42164.365D0)
C                    NOMINAL RADIAL DISTANCE OF SATELLITE (KM)
      REAL*8 AE
           PARAMETER (AE=6378.137D0)
C                    EARTH EQUATORIAL RADIUS (KM)
      REAL*8 FER
           PARAMETER (FER=1.D0-(6356.7533D0/AE))
C                    EARTH FLATTENING COEFFICIENT = 1-(BE/AE)
      REAL*4 AEBE2
           PARAMETER (AEBE2=1.D0/(1.D0-FER)**2)
      REAL*4 AEBE3
           PARAMETER (AEBE3=AEBE2-1.)
      REAL*4 AEBE4
           PARAMETER (AEBE4=(1.D0-FER)**4-1.)
C========================= END DIELCONS ===========================
C     INCLUDE 'ELCOMM.INC'
C=========================== DIELCOMM =============================

C $   (JR)
C $   THIS INCLUDE FILE IS PART OF THE SUPPLIED GVAR NAV SOFTWARE
C $
C $   DESCRIPTION:
C $
C $
C $$  DIELCOMM = INCLUDE

C     COMMON VARIABLES
C
      REAL*8 XS(3)
C                      NORMALIZED S/C POSITION IN ECEF COORDINATES
      REAL*8 BT(3,3)
C                      ECEF TO INSTRUMENT COORDINATES TRANSFORMATION
      REAL*8  Q3
C                      USED IN SUBRTN LPOINT
      REAL*8 PITCH,ROLL,YAW
C                          PITCH,ROLL,YAW ANGLES OF INSTRUMENT (RAD)
      REAL*4 PMA,RMA
C                          PITCH,ROLL MISALIGNMENTS OF INSTRUMENT (RAD)
         COMMON /ELCOMM/ XS,BT,Q3,PITCH,ROLL,YAW,PMA,RMA
C========================= END DIELCOMM ===========================
C     INCLUDE 'INSTCO.INC'
C=========================== DIINSTCO =============================

C $   (JR)
C $   THIS INCLUDE FILE IS PART OF THE SUPPLIED GVAR NAV SOFTWARE
C $
C $   DESCRIPTION:
C $
C $
C $$  DIINSTCO = INCLUDE

      INTEGER*4 INCMAX(2)
C                       NUMBER OF INCREMENTS PER CYCLE
      REAL*4 ELVMAX(2)
C                       BOUNDS IN ELEVATION (RADIANS)
      REAL*4 SCNMAX(2)
C                       BOUNDS IN SCAN ANGLE (RADIANS)
      REAL*4 ELVINC(2)
C                       CHANGE IN ELEVATION ANGLE PER INCREMENT (RAD)
      REAL*4 SCNINC(2)
C                       CHANGE IN SCAN ANGLE PER INCREMENT (RADIANS)
      REAL*4 ELVLN(2)
C                       ELEVATION ANGLE PER DETECTOR LINE (RADIANS)
      REAL*4 SCNPX(2)
C                       SCAN ANGLE PER PIXEL (RADIANS)

      COMMON /INSTCO/ INCMAX,ELVMAX,SCNMAX,
     1   ELVINC,SCNINC,ELVLN,SCNPX
C========================= END DIINSTCO ===========================
C     INCLUDE 'GVRCOD.INC'
C=========================== DIGVRCOD =============================

C $   (JR)
C $   THIS INCLUDE FILE IS PART OF THE GVAR NAVIGATION SUBSYSTEM.
C $
C $   DESCRIPTION:
C $      THESE ARE CONSTANTS USED TO DESCRIBE THE GVAR CODICIL STRUCTURE
C $
C $$  DIGVRCOD = INCLUDE, NAVIGATION, GVAR

C DEFINITIONS OF PARAMETERS
C STTYPE   : POSITION OF SATELLITE TYPE
C ETIME    : POSITION OF EPOCH TIME
C IMCACT   : POSITION OF IMC ACTIVE FLAG
C IMGDAY   : POSITION OF IMAGE DAY VALUE  (YYDDD)
C IMGTM    : POSITION OF IMAGE TIME VALUE (HHMMSS)
C REFLON   : POSITION OF REFERENCE LONGITUDE
C REFDIS   : POSITION OF REFERENCE DISTANCE FROM NOMINAL
C REFLAT   : POSITION OF REFERENCE LATITUDE
C REFYAW   : POSITION OF REFERENCE YAW
C RATROL   : POSITION OF REFERENCE ATTITUDE ROLL
C RATPTC   : POSITION OF REFERENCE ATTITUDE PITCH
C RATYAW   : POSITION OF REFERENCE ATTITUDE YAW
C LDR1-13  : LOCATION OF LONGITUDE DELTA FROM REFERENCE PARAMETERS
C RDDR1-11 : LOCATION OF RADIAL DISTANCE DELTA FROM REFERENCE PARAMETERS
C DGL1-9   : LOCATION OF GEOCENTRIC LATITUDE DELTA PARAMETERS
C DOY1-9   : LOCATION OF ORBIT YAW DELTA PARAMETERS
C EXPTIM   : EXPONENTIAL START TIME FROM EPOCH
C RAAWDS   : LOCATION OF START OF ROLL ATTITUDE ANGLE WORDS
C PAAWDS   : LOCATION OF START OF PITCH ATTITUDE ANGLE WORDS
C YAAWDS   : LOCATION OF START OF YAW ATTITUDE ANGLE WORDS
C RMAWDS   : LOCATION OF START OF ROLL MISALIGNMENT ANGLE WORDS
C PMAWDS   : LOCATION OF START OF PITCH MISALIGNMENT ANGLE WORDS
C EDTIME   : LOCATION OF DELTA FROM EPOCH TIME
C IMCROL   : LOCATION OF IMAGE MOTION COMPENSATION ROLL
C IMCPTC   : LOCATION OF IMAGE MOTION COMPENSATION PITCH
C IMCYAW   : LOCATION OF IMAGE MOTION COMPENSATION YAW
C IMGSND   : LOCATION OF IMAGER/SOUNDER INSTRUMENT FLAG

      INTEGER      ETIME , IMCACT, IMGDAY, IMGTM , RATROL, RATPTC
      INTEGER      RATYAW, LDR1  , LDR2  , LDR3  , LDR4  , LDR5
      INTEGER      LDR6  , LDR7  , LDR8  , LDR9  , LDR10 , LDR11
      INTEGER      LDR12 , LDR13 , RDDR1 , RDDR2 , RDDR3 , RDDR4
      INTEGER      RDDR5 , RDDR6 , RDDR7 , RDDR8 , RDDR9 , RDDR10
      INTEGER      RDDR11, DGL1  , DGL2  , DGL3  , DGL4  , DGL5
      INTEGER      DGL6  , DGL7  , DGL8  , DGL9  , DOY1  , DOY2
      INTEGER      DOY3  , DOY4  , DOY5  , DOY6  , DOY7  , DOY8
      INTEGER      DOY9  , EXPTIM, RAAWDS, PAAWDS, YAAWDS, RMAWDS
      INTEGER      PMAWDS, EDTIME, IMCROL, IMCPTC, IMCYAW, IMGSND
      INTEGER      REFLON, REFDIS, REFLAT, REFYAW, STTYPE, IDNTFR
      INTEGER      IOFNC , IOFNI , IOFEC , IOFEI

      PARAMETER (STTYPE = 1)
      PARAMETER (IDNTFR = 2)
      PARAMETER (IMCACT = 3)
      PARAMETER (REFLON = 6)
      PARAMETER (REFDIS = 7)
      PARAMETER (REFLAT = 8)
      PARAMETER (REFYAW = 9)
      PARAMETER (RATROL = 10)
      PARAMETER (RATPTC = 11)
      PARAMETER (RATYAW = 12)
      PARAMETER (ETIME  = 13)
      PARAMETER (EDTIME = 15)
      PARAMETER (IMCROL = 16)
      PARAMETER (IMCPTC = 17)
      PARAMETER (IMCYAW = 18)
      PARAMETER (LDR1   = 19)
      PARAMETER (LDR2   = 20)
      PARAMETER (LDR3   = 21)
      PARAMETER (LDR4   = 22)
      PARAMETER (LDR5   = 23)
      PARAMETER (LDR6   = 24)
      PARAMETER (LDR7   = 25)
      PARAMETER (LDR8   = 26)
      PARAMETER (LDR9   = 27)
      PARAMETER (LDR10  = 28)
      PARAMETER (LDR11  = 29)
      PARAMETER (LDR12  = 30)
      PARAMETER (LDR13  = 31)
      PARAMETER (RDDR1  = 32)
      PARAMETER (RDDR2  = 33)
      PARAMETER (RDDR3  = 34)
      PARAMETER (RDDR4  = 35)
      PARAMETER (RDDR5  = 36)
      PARAMETER (RDDR6  = 37)
      PARAMETER (RDDR7  = 38)
      PARAMETER (RDDR8  = 39)
      PARAMETER (RDDR9  = 40)
      PARAMETER (RDDR10 = 41)
      PARAMETER (RDDR11 = 42)
      PARAMETER (DGL1   = 43)
      PARAMETER (DGL2   = 44)
      PARAMETER (DGL3   = 45)
      PARAMETER (DGL4   = 46)
      PARAMETER (DGL5   = 47)
      PARAMETER (DGL6   = 48)
      PARAMETER (DGL7   = 49)
      PARAMETER (DGL8   = 50)
      PARAMETER (DGL9   = 51)
      PARAMETER (DOY1   = 52)
      PARAMETER (DOY2   = 53)
      PARAMETER (DOY3   = 54)
      PARAMETER (DOY4   = 55)
      PARAMETER (DOY5   = 56)
      PARAMETER (DOY6   = 57)
      PARAMETER (DOY7   = 58)
      PARAMETER (DOY8   = 59)
      PARAMETER (DOY9   = 60)
      PARAMETER (EXPTIM = 62)
      PARAMETER (RAAWDS = 63)
      PARAMETER (PAAWDS = 130)
      PARAMETER (YAAWDS = 185)
      PARAMETER (RMAWDS = 258)
      PARAMETER (PMAWDS = 313)

      PARAMETER (IMGDAY = 368)
      PARAMETER (IMGTM  = 369)
      PARAMETER (IMGSND = 370)

C THESE FOUR WORDS WERE ADDED 5-26-94 TO COMPLY W/ THE NEW ELUG
C NUMBERING STARTED AT 380 BECAUSE THESE SAME PARAMETERS ARE USED
C IN THE NAV MESSAGE SENT FROM THE INGESTOR TO EVX, AND WE HAD TO
C START SOMEWHERE AFTER THE 378 NAV PARAMETERS
      PARAMETER (IOFNC = 380)
      PARAMETER (IOFEC = 381)
      PARAMETER (IOFNI = 382)
      PARAMETER (IOFEI = 383)

C========================= END DIGVRCOD ===========================
C     INCLUDE 'MXCDSZ.INC'
C=========================== DIMXCDSZ =============================

C $   (JR)
C $   THIS INCLUDE FILE IS PART OF THE NAVIGATION SUBSYSTEM
C $
C $   DESCRIPTION:
C $      THIS IS THE VALUE OF LARGEST NAVIGATION CODICIL SIZE ALLOWED
C $
C $$  DIMXCDSZ = INCLUDE, NAVIGATION

C DEFINITIONS OF PARAMETERS
C MXCDSZ   : LARGEST ALLOWED CODICIL SIZE

      INTEGER      MXCDSZ

      PARAMETER (MXCDSZ = 5*128)

C========================= END DIMXCDSZ ===========================
C     INCLUDE 'OASET.INC'
C=========================== DIOASET =============================

C $   (JR)
C $   THIS INCLUDE FILE IS PART OF THE GVAR NAVIGATION SUBSYSTEM.
C $
C $   DESCRIPTION:
C $      THESE ARE CONSTANTS USED TO DESCRIBE THE GVAR O&A SET STRUCTURE
C $
C $$  DIOASET = INCLUDE, NAVIGATION, GVAR

C DEFINITIONS OF PARAMETERS
C OASIZE   : SIZE OF GVAR ORBIT AND ATTITUDE SET
C PCOEFS   : START OF PITCH COEFFICIENTS
C RMACFS   : START OF RMA COEFFICIENTS
C CUTOF1   : FIRST DIVIDING POINT IN O&A SET (WHERE DIVIDED INTO 128
C            WORD SETS)
C CUTOF2   : SECOND DIVIDING POINT IN O&A SET (WHERE DIVIDED INTO 128
C            WORD SETS)

      INTEGER      OASIZE, PCOEFS, RMACFS, CUTOF1, CUTOF2

      PARAMETER (OASIZE = 336)
      PARAMETER (PCOEFS = 117)
      PARAMETER (RMACFS = 227)
      PARAMETER (CUTOF1 = 116)
      PARAMETER (CUTOF2 = 226)

C========================= END DIOASET ===========================
C     INCLUDE 'B11DOC.INC'
C=========================== DIB11DOC =============================

C $   (JR)
C $   THIS INCLUDE FILE IS PART OF THE DECODER SCHEDULER SUBSYSTEM.
C $
C $   DESCRIPTION:
C $      THESE ARE INDICES INTO THE SOUNDER DOCUMENTATION BLOCK
C $      (INDICES ARE 1-BASED FROM START OF SAD ID)
C $
C $$  DIB11DOC = INCLUDE, DECODER

C DEFINITIONS OF PARAMETERS
C     SSCAN   SOUNDER SCAN STATUS IN TWO BYTES
C     SGBCT   TOTAL NUMBER OF BLOCK 11'S IN THIS SCAN (INCLUDES DOC)
C     TSCLS   TIME OF SCAN LINE START IN 8 BYTES
C     NEWIMG  BIT # IN SSCAN FOR NEW IMAGE FLAG
C     ENDIMG  BIT # IN SSCAN FOR END OF IMAGE FLAG
C     WARN    BIT # IN SSCAN FOR WARNING MODE FLAG
C     DWELL1  BIT # IN SSCAN FOR DWELL MODE 1 FLAG
C     DWELL2  BIT # IN SSCAN FOR DWELL MODE 2 FLAG
C     DWELL4  BIT # IN SSCAN FOR DWELL MODE 4 FLAG
C     IMCFLG  BIT # IN SSCAN FOR IMC ACTIVE FLAG
C     STPLOC  BIT # IN SSCAN FOR STEP MODE FLAG
C     TIMSIZ  NUMBER OF BYTES THAT THE TIME OCCUPIES
C     SNSLN   NORTHERN MOST SCAN LINE OF CURRENT SCAN
C     SWFPX   WESTERNMOST PIXEL # OF CURRENT FRAME
C     SEFPX   EASTERNMOST PIXEL # OF CURRENT FRAME
C     SNFLN   NORTHERNMOST SCAN LINE # OF CURRENT FRAME
C     SSFLN   SOUTHERNMOST SCAN LINE # OF CURRENT FRAME
C     NAVLOC  START OF O&A PARAMETERS
C     NAVEND  END OF O&A PARAMETERS

      INTEGER      SSCAN , SGBCT , TSCLS , NEWIMG, ENDIMG, TIMSIZ
      INTEGER      SWFPX , SEFPX , SNFLN , SSFLN , SDSTA , DWELL1
      INTEGER      NAVLOC, DWELL2, DWELL4, SNSLN , WARN  , STPLOC
      INTEGER      NAVEND, IMCFLG

C LOCATION OF FIELDS WITHIN DOCUMENTATION BLOCK
      PARAMETER (SSCAN  = 31)
      PARAMETER (SDSTA  = 33)
      PARAMETER (SGBCT  = 45)
      PARAMETER (TSCLS  = 67)
      PARAMETER (NEWIMG = 15)
      PARAMETER (ENDIMG = 14)
      PARAMETER (WARN   = 11)
      PARAMETER (STPLOC = 3)
      PARAMETER (DWELL1 = 4)
      PARAMETER (DWELL2 = 5)
      PARAMETER (DWELL4 = 6)
      PARAMETER (IMCFLG = 7)
      PARAMETER (TIMSIZ = 8)
      PARAMETER (SNSLN  = 185)
      PARAMETER (SWFPX  = 187)
      PARAMETER (SEFPX  = 189)
      PARAMETER (SNFLN  = 191)
      PARAMETER (SSFLN  = 193)
      PARAMETER (NAVLOC = 307)
      PARAMETER (NAVEND = 1718)

C========================= END DIB11DOC ===========================

      REAL*8  TE    , PHI   , PSI   , U     , SINU  , COSU
      REAL*8  SINOI , COSOI , SLAT  , ASC   , SINASC, COSASC
      REAL*8  SYAW  , W     , SW    , CW    , S2W   , C2W
      REAL*8  SW1   , CW1   , SW3   , CW3   , GATT
      REAL*8  SUBLAT, SUBLON, SECS  , WA
      REAL    RPARMS(MXCDSZ)
      REAL*8  IMGTIM, EPOCH , TIMEX , TIME50, R     , B(3,3)
      REAL*8  AEC   , TS    , DR    , LAM   , DLAT  , DYAW
      REAL*8  AEBE2C, AEBE3C, AEBE4C, FERC
      INTEGER      IFUNC , LIT   , ITYPE , IMC   , YEAR  , DAY
      INTEGER      IFTOK , HOUR  , MINUTS, INSTR , IPARMS(MXCDSZ)
      INTEGER      COUNT , OFFSET, LOOP  , IPARM2(MXCDSZ), TIME(2)
      LOGICAL      BTEST
      CHARACTER*4  CLIT
      CHARACTER*12 CFZ

      EQUIVALENCE(IPARM2,RPARMS)

      COMMON /GVRCOM/ ITYPE , INSTR , SUBLAT, SUBLON
      COMMON /RADCOM/ AEC   , FERC  , AEBE2C, AEBE3C, AEBE4C
      COMMON /SAVCOM/ B     , DR    , PHI

C     INCLUDE 'GVRDAT.INC'
C=========================== DIGVRDAT =============================

C $   (JR)
C $   THIS INCLUDE FILE IS PART OF THE GVAR NAVIGATION SUBSYSTEM.
C $
C $   DESCRIPTION:
C $      THESE ARE DATA STATEMENTS USED IN THE GVAR NAVIGATION
C $
C $$  DIGVRDAT = INCLUDE, NAVIGATION, GVAR

C DEFINITIONS OF PARAMETERS
C RELLST   : TWO DIMENSIONAL ARRAY OF LOCATION OF REALS IN O&A SET

      INTEGER      I     , J     , RELLST(40,2)

      DATA ((RELLST(I,J),J=1,2),I=1,31) / 5,11,14,64,66,95,99,101,104,
     &        106,109,111,114,116,117,119,121,150,154,156,159,161,164,
     &        166,169,171,172,174,176,205,209,211,214,216,219,221,224,
     &        226,227,229,231,260,264,266,269,271,274,276,279,284,286,
     &        315,319,321,324,326,329,331,334,336,-1,-1 /

C========================= END DIGVRDAT ===========================

C--------------------------- NVXGVR ------------------------------

C INITIALIZE THE EARTH RADIUS CONSTANTS
      AEC    = AE
      FERC   = FER
      AEBE2C = DBLE(AEBE2)
      AEBE3C = DBLE(AEBE3)
      AEBE4C = DBLE(AEBE4)

C IFUNC = 1
      IF (IFUNC.EQ.1) THEN
        IF (IPARMS(STTYPE).NE.LIT('GVAR')) THEN
          GVRINI = -1
          RETURN
        ENDIF
        ITYPE = 1

C COPY CODICIL SO WE CAN EQUIVALENCE IT TO REAL ARRAY
        DO 1 LOOP = 1, MXCDSZ
 1        IPARM2(LOOP) = IPARMS(LOOP)

C SCALE DOWN REAL CODICIL ELEMENTS
        COUNT = 1
        RPARMS(IMGTM) = REAL(IPARM2(IMGTM))/1000.
 2      IF (RELLST(COUNT,1).EQ.-1) GOTO 4
          OFFSET = 1
          IF (RELLST(COUNT,1).GT.CUTOF1) OFFSET = 13
          IF (RELLST(COUNT,1).GT.CUTOF2) OFFSET = 31
          DO 3 LOOP=RELLST(COUNT,1),RELLST(COUNT,2)
            IF (LOOP.EQ.14.OR.LOOP.EQ.61.OR.(MOD(LOOP-8,55).EQ.0.AND.
     &                                               LOOP.NE.8)) THEN
              RPARMS(LOOP+OFFSET) = REAL(IPARM2(LOOP+OFFSET))/100.
            ELSE
              RPARMS(LOOP+OFFSET) = REAL(IPARM2(LOOP+OFFSET))/10000000.
            ENDIF
 3        CONTINUE
          COUNT = COUNT + 1
        GOTO 2

C SEE IF THIS CODICIL IS FOR IMAGER OR SOUNDER
 4      INSTR = IPARMS(IMGSND)

C INITIALIZE CONSTANTS IN COMMON INSTCO
        CALL SETCON(INSTR, IPARMS(IOFNC), IPARMS(IOFNI), IPARMS(IOFEC),
     &                                                  IPARMS(IOFEI))

C GET CONTROL INFO FROM CODICIL
        YEAR    = 1900 + IPARMS(IMGDAY) / 1000
        DAY     = IPARMS(IMGDAY) - IPARMS(IMGDAY) / 1000 * 1000
        HOUR    = RPARMS(IMGTM) / 10000
        MINUTS  = RPARMS(IMGTM) / 100 - HOUR * 100
        SECS    = RPARMS(IMGTM) - REAL(100*MINUTS) - REAL(10000*HOUR)
        IMGTIM  = TIMEX(YEAR, DAY, HOUR, MINUTS, SECS)
        TIME(1) = IFTOK(CFZ(IPARMS(ETIME)))
        TIME(2) = IFTOK(CFZ(IPARMS(ETIME+1)))
        EPOCH   = TIME50(TIME)
        IMC     = 1
        IF (BTEST(IPARMS(IMCACT),IMCFLG)) IMC = 0

C ASSIGN REFERENCE VALUES TO THE SUBSATELLITE LONGITUDE AND
C LATITUDE, THE RADIAL DISTANCE AND THE ORBIT YAW.
        LAM = RPARMS(REFLON)
        DR  = RPARMS(REFDIS)
        PHI = RPARMS(REFLAT)
        PSI = RPARMS(REFYAW)

C ASSIGN REFERENCE VALUES TO THE ATTITUDES AND MISALIGNMENTS
        ROLL  = RPARMS(RATROL)
        PITCH = RPARMS(RATPTC)
        YAW   = RPARMS(RATYAW)
        RMA   = 0.
        PMA   = 0.

C IF IMC IS OFF, COMPUTE CHANGES IN THE SATELLITE ORBIT
        IF (IMC.NE.0) THEN

C SET REFERENCE RADIAL DISTANCE, LATITUDE AND ORBIT YAW TO ZERO
          DR  = 0.
          PHI = 0.
          PSI = 0.

C COMPUTE TIME SINCE EPOCH (IN MINUTES)
          TS = IMGTIM - EPOCH

C COMPUTES ORBIT ANGLE AND THE RELATED TRIGONOMETRIC FUNKTIONS.
C EARTH ROTATIONAL RATE=.729115E-4 (RAD/S)
          W   = 0.729115D - 4 * 60.0D0 * TS
          SW  = DSIN(W)
          CW  = DCOS(W)
          SW1 = DSIN(0.927*W)
          CW1 = DCOS(0.927*W)
          S2W = DSIN(2.*W)
          C2W = DCOS(2.*W)
          SW3 = DSIN(1.9268*W)
          CW3 = DCOS(1.9268*W)

C COMPUTES CHANGE IN THE IMC LONGITUDE FROM THE REFERENCE
          LAM = LAM + RPARMS(LDR1) + (RPARMS(LDR2) + RPARMS(LDR3)*W) * W
     &        + (RPARMS(LDR10)*SW1 + RPARMS(LDR11)*CW1 + RPARMS(LDR4)*SW
     &        + RPARMS(LDR5)*CW + RPARMS(LDR6)*S2W + RPARMS(LDR7)*C2W
     &        + RPARMS(LDR8)*SW3+RPARMS(LDR9)*CW3 + W*(RPARMS(LDR12)*SW
     &        + RPARMS(LDR13)*CW))*2.

C COMPUTES CHANGE IN RADIAL DISTANCE FROM THE REFERENCE (KM)
          DR = DR + RPARMS(RDDR1) + RPARMS(RDDR2)*CW + RPARMS(RDDR3)*SW
     &          + RPARMS(RDDR4)*C2W + RPARMS(RDDR5)*S2W + RPARMS(RDDR6)
     &          * CW3+RPARMS(RDDR7)*SW3 + RPARMS(RDDR8)*CW1
     &          + RPARMS(RDDR9)*SW1 + W*(RPARMS(RDDR10)*CW
     &          + RPARMS(RDDR11)*SW)

C COMPUTES THE SINE OF THE CHANGE IN THE GEOCENTRIC LATITUDE
          DLAT = RPARMS(DGL1) + RPARMS(DGL2)*CW + RPARMS(DGL3)*SW
     &        + RPARMS(DGL4)*C2W + RPARMS(DGL5)*S2W + W*(RPARMS(DGL6)*CW
     &        + RPARMS(DGL7)*SW) + RPARMS(DGL8)*CW1+RPARMS(DGL9)*SW1

C COMPUTES GEOCENTRIC LATITUDE BY USING AN EXPANSION FOR ARCSINE
          PHI = PHI + DLAT * (1. + DLAT * DLAT / 6.)

C COMPUTES SINE OF THE CHANGE IN THE ORBIT YAW
          DYAW = RPARMS(DOY1) + RPARMS(DOY2)*SW + RPARMS(DOY3)*CW
     &         + RPARMS(DOY4)*S2W + RPARMS(DOY5)*C2W
     &         + W*(RPARMS(DOY6)*SW + RPARMS(DOY7)*CW)
     &         + RPARMS(DOY8)*SW1 + RPARMS(DOY9)*CW1

C COMPUTES THE ORBIT YAW BY USING AN EXPANSION FOR ARCSINE.
          PSI = PSI + DYAW * (1. + DYAW * DYAW / 6.)

C CALCULATION OF CHANGES IN THE SATELLITE ORBIT ENDS HERE
        ENDIF

C CONVERSION OF THE IMC LONGITUDE AND ORBIT YAW TO THE SUBSATELLITE
C LONGITUDE AND THE ORBIT INCLINATION (REF: GOES-PCC-TM-2473, INPUTS
C REQUIRED FOR EARTH LOCATION AND GRIDDING BY SPS, JUNE 6, 1988)
        SLAT  = DSIN(PHI)
        SYAW  = DSIN(PSI)
        SINOI = SLAT**2 + SYAW**2
        COSOI = DSQRT(1.-SINOI)
        SINOI = DSQRT(SINOI)

        IF (SLAT .EQ. 0.0D0 .AND. SYAW .EQ. 0.0D0) THEN
          U = 0.0D0
        ELSE
          U = DATAN2(SLAT,SYAW)
        ENDIF

        SINU  = DSIN(U)
        COSU  = DCOS(U)

C COMPUTES LONGITUDE OF THE ASCENDING NODE
        ASC    = LAM-U
        SINASC = DSIN(ASC)
        COSASC = DCOS(ASC)

C COMPUTES THE SUBSATELLITE GEOGRAPHIC LATITUDE
        SUBLAT = DATAN(AEBE2C*DTAN(PHI))

C COMPUTES THE SUBSATELLITE LONGITUDE
        SUBLON = ASC+DATAN2(COSOI*SINU,COSU)

C COMPUTES THE SPACECRAFT TO EARTH FIXED COORDINATES TRANSFORMATION
C MATRIX:
C     (VECTOR IN ECEF COORDINATES) = B * (VECTOR IN S/C COORDINATES)
        B(1,2) = (-SINASC)*SINOI
        B(2,2) =  COSASC*SINOI
        B(3,2) = -COSOI
        B(1,3) = (-COSASC)*COSU+SINASC*SINU*COSOI
        B(2,3) = (-SINASC)*COSU-COSASC*SINU*COSOI
        B(3,3) = -SLAT
        B(1,1) = (-COSASC)*SINU-SINASC*COSU*COSOI
        B(2,1) = (-SINASC)*SINU+COSASC*COSU*COSOI
        B(3,1) =  COSU*SINOI

C COMPUTES THE NORMALIZED SPACECRAFT POSITION VECTOR IN EARTH FIXED
C COORDINATES - XS.
        R     = (NOMORB+DR)/AEC
        XS(1) = (-B(1,3))*R
        XS(2) = (-B(2,3))*R
        XS(3) = (-B(3,3))*R

C PRECOMPUTES Q3 (USED IN LPOINT)
        Q3 = XS(1)**2 + XS(2)**2 + AEBE2C * XS(3)**2 - 1.0

C COMPUTES THE ATTITUDES AND MISALIGNMENTS IF IMC IS OFF
        IF (IMC.NE.0) THEN

C COMPUTES THE SOLAR ORBIT ANGLE
          WA = RPARMS(61)*TS

C COMPUTES THE DIFFERENCE BETWEEN CURRENT TIME, TS, AND THE
C EXPONENTIAL TIME, IPARMS(62). NOTE THAT BOTH TIMES ARE SINCE EPOCH.
          TE = TS - RPARMS(EXPTIM)

C COMPUTES ROLL + ROLL MISALIGNMENT
          ROLL = ROLL + GATT(RAAWDS,RPARMS,IPARMS,WA,TE)

C COMPUTES PITCH + PITCH MISALIGNMENT
          PITCH = PITCH + GATT(PAAWDS,RPARMS,IPARMS,WA,TE)

C COMPUTES YAW
          YAW = YAW + GATT(YAAWDS,RPARMS,IPARMS,WA,TE)

C COMPUTES ROLL MISALIGNMENT
          RMA = GATT(RMAWDS,RPARMS,IPARMS,WA,TE)

C COMPUTES PITCH MISALIGNMENT
          PMA = GATT(PMAWDS,RPARMS,IPARMS,WA,TE)

C APPLY THE EARTH SENSOR COMPENSATION IF NEEDED
          ROLL   = ROLL + RPARMS(IMCROL)
          PITCH  = PITCH + RPARMS(IMCPTC)
          YAW    = YAW + RPARMS(IMCYAW)
        ENDIF

C COMPUTES THE INSTRUMENT TO EARTH FIXED COORDINATES TRANSFORMATION
C MATRIX - BT
        CALL INST2E(ROLL,PITCH,YAW,B,BT)

C IFUNC = 2
      ELSEIF (IFUNC.EQ.2) THEN
        IF (INDEX(CLIT(IPARMS(1)),'LL').NE.0) ITYPE=1
        IF (INDEX(CLIT(IPARMS(1)),'XY').NE.0) ITYPE=2
      ENDIF
      GVRINI=0
      RETURN
      END

      INTEGER FUNCTION GVRSAE(XLIN,XELE,XDUM,XLAT,XLON,Z)

      IMPLICIT NONE

      REAL     XLIN  , XELE  , XDUM  , XLAT  , XLON  , Z
      REAL*8   RL    , RP    , E     , S
      REAL     YLAT  , YLON
      REAL*8   SUBLAT, SUBLON, EVLN  , SCPX  , TMPLAT, TMPLON
      INTEGER      INSTR , ITYPE , STAT

C     INCLUDE 'ELCONS.INC'
C=========================== DIELCONS =============================

C $   (JR)
C $   THIS INCLUDE FILE IS PART OF THE SUPPLIED GVAR NAV SOFTWARE
C $
C $   DESCRIPTION:
C $
C $
C $$  DIELCONS = INCLUDE

      REAL*8 PI
           PARAMETER (PI=3.141592653589793D0)
      REAL*8 DEG
           PARAMETER (DEG=180.D0/PI)
      REAL*8 RAD
           PARAMETER (RAD=PI/180.D0)
C                    DEGREES TO RADIANS CONVERSION PI/180
      REAL*8 NOMORB
           PARAMETER (NOMORB=42164.365D0)
C                    NOMINAL RADIAL DISTANCE OF SATELLITE (KM)
      REAL*8 AE
           PARAMETER (AE=6378.137D0)
C                    EARTH EQUATORIAL RADIUS (KM)
      REAL*8 FER
           PARAMETER (FER=1.D0-(6356.7533D0/AE))
C                    EARTH FLATTENING COEFFICIENT = 1-(BE/AE)
      REAL*4 AEBE2
           PARAMETER (AEBE2=1.D0/(1.D0-FER)**2)
      REAL*4 AEBE3
           PARAMETER (AEBE3=AEBE2-1.)
      REAL*4 AEBE4
           PARAMETER (AEBE4=(1.D0-FER)**4-1.)
C========================= END DIELCONS ===========================
C     INCLUDE 'ELCOMM.INC'
C=========================== DIELCOMM =============================

C $   (JR)
C $   THIS INCLUDE FILE IS PART OF THE SUPPLIED GVAR NAV SOFTWARE
C $
C $   DESCRIPTION:
C $
C $
C $$  DIELCOMM = INCLUDE

C     COMMON VARIABLES
C
      REAL*8 XS(3)
C                      NORMALIZED S/C POSITION IN ECEF COORDINATES
      REAL*8 BT(3,3)
C                      ECEF TO INSTRUMENT COORDINATES TRANSFORMATION
      REAL*8  Q3
C                      USED IN SUBRTN LPOINT
      REAL*8 PITCH,ROLL,YAW
C                          PITCH,ROLL,YAW ANGLES OF INSTRUMENT (RAD)
      REAL*4 PMA,RMA
C                          PITCH,ROLL MISALIGNMENTS OF INSTRUMENT (RAD)
         COMMON /ELCOMM/ XS,BT,Q3,PITCH,ROLL,YAW,PMA,RMA
C========================= END DIELCOMM ===========================
C     INCLUDE 'OASET.INC'
C=========================== DIOASET =============================

C $   (JR)
C $   THIS INCLUDE FILE IS PART OF THE GVAR NAVIGATION SUBSYSTEM.
C $
C $   DESCRIPTION:
C $      THESE ARE CONSTANTS USED TO DESCRIBE THE GVAR O&A SET STRUCTURE
C $
C $$  DIOASET = INCLUDE, NAVIGATION, GVAR

C DEFINITIONS OF PARAMETERS
C OASIZE   : SIZE OF GVAR ORBIT AND ATTITUDE SET
C PCOEFS   : START OF PITCH COEFFICIENTS
C RMACFS   : START OF RMA COEFFICIENTS
C CUTOF1   : FIRST DIVIDING POINT IN O&A SET (WHERE DIVIDED INTO 128
C            WORD SETS)
C CUTOF2   : SECOND DIVIDING POINT IN O&A SET (WHERE DIVIDED INTO 128
C            WORD SETS)

      INTEGER      OASIZE, PCOEFS, RMACFS, CUTOF1, CUTOF2

      PARAMETER (OASIZE = 336)
      PARAMETER (PCOEFS = 117)
      PARAMETER (RMACFS = 227)
      PARAMETER (CUTOF1 = 116)
      PARAMETER (CUTOF2 = 226)

C========================= END DIOASET ===========================
C     INCLUDE 'INSTCO.INC'
C=========================== DIINSTCO =============================

C $   (JR)
C $   THIS INCLUDE FILE IS PART OF THE SUPPLIED GVAR NAV SOFTWARE
C $
C $   DESCRIPTION:
C $
C $
C $$  DIINSTCO = INCLUDE

      INTEGER*4 INCMAX(2)
C                       NUMBER OF INCREMENTS PER CYCLE
      REAL*4 ELVMAX(2)
C                       BOUNDS IN ELEVATION (RADIANS)
      REAL*4 SCNMAX(2)
C                       BOUNDS IN SCAN ANGLE (RADIANS)
      REAL*4 ELVINC(2)
C                       CHANGE IN ELEVATION ANGLE PER INCREMENT (RAD)
      REAL*4 SCNINC(2)
C                       CHANGE IN SCAN ANGLE PER INCREMENT (RADIANS)
      REAL*4 ELVLN(2)
C                       ELEVATION ANGLE PER DETECTOR LINE (RADIANS)
      REAL*4 SCNPX(2)
C                       SCAN ANGLE PER PIXEL (RADIANS)

      COMMON /INSTCO/ INCMAX,ELVMAX,SCNMAX,
     1   ELVINC,SCNINC,ELVLN,SCNPX
C========================= END DIINSTCO ===========================

      COMMON /GVRCOM/ ITYPE, INSTR , SUBLAT, SUBLON

C--------------------------------------------------------------

C TRANSFORM LINE/PIXEL TO GEOGRAPHIC COORDINATES:
C SET INPUT LINE/PIXEL NUMBERS
      RL = XLIN
      RP = XELE

C IF DOING SOUNDER NAV, HAVE TO TRICK ROUTINES INTO THINKING IMAGE IS
C AT RES 1, BECAUSE NAV ROUTINES TAKE SOUNDER RES INTO ACCOUNT
      IF (INSTR.EQ.2) THEN
        RL = (RL+9.)/10.
        RP = (RP+9.)/10.
      ENDIF

C COMPUTE ELEVATION AND SCAN ANGLES (E,S) RELATED TO INPUT
C LINE AND PIXEL NUMBERS
      E = EVLN(INSTR,RL)
      S = SCPX(INSTR,RP)

C TRANSFORM ELEVATION AND SCAN ANGLES TO GEOGRAPHIC COORDINATES
      CALL LPOINT(E,S,TMPLAT,TMPLON,STAT)
      IF (STAT.NE.0) GOTO 900
      TMPLAT = TMPLAT * DEG
      TMPLON = TMPLON * DEG

C PUT LONGITUDE INTO MCIDAS FORM
      TMPLON = -TMPLON

C SEE IF WE HAVE TO CONVERT TO X Y Z COORDINATES
      IF (ITYPE.EQ.2) THEN
        YLAT = TMPLAT
        YLON = TMPLON
        CALL LLCART(YLAT,YLON,XLAT,XLON,Z)
      ELSE
	XLAT = TMPLAT
	XLON = TMPLON
      ENDIF

      GVRSAE=0
      RETURN

 900  CONTINUE
      XLAT = TMPLAT
      XLON = -TMPLON
      GVRSAE=-1
      RETURN
      END

      INTEGER FUNCTION GVREAS(ZLAT,ZLON,Z,XLIN,XELE,XDUM)

      IMPLICIT NONE

      REAL         ZLAT  , ZLON  , Z     , XLIN  , XELE  , XDUM
      REAL         X     , Y
      REAL*8       TMPLIN, TMPELE, E     , S     , TMPLAT, TMPLON
      REAL*8       SUBLAT, SUBLON
      INTEGER      INSTR , ITYPE , IER

C     INCLUDE 'ELCONS.INC'
C=========================== DIELCONS =============================

C $   (JR)
C $   THIS INCLUDE FILE IS PART OF THE SUPPLIED GVAR NAV SOFTWARE
C $
C $   DESCRIPTION:
C $
C $
C $$  DIELCONS = INCLUDE

      REAL*8 PI
           PARAMETER (PI=3.141592653589793D0)
      REAL*8 DEG
           PARAMETER (DEG=180.D0/PI)
      REAL*8 RAD
           PARAMETER (RAD=PI/180.D0)
C                    DEGREES TO RADIANS CONVERSION PI/180
      REAL*8 NOMORB
           PARAMETER (NOMORB=42164.365D0)
C                    NOMINAL RADIAL DISTANCE OF SATELLITE (KM)
      REAL*8 AE
           PARAMETER (AE=6378.137D0)
C                    EARTH EQUATORIAL RADIUS (KM)
      REAL*8 FER
           PARAMETER (FER=1.D0-(6356.7533D0/AE))
C                    EARTH FLATTENING COEFFICIENT = 1-(BE/AE)
      REAL*4 AEBE2
           PARAMETER (AEBE2=1.D0/(1.D0-FER)**2)
      REAL*4 AEBE3
           PARAMETER (AEBE3=AEBE2-1.)
      REAL*4 AEBE4
           PARAMETER (AEBE4=(1.D0-FER)**4-1.)
C========================= END DIELCONS ===========================
C     INCLUDE 'ELCOMM.INC'
C=========================== DIELCOMM =============================

C $   (JR)
C $   THIS INCLUDE FILE IS PART OF THE SUPPLIED GVAR NAV SOFTWARE
C $
C $   DESCRIPTION:
C $
C $
C $$  DIELCOMM = INCLUDE

C     COMMON VARIABLES
C
      REAL*8 XS(3)
C                      NORMALIZED S/C POSITION IN ECEF COORDINATES
      REAL*8 BT(3,3)
C                      ECEF TO INSTRUMENT COORDINATES TRANSFORMATION
      REAL*8  Q3
C                      USED IN SUBRTN LPOINT
      REAL*8 PITCH,ROLL,YAW
C                          PITCH,ROLL,YAW ANGLES OF INSTRUMENT (RAD)
      REAL*4 PMA,RMA
C                          PITCH,ROLL MISALIGNMENTS OF INSTRUMENT (RAD)
         COMMON /ELCOMM/ XS,BT,Q3,PITCH,ROLL,YAW,PMA,RMA
C========================= END DIELCOMM ===========================
C     INCLUDE 'OASET.INC'
C=========================== DIOASET =============================

C $   (JR)
C $   THIS INCLUDE FILE IS PART OF THE GVAR NAVIGATION SUBSYSTEM.
C $
C $   DESCRIPTION:
C $      THESE ARE CONSTANTS USED TO DESCRIBE THE GVAR O&A SET STRUCTURE
C $
C $$  DIOASET = INCLUDE, NAVIGATION, GVAR

C DEFINITIONS OF PARAMETERS
C OASIZE   : SIZE OF GVAR ORBIT AND ATTITUDE SET
C PCOEFS   : START OF PITCH COEFFICIENTS
C RMACFS   : START OF RMA COEFFICIENTS
C CUTOF1   : FIRST DIVIDING POINT IN O&A SET (WHERE DIVIDED INTO 128
C            WORD SETS)
C CUTOF2   : SECOND DIVIDING POINT IN O&A SET (WHERE DIVIDED INTO 128
C            WORD SETS)

      INTEGER      OASIZE, PCOEFS, RMACFS, CUTOF1, CUTOF2

      PARAMETER (OASIZE = 336)
      PARAMETER (PCOEFS = 117)
      PARAMETER (RMACFS = 227)
      PARAMETER (CUTOF1 = 116)
      PARAMETER (CUTOF2 = 226)

C========================= END DIOASET ===========================
C     INCLUDE 'INSTCO.INC'
C=========================== DIINSTCO =============================

C $   (JR)
C $   THIS INCLUDE FILE IS PART OF THE SUPPLIED GVAR NAV SOFTWARE
C $
C $   DESCRIPTION:
C $
C $
C $$  DIINSTCO = INCLUDE

      INTEGER*4 INCMAX(2)
C                       NUMBER OF INCREMENTS PER CYCLE
      REAL*4 ELVMAX(2)
C                       BOUNDS IN ELEVATION (RADIANS)
      REAL*4 SCNMAX(2)
C                       BOUNDS IN SCAN ANGLE (RADIANS)
      REAL*4 ELVINC(2)
C                       CHANGE IN ELEVATION ANGLE PER INCREMENT (RAD)
      REAL*4 SCNINC(2)
C                       CHANGE IN SCAN ANGLE PER INCREMENT (RADIANS)
      REAL*4 ELVLN(2)
C                       ELEVATION ANGLE PER DETECTOR LINE (RADIANS)
      REAL*4 SCNPX(2)
C                       SCAN ANGLE PER PIXEL (RADIANS)

      COMMON /INSTCO/ INCMAX,ELVMAX,SCNMAX,
     1   ELVINC,SCNINC,ELVLN,SCNPX
C========================= END DIINSTCO ===========================

      COMMON /GVRCOM/ ITYPE, INSTR, SUBLAT, SUBLON

C--------------------------------------------------------------

      GVREAS = 0

C IF IN CARTESIAN COORDS, TRANSFORM TO LAT/LON
      IF (ITYPE.EQ.2) THEN
        X = ZLAT
        Y = ZLON
        CALL CARTLL(X,Y,Z,ZLAT,ZLON)
      ENDIF

      IF (ABS(ZLAT).GT.90.) THEN
        GVREAS = -1
        RETURN
      ENDIF

C ELUG ROUTINES USE REAL*8 VARS; USE TEMP VARIABLES
      TMPLAT = ZLAT
      TMPLON = ZLON

C TRANSFORM LAT/LON TO ELEVATION AND SCAN ANGLES
      CALL GPOINT(TMPLAT*RAD,(-TMPLON)*RAD,E,S,IER)

C SEE IF THE LAT/LON GIVEN IS ABLE TO BE SEEN
      IF (IER.NE.0) THEN
        GVREAS = -1
        RETURN
      ENDIF

C CONVERT ELEVATION AND SCAN ANGLES TO LINE/PIXEL COORDINATES
      CALL EVSC2L(INSTR,E,S,TMPLIN,TMPELE)

C CONVERT INTERNAL 8 BYTE VALUES TO 4 BYTES
      XLIN = TMPLIN
      XELE = TMPELE

C IF DOING SOUNDER NAV, CHANGE LIN & ELE RETURNED TO RES 10 VALUES
      IF (INSTR.EQ.2) THEN
        XLIN = XLIN*10.-9.
        XELE = XELE*10.-9.
      ENDIF

      RETURN
      END

      INTEGER FUNCTION GVROPT(IFUNC,XIN,XOUT)

      IMPLICIT NONE

C     INCLUDE 'ELCONS.INC'
C=========================== DIELCONS =============================

C $   (JR)
C $   THIS INCLUDE FILE IS PART OF THE SUPPLIED GVAR NAV SOFTWARE
C $
C $   DESCRIPTION:
C $
C $
C $$  DIELCONS = INCLUDE

      REAL*8 PI
           PARAMETER (PI=3.141592653589793D0)
      REAL*8 DEG
           PARAMETER (DEG=180.D0/PI)
      REAL*8 RAD
           PARAMETER (RAD=PI/180.D0)
C                    DEGREES TO RADIANS CONVERSION PI/180
      REAL*8 NOMORB
           PARAMETER (NOMORB=42164.365D0)
C                    NOMINAL RADIAL DISTANCE OF SATELLITE (KM)
      REAL*8 AE
           PARAMETER (AE=6378.137D0)
C                    EARTH EQUATORIAL RADIUS (KM)
      REAL*8 FER
           PARAMETER (FER=1.D0-(6356.7533D0/AE))
C                    EARTH FLATTENING COEFFICIENT = 1-(BE/AE)
      REAL*4 AEBE2
           PARAMETER (AEBE2=1.D0/(1.D0-FER)**2)
      REAL*4 AEBE3
           PARAMETER (AEBE3=AEBE2-1.)
      REAL*4 AEBE4
           PARAMETER (AEBE4=(1.D0-FER)**4-1.)
C========================= END DIELCONS ===========================
C     INCLUDE 'ELCOMM.INC'
C=========================== DIELCOMM =============================

C $   (JR)
C $   THIS INCLUDE FILE IS PART OF THE SUPPLIED GVAR NAV SOFTWARE
C $
C $   DESCRIPTION:
C $
C $
C $$  DIELCOMM = INCLUDE

C     COMMON VARIABLES
C
      REAL*8 XS(3)
C                      NORMALIZED S/C POSITION IN ECEF COORDINATES
      REAL*8 BT(3,3)
C                      ECEF TO INSTRUMENT COORDINATES TRANSFORMATION
      REAL*8  Q3
C                      USED IN SUBRTN LPOINT
      REAL*8 PITCH,ROLL,YAW
C                          PITCH,ROLL,YAW ANGLES OF INSTRUMENT (RAD)
      REAL*4 PMA,RMA
C                          PITCH,ROLL MISALIGNMENTS OF INSTRUMENT (RAD)
         COMMON /ELCOMM/ XS,BT,Q3,PITCH,ROLL,YAW,PMA,RMA
C========================= END DIELCOMM ===========================

      INTEGER      IFUNC , LIT   , ITYPE , INSTR , JDAY  , JTIME
      INTEGER      IROUND, ITIMC, LASDAY, LASTIM
      REAL         FLAT  , FLON  , GHA   , XIN(*)
      REAL        DEC   , XLAT  , XLON  , XOUT(*)
      REAL*8 AEC , FERC  , AEBE2C, AEBE3C, AEBE4C, B(3,3)
      REAL*8 SUBLAT, SUBLON, R     , DR    , PHI

      COMMON /GVRCOM/ ITYPE , INSTR , SUBLAT, SUBLON
      COMMON /RADCOM/ AEC   , FERC  , AEBE2C, AEBE3C, AEBE4C
      COMMON /SAVCOM/ B     , DR    , PHI

      DATA LASDAY/-1/,LASTIM/-1/

C-----------------------------------------------------------------

      GVROPT = -1

      IF (IFUNC.EQ.LIT('SPOS')) THEN
        XOUT(1) = SUBLAT * DEG
        XOUT(2) = (-SUBLON) * DEG
        XOUT(2) = AMOD(XOUT(2),360.)
        GVROPT = 0

      ELSEIF (IFUNC.EQ.LIT('HGT ')) THEN
        AEC    = AE + DBLE(XIN(1))
        FERC   = 1.D0 - ((6356.7533D0 + DBLE(XIN(1))) / AEC)
        AEBE2C = 1.D0 / ((1.D0 - FERC)*(1.D0 - FERC))
        AEBE3C = AEBE2C - 1.D0
        AEBE4C = (1.D0-FERC)*(1.D0-FERC)*(1.D0-FERC)*(1.D0-FERC) - 1.D0

C RECOMPUTE VALUES THAT DEPEND ON EARTH RADIUS
        SUBLAT = DATAN(AEBE2C*DTAN(PHI))
        R      = (NOMORB+DR)/AEC
        XS(1)  = (-B(1,3))*R
        XS(2)  = (-B(2,3))*R
        XS(3)  = (-B(3,3))*R
        Q3     = XS(1)**2 + XS(2)**2 + AEBE2C * XS(3)**2 - 1.0
        GVROPT = 0

      ELSEIF (IFUNC.EQ.LIT('ANG ')) THEN
        JDAY  = IROUND(XIN(1))
        JTIME = ITIMC(XIN(2))
        FLAT  = XIN(3)
        FLON  = XIN(4)
        IF (JDAY.NE.LASDAY.OR.JTIME.NE.LASTIM) THEN
          CALL SOLARP(JDAY,JTIME,GHA,DEC,XLAT,XLON)
          LASDAY = JDAY
          LASTIM = JTIME
        ENDIF
        CALL GVRANG(JDAY, JTIME, FLAT, FLON, GHA, DEC, XOUT(1), XOUT(2),
     &                                                          XOUT(3))
        GVROPT = 0
      ENDIF

      RETURN
      END

C THE FOLLOWING ARE SUBSIDIARY ROUTINES

      SUBROUTINE SETCON(INSTR, NADNSC, NADNSI, NADEWC, NADEWI)

      IMPLICIT NONE

      INTEGER INSTR, NADNSC, NADNSI, NADEWC, NADEWI
C      character*12 cfg

C     INCLUDE 'ELCONS.INC'
C=========================== DIELCONS =============================

C $   (JR)
C $   THIS INCLUDE FILE IS PART OF THE SUPPLIED GVAR NAV SOFTWARE
C $
C $   DESCRIPTION:
C $
C $
C $$  DIELCONS = INCLUDE

      REAL*8 PI
           PARAMETER (PI=3.141592653589793D0)
      REAL*8 DEG
           PARAMETER (DEG=180.D0/PI)
      REAL*8 RAD
           PARAMETER (RAD=PI/180.D0)
C                    DEGREES TO RADIANS CONVERSION PI/180
      REAL*8 NOMORB
           PARAMETER (NOMORB=42164.365D0)
C                    NOMINAL RADIAL DISTANCE OF SATELLITE (KM)
      REAL*8 AE
           PARAMETER (AE=6378.137D0)
C                    EARTH EQUATORIAL RADIUS (KM)
      REAL*8 FER
           PARAMETER (FER=1.D0-(6356.7533D0/AE))
C                    EARTH FLATTENING COEFFICIENT = 1-(BE/AE)
      REAL*4 AEBE2
           PARAMETER (AEBE2=1.D0/(1.D0-FER)**2)
      REAL*4 AEBE3
           PARAMETER (AEBE3=AEBE2-1.)
      REAL*4 AEBE4
           PARAMETER (AEBE4=(1.D0-FER)**4-1.)
C========================= END DIELCONS ===========================
C     INCLUDE 'INSTCO.INC'
C=========================== DIINSTCO =============================

C $   (JR)
C $   THIS INCLUDE FILE IS PART OF THE SUPPLIED GVAR NAV SOFTWARE
C $
C $   DESCRIPTION:
C $
C $
C $$  DIINSTCO = INCLUDE

      INTEGER*4 INCMAX(2)
C                       NUMBER OF INCREMENTS PER CYCLE
      REAL*4 ELVMAX(2)
C                       BOUNDS IN ELEVATION (RADIANS)
      REAL*4 SCNMAX(2)
C                       BOUNDS IN SCAN ANGLE (RADIANS)
      REAL*4 ELVINC(2)
C                       CHANGE IN ELEVATION ANGLE PER INCREMENT (RAD)
      REAL*4 SCNINC(2)
C                       CHANGE IN SCAN ANGLE PER INCREMENT (RADIANS)
      REAL*4 ELVLN(2)
C                       ELEVATION ANGLE PER DETECTOR LINE (RADIANS)
      REAL*4 SCNPX(2)
C                       SCAN ANGLE PER PIXEL (RADIANS)

      COMMON /INSTCO/ INCMAX,ELVMAX,SCNMAX,
     1   ELVINC,SCNINC,ELVLN,SCNPX
C========================= END DIINSTCO ===========================

      INCMAX(1)  = 6136
      INCMAX(2)  = 2805
      ELVINC(1)  = 8.E-6
      ELVINC(2)  = 17.5E-6
      SCNINC(1)  = 16.E-6
      SCNINC(2)  = 35.E-6
      ELVLN(1)   = 28.E-6
      ELVLN(2)   = 280.E-6
      SCNPX(1)   = 16.E-6
      SCNPX(2)   = 280.E-6
      ELVMAX(1)  = 0.220896
      ELVMAX(2)  = 0.22089375
      SCNMAX(1)  = 0.24544
      SCNMAX(2)  = 0.2454375

C new code because of change to elug; nadir position is available
C in the signal, so should compute 4 values above using them
C instead of having them hardwired
C
C*  These ddest() calls appear to just output some text - don't need them
C*  -Cowie
C
C      call ddest('nadnsc: ', nadnsc)
C      call ddest('nadnsi: ', nadnsi)
C      call ddest('nadewc: ', nadewc)
C      call ddest('nadewi: ', nadewi)
c
c--- new code from Kathy Kelly for sounder nav - 10/27
c
      if (nadnsc .ne. 0 .and. nadnsi .ne. 0 .and. nadewc .ne. 0
     &                                 .and. nadewi .ne. 0) then
	IF (INSTR .EQ. 1) THEN

           ELVMAX(INSTR) = (NADNSC*INCMAX(INSTR)+NADNSI)
     &                                             *ELVINC(INSTR)

        ELSE

           ELVMAX(INSTR) = ((9-NADNSC)*INCMAX(INSTR)-NADNSI)
     &                                             *ELVINC(INSTR)

	ENDIF

        SCNMAX(INSTR) = (NADEWC*INCMAX(INSTR)+NADEWI)
     &                                            *SCNINC(INSTR)
      endif
c
c--- end of new code from Kathy Kelly for sounder nav - 10/27
c

C      call ddest('elvmax(1): '//cfg(elvmax(1)), 0)
C      call ddest('elvmax(2): '//cfg(elvmax(2)), 0)
C      call ddest('scnmax(1): '//cfg(scnmax(1)), 0)
C      call ddest('scnmax(2): '//cfg(scnmax(2)), 0)

      RETURN
      END

      FUNCTION TIMEX(NY,ND,NH,NM,S)

      IMPLICIT NONE

      INTEGER      ND    , NH    , NM    , J     , NY
      REAL*8       TIMEX , S

C-----------------------------------------------------------------------

      J = ND + 1461 * (NY + 4799) / 4 - 3 * ((NY + 4899) / 100) / 4
     &                                                        - 2465022

C COMPUTE ACTUAL TIME IN MINUTES FROM JANUARY 1.0, 1950
      TIMEX = J * 1440.D0 + NH * 60.D0 + NM + S / 60.D0

      RETURN
      END

      FUNCTION TIME50(I)

      IMPLICIT NONE

      INTEGER      NY    , ND    , NH    , NM    , J     , I(2)
      INTEGER      IAA   , IAB   , NBC   , IAC   , DEF
      REAL*8       TIME50, S

C-----------------------------------------------------------------

      NY  = I(1) / 10000
      IAA = I(1) - (NY*10000)

      ND  = (I(1) - (NY*10000)) * .1
      IAB = (IAA - (ND * 10)) * 10
      NBC = I(2) / 10000000.
      IAC = I(2) - (NBC * 10000000)
      NH  = IAB + NBC
      DEF = I(2) - IAC
      NM  = IAC * .00001
      S   = (I(2) - (DEF + (NM * 100000))) * .001

C HERE WE CONVERT INTEGER YEAR AND DAY OF YEAR TO NUMBER OF
C DAYS FROM 0 HOUR UT, 1950 JAN. 1.0
C THIS CONVERTION IS BASED ON AN ALGORITHM BY FLIEGEL AND VAN
C FLANDERN, COMM. OF ACM, VOL.11, NO. 10, OCT. 1968 (P.657)

      J = ND + 1461 * (NY + 4799) / 4 - 3 * ((NY + 4899) / 100) / 4
     &                                                       - 2465022

C COMPUTE TIME IN MINUTES FROM JANUARY 1.0, 1950
      TIME50 = J * 1440.D0 + NH * 60.D0 + NM + S / 60.D0

      RETURN
      END

      REAL*8 FUNCTION GATT(K0,RPARMS,IPARMS,WA,TE)

      IMPLICIT NONE

C     INCLUDE 'MXCDSZ.INC'
C=========================== DIMXCDSZ =============================

C $   (JR)
C $   THIS INCLUDE FILE IS PART OF THE NAVIGATION SUBSYSTEM
C $
C $   DESCRIPTION:
C $      THIS IS THE VALUE OF LARGEST NAVIGATION CODICIL SIZE ALLOWED
C $
C $$  DIMXCDSZ = INCLUDE, NAVIGATION

C DEFINITIONS OF PARAMETERS
C MXCDSZ   : LARGEST ALLOWED CODICIL SIZE

      INTEGER      MXCDSZ

      PARAMETER (MXCDSZ = 5*128)

C========================= END DIMXCDSZ ===========================

C CALLING PARAMETERS
      INTEGER K0
C STARTING POSITION OF A PARAMETER SUBSET IN THE REAL O&A SET
      REAL*4 RPARMS(MXCDSZ)
C INPUT O&A PARAMETER SET
      INTEGER IPARMS(MXCDSZ)
C INPUT O&A PARAMETER SET
      REAL*8 WA
C INPUT SOLAR ORBIT ANGLE IN RADIANS
      REAL*8 TE
C INPUT EXPONENTIAL TIME DELAY FROM EPOCH (MINUTES)

C LOCAL VARIABLES
      INTEGER      KKK   , I     , J     , M     , L     , LL    , K
      REAL*8       IR    , JR    , MR    , ATT

C------------------------------------------------------------------

C CONSTANT COMPONENT
      K   = K0
      ATT = RPARMS(K+2)

C COMPUTES THE EXPONENTIAL TERM
      IF (TE.GE.0) ATT = ATT + RPARMS(K) * DEXP((-TE) / RPARMS(K+1))

C EXTRACTS THE NUMBER OF SINUSOIDS
      IR = REAL(IPARMS(K+3))
      I  = NINT(IR)

C CALCULATION OF SINUSOIDS
      DO 10 L=1,I
        ATT = ATT + RPARMS(K+2*L+2) * DCOS(WA*L+RPARMS(K+2*L+3))
 10   CONTINUE

C POINTER TO THE NUMBER OF MONOMIAL SINUSOIDS
      K = K + 34

C EXTACTS NUMBER OF MONOMIAL SINUSOIDS
      IR  = REAL(IPARMS(K))
      KKK = IPARMS(K)

C COMPUTES MONOMIAL SINUSOIDS
      DO 20 L=1,KKK
        LL = K + 5 * L

C ORDER OF SINUSOID
        JR = REAL(IPARMS(LL-4))

C ORDER OF MONOMIAL SINUSOID
        MR  = REAL(IPARMS(LL-3))
        J   = NINT(JR)
        M   = NINT(MR)
        ATT = ATT + RPARMS(LL-2) * ((WA - RPARMS(LL))**MR)
     &                                      * DCOS(JR*WA+RPARMS(LL-1))
 20   CONTINUE

      GATT = ATT
      RETURN
      END

      SUBROUTINE EVSC2L(INSTR,ELEV,SCAN,RL,RP)

      IMPLICIT NONE

C CALLING PARAMETERS
      INTEGER INSTR
C INSTRUMENT CODE (1-IMAGER, 2-SOUNDER)
      REAL*8 ELEV
C ELEVATION ANGLE IN RADIANS
      REAL*8 SCAN
C SCAN ANGLE IN RADIANS
      REAL*8 RL
C LINE NUMBER
      REAL*8 RP
C PIXEL NUMBER

C     INCLUDE 'INSTCO.INC'
C=========================== DIINSTCO =============================

C $   (JR)
C $   THIS INCLUDE FILE IS PART OF THE SUPPLIED GVAR NAV SOFTWARE
C $
C $   DESCRIPTION:
C $
C $
C $$  DIINSTCO = INCLUDE

      INTEGER*4 INCMAX(2)
C                       NUMBER OF INCREMENTS PER CYCLE
      REAL*4 ELVMAX(2)
C                       BOUNDS IN ELEVATION (RADIANS)
      REAL*4 SCNMAX(2)
C                       BOUNDS IN SCAN ANGLE (RADIANS)
      REAL*4 ELVINC(2)
C                       CHANGE IN ELEVATION ANGLE PER INCREMENT (RAD)
      REAL*4 SCNINC(2)
C                       CHANGE IN SCAN ANGLE PER INCREMENT (RADIANS)
      REAL*4 ELVLN(2)
C                       ELEVATION ANGLE PER DETECTOR LINE (RADIANS)
      REAL*4 SCNPX(2)
C                       SCAN ANGLE PER PIXEL (RADIANS)

      COMMON /INSTCO/ INCMAX,ELVMAX,SCNMAX,
     1   ELVINC,SCNINC,ELVLN,SCNPX
C========================= END DIINSTCO ===========================

C--------------------------------------------------------------

C COMPUTE FRACTIONAL LINE NUMBER
      RL = (ELVMAX(INSTR) - ELEV) / ELVLN(INSTR)
      IF (INSTR.EQ.1) THEN
        RL = RL + 4.5
      ELSE
        RL = RL + 2.5
      ENDIF

C COMPUTE FRACTIONAL PIXEL NUMBER
      RP = (SCNMAX(INSTR) + SCAN) / SCNPX(INSTR) + 1.

      RETURN
      END

      REAL*8 FUNCTION EVLN(INSTR,RLINE)

      IMPLICIT NONE

C CALLING PARAMETERS
      INTEGER INSTR
C INSTRUMENT CODE (1-IMAGER, 2-SOUNDER)
      REAL*8  RLINE
C FRACTIONAL LINE  NUMBER

C     INCLUDE 'INSTCO.INC'
C=========================== DIINSTCO =============================

C $   (JR)
C $   THIS INCLUDE FILE IS PART OF THE SUPPLIED GVAR NAV SOFTWARE
C $
C $   DESCRIPTION:
C $
C $
C $$  DIINSTCO = INCLUDE

      INTEGER*4 INCMAX(2)
C                       NUMBER OF INCREMENTS PER CYCLE
      REAL*4 ELVMAX(2)
C                       BOUNDS IN ELEVATION (RADIANS)
      REAL*4 SCNMAX(2)
C                       BOUNDS IN SCAN ANGLE (RADIANS)
      REAL*4 ELVINC(2)
C                       CHANGE IN ELEVATION ANGLE PER INCREMENT (RAD)
      REAL*4 SCNINC(2)
C                       CHANGE IN SCAN ANGLE PER INCREMENT (RADIANS)
      REAL*4 ELVLN(2)
C                       ELEVATION ANGLE PER DETECTOR LINE (RADIANS)
      REAL*4 SCNPX(2)
C                       SCAN ANGLE PER PIXEL (RADIANS)

      COMMON /INSTCO/ INCMAX,ELVMAX,SCNMAX,
     1   ELVINC,SCNINC,ELVLN,SCNPX
C========================= END DIINSTCO ===========================

C----------------------------------------------------------------

      IF (INSTR.EQ.1) THEN
        EVLN = ELVMAX(INSTR) - (RLINE - 4.5) * ELVLN(INSTR)
      ELSE
        EVLN = ELVMAX(INSTR) - (RLINE - 2.5) * ELVLN(INSTR)
      ENDIF

      RETURN
      END

      REAL*8 FUNCTION SCPX(INSTR,PIX)

      IMPLICIT NONE

C CALLING PARAMETERS
      INTEGER INSTR
C INSTRUMENT CODE (1-IMAGER, 2-SOUNDER)
      REAL*8 PIX
C FRACTIONAL PIXEL NUMBER

C     INCLUDE 'INSTCO.INC'
C=========================== DIINSTCO =============================

C $   (JR)
C $   THIS INCLUDE FILE IS PART OF THE SUPPLIED GVAR NAV SOFTWARE
C $
C $   DESCRIPTION:
C $
C $
C $$  DIINSTCO = INCLUDE

      INTEGER*4 INCMAX(2)
C                       NUMBER OF INCREMENTS PER CYCLE
      REAL*4 ELVMAX(2)
C                       BOUNDS IN ELEVATION (RADIANS)
      REAL*4 SCNMAX(2)
C                       BOUNDS IN SCAN ANGLE (RADIANS)
      REAL*4 ELVINC(2)
C                       CHANGE IN ELEVATION ANGLE PER INCREMENT (RAD)
      REAL*4 SCNINC(2)
C                       CHANGE IN SCAN ANGLE PER INCREMENT (RADIANS)
      REAL*4 ELVLN(2)
C                       ELEVATION ANGLE PER DETECTOR LINE (RADIANS)
      REAL*4 SCNPX(2)
C                       SCAN ANGLE PER PIXEL (RADIANS)

      COMMON /INSTCO/ INCMAX,ELVMAX,SCNMAX,
     1   ELVINC,SCNINC,ELVLN,SCNPX
C========================= END DIINSTCO ===========================

C--------------------------------------------------------------------

      SCPX = (PIX - 1.) * SCNPX(INSTR) - SCNMAX(INSTR)

      RETURN
      END

      SUBROUTINE INST2E(R,P,Y,A,AT)

      IMPLICIT NONE

C CALLING PARAMETERS
      REAL*8 R
C ROLL ANGLE IN RADIANS
      REAL*8 P
C PITCH ANGLE IN RADIANS
      REAL*8 Y
C YAW ANGLE IN RADIANS
      REAL*8 A(3,3)
C SPACECRAFT TO ECEF COORDINATES TRANSFORMATION MATRIX
      REAL*8 AT(3,3)
C INSTRUMENT TO ECEF COORDINATES TRANSFORMATION MATRIX

      REAL*8       RPY(3,3)
      INTEGER*4    I     , J

C-----------------------------------------------------------------------

C WE COMPUTE INSTRUMENT TO BODY COORDINATES TRANSFORMATION
C MATRIX BY USING A SMALL ANGLE APPROXIMATION OF TRIGONOMETRIC
C FUNKTIONS OF THE ROLL, PITCH AND YAW.
      RPY(1,1) = 1. - 0.5 * (P * P + Y * Y)
      RPY(1,2) = -Y
      RPY(1,3) = P
      RPY(2,1) = Y + P * R
      RPY(2,2) = 1. - 0.5 * (Y * Y + R * R)
      RPY(2,3) = -R
      RPY(3,1) = -P + R * Y
      RPY(3,2) = R + P * Y
      RPY(3,3) = 1. - 0.5 * (P * P + R * R)

C MULTIPLICATION OF MATRICES A AND RPY
      DO 20 I=1,3
        DO 10 J=1,3
          AT(I,J) = A(I,1) * RPY(1,J) + A(I,2) * RPY(2,J) + A(I,3)
     &                                                    * RPY(3,J)
 10     CONTINUE
 20   CONTINUE

      RETURN
      END

      SUBROUTINE LPOINT(ALPHA,ZETA,RLAT,RLON,IERR)

      IMPLICIT NONE

C CALLING PARAMETERS
      REAL*8   ALPHA
C ELEVATION ANGLE (RAD)
      REAL*8   ZETA
C SCAN ANGLE (RAD)
      REAL*8   RLAT
C LATITUDE IN RADIANS (OUTPUT)
      REAL*8   RLON
C LONGITUDE IN RADIANS (OUTPUT)
      INTEGER IERR
C OUTPUT STATUS; 0 - POINT ON THE EARTH FOUND,
C                1 - INSTRUMENT POINTS OFF EARTH

      REAL*8       AEC   , FERC  , Q1    , Q2    , D     , H     , G1(3)
      REAL*8       SA    , CA    , DA    , DZ    , D1    , CZ    , G(3)
      REAL*8       U(3)
      REAL*8       AEBE2C, AEBE3C, AEBE4C

C     INCLUDE 'ELCONS.INC'
C=========================== DIELCONS =============================

C $   (JR)
C $   THIS INCLUDE FILE IS PART OF THE SUPPLIED GVAR NAV SOFTWARE
C $
C $   DESCRIPTION:
C $
C $
C $$  DIELCONS = INCLUDE

      REAL*8 PI
           PARAMETER (PI=3.141592653589793D0)
      REAL*8 DEG
           PARAMETER (DEG=180.D0/PI)
      REAL*8 RAD
           PARAMETER (RAD=PI/180.D0)
C                    DEGREES TO RADIANS CONVERSION PI/180
      REAL*8 NOMORB
           PARAMETER (NOMORB=42164.365D0)
C                    NOMINAL RADIAL DISTANCE OF SATELLITE (KM)
      REAL*8 AE
           PARAMETER (AE=6378.137D0)
C                    EARTH EQUATORIAL RADIUS (KM)
      REAL*8 FER
           PARAMETER (FER=1.D0-(6356.7533D0/AE))
C                    EARTH FLATTENING COEFFICIENT = 1-(BE/AE)
      REAL*4 AEBE2
           PARAMETER (AEBE2=1.D0/(1.D0-FER)**2)
      REAL*4 AEBE3
           PARAMETER (AEBE3=AEBE2-1.)
      REAL*4 AEBE4
           PARAMETER (AEBE4=(1.D0-FER)**4-1.)
C========================= END DIELCONS ===========================
C     INCLUDE 'ELCOMM.INC'
C=========================== DIELCOMM =============================

C $   (JR)
C $   THIS INCLUDE FILE IS PART OF THE SUPPLIED GVAR NAV SOFTWARE
C $
C $   DESCRIPTION:
C $
C $
C $$  DIELCOMM = INCLUDE

C     COMMON VARIABLES
C
      REAL*8 XS(3)
C                      NORMALIZED S/C POSITION IN ECEF COORDINATES
      REAL*8 BT(3,3)
C                      ECEF TO INSTRUMENT COORDINATES TRANSFORMATION
      REAL*8  Q3
C                      USED IN SUBRTN LPOINT
      REAL*8 PITCH,ROLL,YAW
C                          PITCH,ROLL,YAW ANGLES OF INSTRUMENT (RAD)
      REAL*4 PMA,RMA
C                          PITCH,ROLL MISALIGNMENTS OF INSTRUMENT (RAD)
         COMMON /ELCOMM/ XS,BT,Q3,PITCH,ROLL,YAW,PMA,RMA
C========================= END DIELCOMM ===========================

      COMMON /RADCOM/ AEC   , FERC  , AEBE2C, AEBE3C, AEBE4C

C-----------------------------------------------------------------------

      IERR=1

C COMPUTES TRIGONOMETRIC FUNKTIONS OF THE SCAN AND ELEVATION
C ANGLES CORRECTED FOR THE ROLL AND PITCH MISALIGNMENTS
      CA = DCOS(ALPHA)
      SA = DSIN(ALPHA)
      CZ = DCOS(ZETA)
      DA = ALPHA-PMA*SA*(1.0D0/CZ+DTAN(ZETA))-RMA*(1.0D0-CA/CZ)
      DZ = ZETA + RMA * SA

C CORRECTED SCAN ANGLE
      CZ = DCOS(DZ)

C COMPUTES POINTING VECTOR IN INSTRUMENT COORDINATES
      G(1) = DSIN(DZ)
      G(2) = (-CZ) * DSIN(DA)
      G(3) = CZ * DCOS(DA)

C TRANSFORMS THE POINTING VECTOR TO EARTH FIXED COORDINATES
      G1(1) = BT(1,1) * G(1) + BT(1,2) * G(2) + BT(1,3) * G(3)
      G1(2) = BT(2,1) * G(1) + BT(2,2) * G(2) + BT(2,3) * G(3)
      G1(3) = BT(3,1) * G(1) + BT(3,2) * G(2) + BT(3,3) * G(3)

C COMPUTES COEFFICIENTS AND SOLVES A QUADRATIC EQUATION TO
C FIND THE INTERSECT OF THE POINTING VECTOR WITH THE EARTH
C SURFACE
      Q1 = G1(1)**2 + G1(2)**2 + AEBE2C * G1(3)**2
      Q2 = XS(1) * G1(1) + XS(2) * G1(2) + AEBE2C * XS(3) * G1(3)
      D  = Q2 * Q2 - Q1 * Q3
      IF (DABS(D).LT.1.D-9) D=0.

C IF THE DISCRIMINANT OF THE EQUATION, D, IS NEGATIVE, THE
C INSTRUMENT POINTS OFF THE EARTH
      IF (D.LT.0) THEN
        RLAT = 999999.
        RLON = 999999.
        RETURN
      ENDIF
      D = DSQRT(D)

C SLANT DISTANCE FROM THE SATELLITE TO THE EARTH POINT
      H = (-(Q2 + D)) / Q1

C CARTESIAN COORDINATES OF THE EARTH POINT
      U(1) = XS(1) + H * G1(1)
      U(2) = XS(2) + H * G1(2)
      U(3) = XS(3) + H * G1(3)

C SINUS OF GEOCENTRIC LATITUDE
      D1 = U(3) / DSQRT(U(1)**2 + U(2)**2 + U(3)**2)

C GEOGRAPHIC (GEODETIC) COORDINATES OF THE POINT
      RLAT = DATAN(AEBE2C * D1 / DSQRT(1. - D1 * D1))
      RLON = DATAN2(U(2),U(1))
      IERR = 0

      RETURN
      END

      SUBROUTINE GPOINT(RLAT,RLON,ALF,GAM,IERR)

      IMPLICIT NONE

C CALLING PARAMETERS
      REAL*8   RLAT
C GEOGRAPHIC LATITUDE IN RADIANS (INPUT)
      REAL*8   RLON
C GEOGRAPHIC LONGITUDE IN RADIANS (INPUT)
      REAL*8   ALF
C ELEVATION ANGLE IN RADIANS (OUTPUT)
      REAL*8   GAM
C SCAN ANGLE IN RADIANS (OUTPUT)
      INTEGER IERR
C OUTPUT STATUS; 0 - SUCCESSFUL COMPLETION,
C                1 - POINT WITH GIVEN LAT/LON IS INVISIBLE

      REAL*8       F(3)  , AEC   , FERC  , AEBE2C, AEBE3C, AEBE4C
      REAL*8       FT(3) , U(3)  , SING  , SLAT  , W1    , W2

C     INCLUDE 'ELCONS.INC'
C=========================== DIELCONS =============================

C $   (JR)
C $   THIS INCLUDE FILE IS PART OF THE SUPPLIED GVAR NAV SOFTWARE
C $
C $   DESCRIPTION:
C $
C $
C $$  DIELCONS = INCLUDE

      REAL*8 PI
           PARAMETER (PI=3.141592653589793D0)
      REAL*8 DEG
           PARAMETER (DEG=180.D0/PI)
      REAL*8 RAD
           PARAMETER (RAD=PI/180.D0)
C                    DEGREES TO RADIANS CONVERSION PI/180
      REAL*8 NOMORB
           PARAMETER (NOMORB=42164.365D0)
C                    NOMINAL RADIAL DISTANCE OF SATELLITE (KM)
      REAL*8 AE
           PARAMETER (AE=6378.137D0)
C                    EARTH EQUATORIAL RADIUS (KM)
      REAL*8 FER
           PARAMETER (FER=1.D0-(6356.7533D0/AE))
C                    EARTH FLATTENING COEFFICIENT = 1-(BE/AE)
      REAL*4 AEBE2
           PARAMETER (AEBE2=1.D0/(1.D0-FER)**2)
      REAL*4 AEBE3
           PARAMETER (AEBE3=AEBE2-1.)
      REAL*4 AEBE4
           PARAMETER (AEBE4=(1.D0-FER)**4-1.)
C========================= END DIELCONS ===========================
C     INCLUDE 'ELCOMM.INC'
C=========================== DIELCOMM =============================

C $   (JR)
C $   THIS INCLUDE FILE IS PART OF THE SUPPLIED GVAR NAV SOFTWARE
C $
C $   DESCRIPTION:
C $
C $
C $$  DIELCOMM = INCLUDE

C     COMMON VARIABLES
C
      REAL*8 XS(3)
C                      NORMALIZED S/C POSITION IN ECEF COORDINATES
      REAL*8 BT(3,3)
C                      ECEF TO INSTRUMENT COORDINATES TRANSFORMATION
      REAL*8  Q3
C                      USED IN SUBRTN LPOINT
      REAL*8 PITCH,ROLL,YAW
C                          PITCH,ROLL,YAW ANGLES OF INSTRUMENT (RAD)
      REAL*4 PMA,RMA
C                          PITCH,ROLL MISALIGNMENTS OF INSTRUMENT (RAD)
         COMMON /ELCOMM/ XS,BT,Q3,PITCH,ROLL,YAW,PMA,RMA
C========================= END DIELCOMM ===========================

      COMMON /RADCOM/ AEC   , FERC  , AEBE2C, AEBE3C, AEBE4C

C-----------------------------------------------------------------------

C COMPUTES SINUS OF GEOGRAPHIC (GEODETIC) LATITUDE
      SING = DSIN(RLAT)
      W1   = AEBE4C * SING * SING

C SINUS OF THE GEOCENTRIC LATITUDE
      SLAT = ((0.375 * W1 - 0.5) * W1 + 1.) * SING / AEBE2C

C COMPUTES LOCAL EARTH RADIUS AT SPECIFIED POINT
      W2 = SLAT * SLAT
      W1 = AEBE3C * W2
      W1 = (0.375 * W1 - 0.5) * W1 + 1.

C COMPUTES CARTESIAN COORDINATES OF THE POINT
      U(3) = SLAT * W1
      W2   = W1 * DSQRT(1. - W2)
      U(1) = W2 * DCOS(RLON)
      U(2) = W2 * DSIN(RLON)

C POINTING VECTOR FROM SATELLITE TO THE EARTH POINT
      F(1) = U(1) - XS(1)
      F(2) = U(2) - XS(2)
      F(3) = U(3) - XS(3)
      W2   = U(1) * SNGL(F(1)) + U(2) * SNGL(F(2))+ U(3) * SNGL(F(3))
     &                                                   * AEBE2C

C VERIFIES VISIBILITY OF THE POINT
      IF (W2.GT.0.) THEN

C INVISIBLE POINT ON THE EARTH
        IERR = 1
        ALF  = 99999.
        GAM  = 99999.
        RETURN
      ENDIF

C CONVERTS POINTING VECTOR TO INSTRUMENT COORDINATES
      FT(1) = BT(1,1) * F(1) + BT(2,1) * F(2) + BT(3,1) * F(3)
      FT(2) = BT(1,2) * F(1) + BT(2,2) * F(2) + BT(3,2) * F(3)
      FT(3) = BT(1,3) * F(1) + BT(2,3) * F(2) + BT(3,3) * F(3)

C CONVERTS POINTING VECTOR TO SCAN AND ELEVATION ANGLES AND
C CORRECTS FOR THE ROLL AND PITCH MISALIGNMENTS
      GAM  = DATAN(FT(1) / DSQRT(FT(2)**2 + FT(3)**2))
      ALF  = -DATAN(FT(2) / FT(3))
      W1   = DSIN(ALF)
      W2   = DCOS(GAM)
      ALF  = ALF + RMA * (1. - DCOS(ALF) / W2) + PMA * W1 *
     &                                   (1. / W2 + DTAN(GAM))
      GAM  = GAM - RMA * W1
      IERR = 0

      RETURN
      END
c
c--- an in-line version of gvrang (originally angles.for) must be
c--- included because of the common block, elcomm.
c
c--- this is required because the preprocessor, convdlm, renames
c--- all of the common blocks, appending the nav type.  the common
c--- block in angles.for was not being renamed, and therefor not
c--- being accessed.  one of the big problems with this solution
c--- is that applications cannot currently call angles directly
c--- for goes-8 navigation.  applications should use GVROPT('ANG '...)
c--- instead.
c

      SUBROUTINE GVRANG(JDAY,JTIME,XLAT,XLON,GHA,DEC,SATANG,SUNANG,RELAN
     *G)
C
C   THIS IS SSEC PROPRIETARY SOFTWARE - ITS USE IS RESTRICTED.
C $ SUBROUTINE ANGLES(JDAY,JTIME,XLAT,XLON,GHA,DEC,SATANG,SUNANG,RELANG)
C $ ANGLES - computes zenith angles of sun and satellite and relative
C $   azimuth angle              (DAS)
C $ INPUT:
C $   JDAY = (I) picture day (YYDDD)
C $   JTIME = (I) picture start time
C $   XLAT = (R) latitude of point
C $   XLON = (R) longitude of point
C $   GHA = (R) Greenwich hour angle of sun
C $   DEC = (R) declination of sun
C $ OUTPUT:
C $   SATANG = (R) zenith angle of satellite
C $   SUNANG = (R) zenith angle of sun
C $   RELANG = (R) relative angle
C $$ ANGLES = COMPUTATION, NAVIGATION
C ANGLES MOSHER 1074 WINLIB  ZENITH ANGLES TO SAT,SUN,AND REL AZIMUTH AN
C
C     INCLUDE 'ELCONS.INC'
C=========================== DIELCONS =============================

C $   (JR)
C $   THIS INCLUDE FILE IS PART OF THE SUPPLIED GVAR NAV SOFTWARE
C $
C $   DESCRIPTION:
C $
C $
C $$  DIELCONS = INCLUDE

      REAL*8 PI
           PARAMETER (PI=3.141592653589793D0)
      REAL*8 DEG
           PARAMETER (DEG=180.D0/PI)
      REAL*8 RAD
           PARAMETER (RAD=PI/180.D0)
C                    DEGREES TO RADIANS CONVERSION PI/180
      REAL*8 NOMORB
           PARAMETER (NOMORB=42164.365D0)
C                    NOMINAL RADIAL DISTANCE OF SATELLITE (KM)
      REAL*8 AE
           PARAMETER (AE=6378.137D0)
C                    EARTH EQUATORIAL RADIUS (KM)
      REAL*8 FER
           PARAMETER (FER=1.D0-(6356.7533D0/AE))
C                    EARTH FLATTENING COEFFICIENT = 1-(BE/AE)
      REAL*4 AEBE2
           PARAMETER (AEBE2=1.D0/(1.D0-FER)**2)
      REAL*4 AEBE3
           PARAMETER (AEBE3=AEBE2-1.)
      REAL*4 AEBE4
           PARAMETER (AEBE4=(1.D0-FER)**4-1.)
C========================= END DIELCONS ===========================
C     INCLUDE 'ELCOMM.INC'
C=========================== DIELCOMM =============================

C $   (JR)
C $   THIS INCLUDE FILE IS PART OF THE SUPPLIED GVAR NAV SOFTWARE
C $
C $   DESCRIPTION:
C $
C $
C $$  DIELCOMM = INCLUDE

C     COMMON VARIABLES
C
      REAL*8 XS(3)
C                      NORMALIZED S/C POSITION IN ECEF COORDINATES
      REAL*8 BT(3,3)
C                      ECEF TO INSTRUMENT COORDINATES TRANSFORMATION
      REAL*8  Q3
C                      USED IN SUBROUTINE LPOINT
      REAL*8 PITCH,ROLL,YAW
C                          PITCH,ROLL,YAW ANGLES OF INSTRUMENT (RAD)
      REAL*4 PMA,RMA
C                          PITCH,ROLL MISALIGNMENTS OF INSTRUMENT (RAD)
         COMMON /ELCOMM/ XS,BT,Q3,PITCH,ROLL,YAW,PMA,RMA
C========================= END DIELCOMM ===========================
      DATA IDAY/0/
      DATA R/6371.221/
      RDPDG=PI/180.0
      IF(IDAY.EQ.JDAY)GO TO 1
      IDAY=JDAY
      INORB=0
 1    PICTIM=FTIME(JTIME)
C
C   DETERMINE SATELLITE POSITION
C
      XSAT = XS(1) * AE
      YSAT = XS(2) * AE
      ZSAT = XS(3) * AE

      HEIGHT=SQRT(XSAT**2+YSAT**2+ZSAT**2)
      YLAT=RDPDG*XLAT
      YLAT=GEOLAT(YLAT,1)
      YLON=RDPDG*XLON
      SLAT=SIN(YLAT)
      CLAT=COS(YLAT)
      SLON=SIN(YLON)
      CLON=COS(YLON)
      XSAM=R*CLAT*CLON
      YSAM=R*CLAT*SLON
      ZSAM=R*SLAT
C
C   DETERMINE ZENITH ANGLE OF SUN
C
      SNLG=(-PICTIM)*PI/12.0-RDPDG*GHA
      SNDC=RDPDG*DEC
      COSDEC=COS(SNDC)
      US=COS(SNLG)*COSDEC
      VS=SIN(SNLG)*COSDEC
      WS=SIN(SNDC)
      SUNANG=ACOS((US*XSAM+VS*YSAM+WS*ZSAM)/R)/RDPDG
C
C   DETERMINE ZENITH ANGLE OF SATELLITE
C
      XVEC=XSAT-XSAM
      YVEC=YSAT-YSAM
      ZVEC=ZSAT-ZSAM
      XFACT=SQRT(XVEC**2+YVEC**2+ZVEC**2)
      SATANG=ACOS((XVEC*XSAM+YVEC*YSAM+ZVEC*ZSAM)/(R*XFACT))/RDPDG
C
C   DETERMINE RELATIVE ANGLE
C
      X1=CLAT*CLON
      Y1=CLAT*SLON
      Z1=SLAT
      X2=SLON
      Y2=-CLON
      X3=(-SLAT)*CLON
      Y3=(-SLAT)*SLON
      Z3=CLAT
      XC1=US-X1
      YC1=VS-Y1
      ZC1=WS-Z1
      XC2=XSAT/HEIGHT-X1
      YC2=YSAT/HEIGHT-Y1
      ZC2=ZSAT/HEIGHT-Z1
      XAN1=XC1*X3+YC1*Y3+ZC1*Z3
      XAN2=XC2*X3+YC2*Y3+ZC2*Z3
      YAN1=XC1*X2+YC1*Y2
      YAN2=XC2*X2+YC2*Y2
      XAN3=XAN1*XAN2+YAN1*YAN2
      YAN3=(-YAN1)*XAN2+XAN1*YAN2
      RELANG=ATAN2(YAN3,XAN3)/RDPDG
      RELANG=ABS(RELANG)
      RETURN
      END


! { dg-do compile }
! Bug 45794 - ICE: Segmentation fault in gfc_conv_procedure_call
subroutine foo (vector, mask)
  real :: vector(:)
  logical, optional :: mask(:)
  integer :: loc(1)
  if (present(mask)) then
    loc = maxloc(vector, mask)
  end if
end subroutine


text:
./top.sv
./memory.sv

ext: f


text:
        WRITE (6, 100)
100     FORMAT (1H ,11HHello World)
        CALL EXIT
        END

max_stars_repo_name: realBjornRoden/unix

max_issues_repo_name: realBjornRoden/unix

max_forks_repo_name: realBjornRoden/unix


text:
*DECK XGETUA
      SUBROUTINE XGETUA (IUNITA, N)
C***BEGIN PROLOGUE  XGETUA
C***PURPOSE  Return unit number(s) to which error messages are being
C            sent.
C***LIBRARY   SLATEC (XERROR)
C***CATEGORY  R3C
C***TYPE      ALL (XGETUA-A)
C***KEYWORDS  ERROR, XERROR
C***AUTHOR  Jones, R. E., (SNLA)
C***DESCRIPTION
C
C     Abstract
C        XGETUA may be called to determine the unit number or numbers
C        to which error messages are being sent.
C        These unit numbers may have been set by a call to XSETUN,
C        or a call to XSETUA, or may be a default value.
C
C     Description of Parameters
C      --Output--
C        IUNIT - an array of one to five unit numbers, depending
C                on the value of N.  A value of zero refers to the
C                default unit, as defined by the I1MACH machine
C                constant routine.  Only IUNIT(1),...,IUNIT(N) are
C                defined by XGETUA.  The values of IUNIT(N+1),...,
C                IUNIT(5) are not defined (for N .LT. 5) or altered
C                in any way by XGETUA.
C        N     - the number of units to which copies of the
C                error messages are being sent.  N will be in the
C                range from 1 to 5.
C
C***REFERENCES  R. E. Jones and D. K. Kahaner, XERROR, the SLATEC
C                 Error-handling Package, SAND82-0800, Sandia
C                 Laboratories, 1982.
C***ROUTINES CALLED  J4SAVE
C***REVISION HISTORY  (YYMMDD)
C   790801  DATE WRITTEN
C   861211  REVISION DATE from Version 3.2
C   891214  Prologue converted to Version 4.0 format.  (BAB)
C   920501  Reformatted the REFERENCES section.  (WRB)
C***END PROLOGUE  XGETUA
      DIMENSION IUNITA(5)
C***FIRST EXECUTABLE STATEMENT  XGETUA
      N = J4SAVE(5,0,.FALSE.)
      DO 30 I=1,N
         INDEX = I+4
         IF (I.EQ.1) INDEX = 3
         IUNITA(I) = J4SAVE(INDEX,0,.FALSE.)
   30 CONTINUE
      RETURN
      END


text:
C    @(#)swpepc.f	20.3 2/13/96
      subroutine swpepc (i, j)
C
      include 'ipfinc/epridc.inc'

      character temp * 120

      temp = epctl(i)
      epctl(i) = epctl(j)
      epctl(j) = temp
      return
      end

max_stars_repo_name: mbheinen/bpa-ipf-tsp

max_issues_repo_name: cuihantao/bpa-ipf-tsp

max_forks_repo_name: mbheinen/bpa-ipf-tsp


text:
Just a start here

Davis is hot and sunny from May to September, so sun protection is important.  Here are some ways to protect yourself from those UV rays.

Sunscreen

SPF (Sun Protection Factor) is A number on a scale for rating sunscreens. The SPF numbers on a package can range from as low as 2 to as high as 60. These numbers refer to the products ability to block out the suns rays.

The SPF rating is calculated by comparing the amount of time needed to produce a sunburn on protected skin to the amount of time needed to cause a sunburn on unprotected skin.

An SPF of 15 or higher is generally considered helpful in protecting the skin from UV rays.

recommendations:
  I like Coppertone Sport Ultra Sweatproof/Waterproof  it sprays on and drys instantly.  Complete coverage without getting your hands full of lotion. Also, the baby faces stick with 50+ SPF is great.  Users/JanelleAlvstadMattson Janelle

Hats

Broad brimmed hats are especially useful, keeping the sun off both your face and your neck.

Clothing

Most clothing has some sun protection factor, but you can also purchase clothing specifically made to have 30+ SPF.

Sunglasses

If you already have glasses, you can get prescription sunglasses or those interestinglooking temporary sunglasses (which look normal inside, but get darker in direct sunlight).

Shade

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
program pgmave

  use iso_varying_string
  use cmdline_arguments, only: get_options, bad_options, have_args, option_exists, &
      has_value, get_value, assignment(=), next_arg, num_args
  use string_functions, only: join, real
  use file_functions, only: exists, freeunit, stderr, stdout
  use pnm_class
  use precision

  implicit none

  integer :: i, error, nx, ny, missing_value

  character(len=2000) :: fname, outname

  logical :: initialised = .FALSE., fixed_name = .FALSE., have_missing = .FALSE.

  type (pnm_object) :: pnm
  type (varying_string) :: myoptions(3)

  integer(i8_kind), dimension(:,:), allocatable :: total

  integer, dimension(:,:), allocatable :: data, number

  ! These are our accepted command line options (see subroutine usage for
  ! an explanation)

  myoptions(1) = 'help'
  myoptions(2) = 'outfile'
  myoptions(3) = 'missing'

  ! This call parses the command line arguments for command line options
  call get_options(myoptions, error)

  ! Check we weren't passed duff options -- spit the dummy if we were
  if (error > 0) then
     write(stderr,*) 'ERROR! Unknown options: ',join(bad_options()," ")
     call usage
     STOP
  end if

  ! Check if we just want to print the usage
  if (option_exists('help')) then
     call usage
     STOP
  end if

  if (option_exists('outfile')) then
     ! We have specified an output filename as command line option
     ! Make sure we have a value
     if (.NOT. has_value('outfile')) then
        write(stderr,*) 'Option outfile must specify a filename!'
        call usage
        stop
     end if
     outname = ""
     outname = get_value('outfile')
     fixed_name = .TRUE.
  end if

  ! See if we have specified a maximum we wish to scale the data to
  if (option_exists('missing')) then
     ! Make sure we have a value
     if (.NOT. has_value('missing')) then
        write(stderr,*) 'Option missing must have a value!'
        call usage
        stop
     end if
     ! Get the maximum value we will scale to
     missing_value = get_value('missing')
     have_missing = .TRUE.
  end if

  if (num_args() < 1) then
     write(stderr,*) 'ERROR! Must supply pgm file(s) as command-line arguments'
     call usage
     STOP
  end if

  do while (have_args())

     ! The files are on the command line
     fname = next_arg()

     call read(pnm,trim(fname))

     if (.NOT. initialised) then

        nx = size(pnm,1)
        ny = size(pnm,2)

        ! Allocate the space required for the data
        allocate(data(nx,ny), total(nx,ny), number(nx,ny))

        initialised = .TRUE.

        number = 0

     end if

     data = pnm

     total = total + data

     if (have_missing) then
        where(data /= missing_value) number = number + 1
     else
        number = number + 1
     end if

  end do

  total = total / real(number)

  if (maxval(total) > (2**16 - 1)) then
     total = (real(total)/real(maxval(total))) * (2**16 - 1)
  end if

  pnm = int(total)

  if (.NOT. fixed_name) then
     ! Use the last input filename as a template for the output filename.
     ! Search for a ".pgm" suffix .. if we find one than
     ! replace it with "_ave.pgm", otherwise just slap "_ave.pgm" on the end
     i = index(fname,".pgm")
     if (i == 0) then
        ! Couldn't find a matching suffix, so look for the end of the
        ! filename string in the character variable
        i = verify(fname," ",back=.TRUE.) + 1
     end if
     outname = fname
     outname(i:) = "_ave.pgm"
  end if

  ! Write the pgm file.
  call write(pnm, outname)

  deallocate(data,total)

contains

  subroutine usage

    write(stderr,*)
    write(stderr,*) 'Usage: pgmave --help --outfile=<filename>'
    write(stderr,*)
    write(stderr,*) '  --help    - print this message.'
    write(stderr,*) '  --outfile - output file (default is replace ".pgm" with "_ave.pgm").'
    write(stderr,*) '  --missing - missing value (no data) which will be ignored when averaging.'
    write(stderr,*)

  end subroutine usage

end program pgmave

max_stars_repo_name: Welberry/pgmutils

max_issues_repo_name: Welberry/pgmutils

max_forks_repo_name: Welberry/pgmutils


text:
! { dg-do compile }
! PR20923 gfortran slow for large array constructors.
! Test case prepared from PR by Jerry DeLisle <jvdelisle@gcc.gnu.org>
program sel
    implicit none
    integer(kind=4),parameter      :: n=1000
    integer(kind=4)                :: i,j
    real(kind=4),dimension(n*n)    :: vect
    vect(:) = (/ ((( (i+j+3)),i=1,n),j=1,n) /)
end

max_stars_repo_name: maurizioabba/rose

max_issues_repo_name: sujankh/rose-matlab

max_forks_repo_name: sujankh/rose-matlab


text:
A Computer Engineering undergrad student.

welcome to the wiki are you a fan of Dtaq? Daubert


        !COMPILER-GENERATED INTERFACE MODULE: Fri May 22 16:58:29 2020
        ! This source file is for reference only and may not completely
        ! represent the generated interface used by the compiler.
        MODULE COMPUTE_BIELAK_CONDITIONS__genmod
          INTERFACE
            SUBROUTINE COMPUTE_BIELAK_CONDITIONS(COORD,IGLOB,NGLOB,IT,  &
     &DELTAT,DXUX,DXUZ,DZUX,DZUZ,VELOC_HORIZ,VELOC_VERT,X0_SOURCE,      &
     &Z0_SOURCE,A_PLANE,B_PLANE,C_PLANE,ANGLESOURCE,ANGLESOURCE_REFL,   &
     &C_INC,C_REFL,TIME_OFFSET,F0)
              INTEGER(KIND=4), INTENT(IN) :: NGLOB
              REAL(KIND=8), INTENT(IN) :: COORD(2,NGLOB)
              INTEGER(KIND=4), INTENT(IN) :: IGLOB
              INTEGER(KIND=4), INTENT(IN) :: IT
              REAL(KIND=8), INTENT(IN) :: DELTAT
              REAL(KIND=8), INTENT(OUT) :: DXUX
              REAL(KIND=8), INTENT(OUT) :: DXUZ
              REAL(KIND=8), INTENT(OUT) :: DZUX
              REAL(KIND=8), INTENT(OUT) :: DZUZ
              REAL(KIND=8), INTENT(OUT) :: VELOC_HORIZ
              REAL(KIND=8), INTENT(OUT) :: VELOC_VERT
              REAL(KIND=8) :: X0_SOURCE
              REAL(KIND=8) :: Z0_SOURCE
              REAL(KIND=8) :: A_PLANE(2)
              REAL(KIND=8) :: B_PLANE(2)
              REAL(KIND=8) :: C_PLANE(2)
              REAL(KIND=8) :: ANGLESOURCE
              REAL(KIND=8) :: ANGLESOURCE_REFL
              REAL(KIND=8) :: C_INC
              REAL(KIND=8) :: C_REFL
              REAL(KIND=8) :: TIME_OFFSET
              REAL(KIND=8) :: F0
            END SUBROUTINE COMPUTE_BIELAK_CONDITIONS
          END INTERFACE
        END MODULE COMPUTE_BIELAK_CONDITIONS__genmod


text:
C----------------------------------------------------------------------
        subroutine ftgtbh(iunit,ncols,nrows,nfield,ttype,tbcol,
     &                    tform,tunit,extnam,status)

C       OBSOLETE routine: should call ftghtb instead

        integer iunit,ncols,nrows,nfield,status,tbcol(*)
        character*(*) ttype(*),tform(*),tunit(*),extnam

        call ftghtb(iunit,0,ncols,nrows,nfield,ttype,
     &                    tbcol,tform,tunit,extnam,status)
        end

max_stars_repo_name: iraf-community/stsdas

max_issues_repo_name: kirxkirx/iraf

max_forks_repo_name: kirxkirx/iraf


text:
*DECK C1MERG
      SUBROUTINE C1MERG (TCOS, I1, M1, I2, M2, I3)
C***BEGIN PROLOGUE  C1MERG
C***SUBSIDIARY
C***PURPOSE  Merge two strings of complex numbers.  Each string is
C            ascending by the real part.
C***LIBRARY   SLATEC
C***TYPE      COMPLEX (S1MERG-S, D1MERG-D, C1MERG-C, I1MERG-I)
C***AUTHOR  (UNKNOWN)
C***DESCRIPTION
C
C   This subroutine merges two ascending strings of numbers in the
C   array TCOS.  The first string is of length M1 and starts at
C   TCOS(I1+1).  The second string is of length M2 and starts at
C   TCOS(I2+1).  The merged string goes into TCOS(I3+1).  The ordering
C   is on the real part.
C
C***SEE ALSO  CMGNBN
C***ROUTINES CALLED  CCOPY
C***REVISION HISTORY  (YYMMDD)
C   801001  DATE WRITTEN
C   891214  Prologue converted to Version 4.0 format.  (BAB)
C   900402  Added TYPE section.  (WRB)
C   910408  Modified to use IF-THEN-ELSE.  Make it look like MERGE
C           which was modified earlier due to compiler problems on
C           the IBM RS6000.  (RWC)
C   920130  Code name changed from CMPMRG to C1MERG.  (WRB)
C***END PROLOGUE  C1MERG
      INTEGER I1, I2, I3, M1, M2
      COMPLEX TCOS(*)
C
      INTEGER J1, J2, J3
C
C***FIRST EXECUTABLE STATEMENT  C1MERG
      IF (M1.EQ.0 .AND. M2.EQ.0) RETURN
C
      IF (M1.EQ.0 .AND. M2.NE.0) THEN
         CALL CCOPY (M2, TCOS(I2+1), 1, TCOS(I3+1), 1)
         RETURN
      ENDIF
C
      IF (M1.NE.0 .AND. M2.EQ.0) THEN
         CALL CCOPY (M1, TCOS(I1+1), 1, TCOS(I3+1), 1)
         RETURN
      ENDIF
C
      J1 = 1
      J2 = 1
      J3 = 1
C
   10 IF (REAL(TCOS(J1+I1)) .LE. REAL(TCOS(I2+J2))) THEN
         TCOS(I3+J3) = TCOS(I1+J1)
         J1 = J1+1
         IF (J1 .GT. M1) THEN
            CALL CCOPY (M2-J2+1, TCOS(I2+J2), 1, TCOS(I3+J3+1), 1)
            RETURN
         ENDIF
      ELSE
         TCOS(I3+J3) = TCOS(I2+J2)
         J2 = J2+1
         IF (J2 .GT. M2) THEN
            CALL CCOPY (M1-J1+1, TCOS(I1+J1), 1, TCOS(I3+J3+1), 1)
            RETURN
         ENDIF
      ENDIF
      J3 = J3+1
      GO TO 10
      END


module diagnostics_c

  implicit none

  private

contains

  !!!######################################################################
  subroutine enter_exit_c(sub_name, sub_name_len, state) bind(C, name="enter_exit_c")
    use iso_c_binding, only: c_ptr
    use utils_c, only: strncpy
    use diagnostics, only: enter_exit
    use other_consts, only: MAX_STRING_LEN
    implicit none
    integer,intent(in) :: state, sub_name_len
    type(c_ptr), value, intent(in) :: sub_name
    character(len=MAX_STRING_LEN) :: sub_name_f

    call strncpy(sub_name_f, sub_name, sub_name_len)
#if defined _WIN32 && defined __INTEL_COMPILER
    call so_enter_exit(sub_name_f, state)
#else
    call enter_exit(sub_name_f, state)
#endif

  end subroutine enter_exit_c

  !!!######################################################################
  subroutine set_diagnostics_on_c(state) bind(C, name="set_diagnostics_on_c")
    use diagnostics, only: set_diagnostics_on
    implicit none

    logical, intent(in) :: state

#if defined _WIN32 && defined __INTEL_COMPILER
    call so_set_diagnostics_on(state)
#else
    call set_diagnostics_on(state)
#endif

  end subroutine set_diagnostics_on_c

end module diagnostics_c


text:
      DOUBLE PRECISION FUNCTION DT_DAMG(It)

      IMPLICIT NONE
      DOUBLE PRECISION aam , dam , dgauni , DT_RNDM , gasuni , gauno ,
     &                 gaunon , uniga , v1 , vo , vv
      INTEGER i , io , It , nstab
      SAVE

C particle properties (BAMJET index convention),
C (dublicate of DTPART for HADRIN)
      INCLUDE 'inc/hnablt'

      DIMENSION gasuni(14)
      DATA gasuni/ - 1.D0 , -.98D0 , -.95D0 , -.87D0 , -.72D0 , -.48D0 ,
     &     -.17D0 , .17D0 , .48D0 , .72D0 , .87D0 , .95D0 , .98D0 ,
     &     1.D0/
      DATA gauno/2.352D0/
      DATA gaunon/2.4D0/
      DATA io/14/
      DATA nstab/23/

      i = 1
      IF ( It.LE.0 ) THEN
         DT_DAMG = 0.0D0
         GOTO 99999
      ELSE IF ( It.GT.nstab ) THEN
         dgauni = gauno*gaunon/DBLE(io-1)
         vv = DT_RNDM(dgauni)
         vv = vv*2.0D0 - 1.0D0 + 1.D-16
 50      vo = gasuni(i)
         i = i + 1
         v1 = gasuni(i)
         IF ( vv.GT.v1 ) GOTO 50
         uniga = dgauni*(DBLE(i)-2.0D0+(vv-vo+1.D-16)/(v1-vo)
     &           -(DBLE(io)-1.0D0)*0.5D0)
         dam = GAH(It)*uniga/gauno
         aam = AMH(It) + dam
         DT_DAMG = aam
         RETURN
      END IF
      DT_DAMG = AMH(It)
99999 END FUNCTION

max_stars_repo_name: pzhristov/DPMJET

max_issues_repo_name: pzhristov/DPMJET

max_forks_repo_name: pzhristov/DPMJET


text:
      SUBROUTINE I7COPY(P, Y, X)
C
C  ***  SET Y = X, WHERE X AND Y ARE INTEGER P-VECTORS  ***
C
      INTEGER P
      INTEGER X(P), Y(P)
C
      INTEGER I
C
      DO 10 I = 1, P
 10      Y(I) = X(I)
 999  RETURN
      END

max_stars_repo_name: sbwhitecap/clocc-hg

max_issues_repo_name: sbwhitecap/clocc-hg

max_forks_repo_name: sbwhitecap/clocc-hg


text:
C$Procedure ZZGFRRQ ( Private - GF, range rate between objects )

      SUBROUTINE ZZGFRRQ ( ET, TARG, OBS, ABCORR, VALUE )

C$ Abstract
C
C     SPICE Private routine intended solely for the support of SPICE
C     routines.  Users should not call this routine directly due
C     to the volatile nature of this routine.
C
C     Compute the apparent range rate between two ephemeris objects.
C
C$ Disclaimer
C
C     THIS SOFTWARE AND ANY RELATED MATERIALS WERE CREATED BY THE
C     CALIFORNIA INSTITUTE OF TECHNOLOGY (CALTECH) UNDER A U.S.
C     GOVERNMENT CONTRACT WITH THE NATIONAL AERONAUTICS AND SPACE
C     ADMINISTRATION (NASA). THE SOFTWARE IS TECHNOLOGY AND SOFTWARE
C     PUBLICLY AVAILABLE UNDER U.S. EXPORT LAWS AND IS PROVIDED "AS-IS"
C     TO THE RECIPIENT WITHOUT WARRANTY OF ANY KIND, INCLUDING ANY
C     WARRANTIES OF PERFORMANCE OR MERCHANTABILITY OR FITNESS FOR A
C     PARTICULAR USE OR PURPOSE (AS SET FORTH IN UNITED STATES UCC
C     SECTIONS 2312-2313) OR FOR ANY PURPOSE WHATSOEVER, FOR THE
C     SOFTWARE AND RELATED MATERIALS, HOWEVER USED.
C
C     IN NO EVENT SHALL CALTECH, ITS JET PROPULSION LABORATORY, OR NASA
C     BE LIABLE FOR ANY DAMAGES AND/OR COSTS, INCLUDING, BUT NOT
C     LIMITED TO, INCIDENTAL OR CONSEQUENTIAL DAMAGES OF ANY KIND,
C     INCLUDING ECONOMIC DAMAGE OR INJURY TO PROPERTY AND LOST PROFITS,
C     REGARDLESS OF WHETHER CALTECH, JPL, OR NASA BE ADVISED, HAVE
C     REASON TO KNOW, OR, IN FACT, SHALL KNOW OF THE POSSIBILITY.
C
C     RECIPIENT BEARS ALL RISK RELATING TO QUALITY AND PERFORMANCE OF
C     THE SOFTWARE AND ANY RELATED MATERIALS, AND AGREES TO INDEMNIFY
C     CALTECH AND NASA FOR ALL THIRD-PARTY CLAIMS RESULTING FROM THE
C     ACTIONS OF RECIPIENT IN THE USE OF THE SOFTWARE.
C
C$ Required_Reading
C
C     None.
C
C$ Keywords
C
C     RANGE RATE
C     EPHEMERIS
C     GEOMETRY
C     SEARCH
C
C$ Declarations

      IMPLICIT NONE

      DOUBLE PRECISION      ET
      INTEGER               TARG
      INTEGER               OBS
      CHARACTER*(*)         ABCORR
      DOUBLE PRECISION      VALUE

C$ Brief_I/O
C
C     VARIABLE  I/O  DESCRIPTION
C     --------  ---  --------------------------------------------------
C     ET         I   Ephemeris seconds past J2000 TDB
C     TARG       I   Target body ID
C     OBS        I   Observer body ID
C     ABCORR     I   Aberration correction flag
C     REF        I   Reference frame of the range rate
C     VALUE      O   Value of range rate between objects
C
C$ Detailed_Input
C
C     ET       is the time in ephemeris seconds past J2000 TDB at
C              which the range rate is to be measured.
C
C     TARG     the SPICE integer ID for the target body.
C
C     OBS      the SPICE integer ID for the observer.
C
C     ABCORR   the string description of the aberration corrections to
C              apply to the state evaluations to account for one-way
C              light time and stellar aberration.
C
C              Any aberration correction accepted by the SPICE
C              routine SPKEZR is accepted here. See the header
C              of SPKEZR for a detailed description of the
C              aberration correction options. For convenience,
C              the options are listed below:
C
C                 'NONE'     Apply no correction. Returns the "true"
C                            geometric state.
C
C                 'LT'       "Reception" case:  correct for
C                            one-way light time using a Newtonian
C                            formulation.
C
C                 'LT+S'     "Reception" case:  correct for
C                            one-way light time and stellar
C                            aberration using a Newtonian
C                            formulation.
C
C                 'CN'       "Reception" case:  converged
C                            Newtonian light time correction.
C
C                 'CN+S'     "Reception" case:  converged
C                            Newtonian light time and stellar
C                            aberration corrections.
C
C                 'XLT'      "Transmission" case:  correct for
C                            one-way light time using a Newtonian
C                            formulation.
C
C                 'XLT+S'    "Transmission" case:  correct for
C                            one-way light time and stellar
C                            aberration using a Newtonian
C                            formulation.
C
C                 'XCN'      "Transmission" case:  converged
C                            Newtonian light time correction.
C
C                 'XCN+S'    "Transmission" case:  converged
C                            Newtonian light time and stellar
C                            aberration corrections.
C
C              The ABCORR string lacks sensitivity to case, leading
C              and trailing blanks.
C
C$ Detailed_Output
C
C     VALUE    is the optionally light-time corrected range
C              rate of TARG observed from OBS.
C
C$ Parameters
C
C     None.
C
C$ Exceptions
C
C     None.
C
C$ Files
C
C     None.
C
C$ Particulars
C
C     This routine determines the apparent range rate of a target,
C     TARG, as seen from an observer, OBS, at epoch ET.
C
C$ Examples
C
C     None.
C
C$ Restrictions
C
C     None.
C
C$ Literature_References
C
C     None.
C
C$ Author_and_Institution
C
C     W.L. Taber     (JPL)
C     L.S. Elson     (JPL)
C
C$ Version
C
C-    SPICELIB version 1.0.1 08-JUL-2010 (EDW)
C
C        Minor typo correction to comments.
C
C-    SPICELIB version 1.0.0 09-JUN-2009 (NJB)(EDW)
C
C-&

C$ Index_Entries
C
C   compute the range rate between two objects.
C
C-&

C
C     SPICELIB functions.
C
      DOUBLE PRECISION      DVNORM

      LOGICAL               FAILED
      LOGICAL               RETURN

C
C     Local Variables.
C
      DOUBLE PRECISION      LT
      DOUBLE PRECISION      STATE    ( 6 )

      CHARACTER*(5)         REF

C
C     Standard SPICE error handling.
C
      IF ( RETURN () ) THEN
         RETURN
      ELSE
         CALL CHKIN  ( 'ZZGFRRQ' )
      END IF

C
C     We just want the range rate of TARG relative to OBS.
C     This calculation is invariant with respect to reference
C     frame; we use 'J2000'.
C
      REF = 'J2000'
      CALL SPKEZ ( TARG, ET, REF, ABCORR, OBS, STATE, LT )

      IF (  FAILED() ) THEN
          CALL CHKOUT (  'ZZGFRRQ' )
          RETURN
      END IF

C
C     Calculate the derivative from the STATE vector.
C
      VALUE = DVNORM( STATE )

C
C     All done.
C
      CALL CHKOUT (  'ZZGFRRQ' )
      RETURN

      END

max_stars_repo_name: agforero/FTFramework

max_issues_repo_name: agforero/fortran-testing-framework

max_forks_repo_name: agforero/fortran-testing-framework


FUNCTION gammln_s(xx)
  USE mo_kind, only: sp, dp
  USE mo_nrutil, ONLY : arth,assert
  IMPLICIT NONE
  REAL(SP), INTENT(IN) :: xx
  REAL(SP) :: gammln_s
  REAL(DP) :: tmp,x
  REAL(DP) :: stp = 2.5066282746310005_dp
  REAL(DP), DIMENSION(6) :: coef = (/76.18009172947146_dp,&
       -86.50532032941677_dp,24.01409824083091_dp,&
       -1.231739572450155_dp,0.1208650973866179e-2_dp,&
       -0.5395239384953e-5_dp/)
  call assert(xx > 0.0, 'gammln_s arg')
  x=xx
  tmp=x+5.5_dp
  tmp=(x+0.5_dp)*log(tmp)-tmp
  gammln_s=tmp+log(stp*(1.000000000190015_dp+&
       sum(coef(:)/arth(x+1.0_dp,1.0_dp,size(coef))))/x)
END FUNCTION gammln_s

FUNCTION gammln_v(xx)
  USE mo_kind, only: sp, dp, i4
  USE mo_nrutil, ONLY: assert
  IMPLICIT NONE
  INTEGER(I4) :: i
  REAL(SP), DIMENSION(:), INTENT(IN) :: xx
  REAL(SP), DIMENSION(size(xx)) :: gammln_v
  REAL(DP), DIMENSION(size(xx)) :: ser,tmp,x,y
  REAL(DP) :: stp = 2.5066282746310005_dp
  REAL(DP), DIMENSION(6) :: coef = (/76.18009172947146_dp,&
       -86.50532032941677_dp,24.01409824083091_dp,&
       -1.231739572450155_dp,0.1208650973866179e-2_dp,&
       -0.5395239384953e-5_dp/)
  if (size(xx) == 0) RETURN
  call assert(all(xx > 0.0), 'gammln_v arg')
  x=xx
  tmp=x+5.5_dp
  tmp=(x+0.5_dp)*log(tmp)-tmp
  ser=1.000000000190015_dp
  y=x
  do i=1,size(coef)
     y=y+1.0_dp
     ser=ser+coef(i)/y
  end do
  gammln_v=tmp+log(stp*ser/x)
END FUNCTION gammln_v

FUNCTION dgammln_s(z)
  !  Uses Lanczos-type approximation to ln(gamma) for z > 0.
  !  Reference:
  !       Lanczos, C. ''A precision approximation of the gamma
  !               function'', J. SIAM Numer. Anal., B, 1, 86-96, 1964.
  !  Accuracy: About 14 significant digits except for small regions
  !            in the vicinity of 1 and 2.
  !  Programmer: Alan Miller
  !              1 Creswick Street, Brighton, Vic. 3187, Australia
  !  e-mail: amiller @ bigpond.net.au
  !  Latest revision - 14 October 1996
  USE mo_kind, only: dp
  USE mo_nrutil, ONLY : assert
  IMPLICIT NONE
  REAL(DP), INTENT(IN) :: z
  REAL(DP) :: dgammln_s
  ! Local variables
  REAL(DP)  :: a(9) = (/ 0.9999999999995183_dp, 676.5203681218835_dp, &
       -1259.139216722289_dp, 771.3234287757674_dp, &
       -176.6150291498386_dp, 12.50734324009056_dp, &
       -0.1385710331296526_dp, 0.9934937113930748D-05, &
       0.1659470187408462D-06 /)
  REAL(DP)  :: zero = 0.0_dp,   &
       one = 1.0_dp, &
       lnsqrt2pi =  0.9189385332046727_dp, &
       half = 0.5_dp, &
       sixpt5 = 6.5_dp, &
       seven = 7.0_dp
  REAL(DP)  :: tmp
  INTEGER   :: j
  call assert(z > zero, 'dgammln_s arg')
  dgammln_s = zero
  tmp = z + seven
  DO j = 9, 2, -1
     dgammln_s = dgammln_s + a(j)/tmp
     tmp = tmp - one
  END DO
  dgammln_s = dgammln_s + a(1)
  dgammln_s = LOG(dgammln_s) + lnsqrt2pi - (z + sixpt5) + (z - half)*LOG(z + sixpt5)
  RETURN
END FUNCTION dgammln_s

FUNCTION dgammln_v(z)
  !  Uses Lanczos-type approximation to ln(gamma) for z > 0.
  !  Reference:
  !       Lanczos, C. ''A precision approximation of the gamma
  !               function'', J. SIAM Numer. Anal., B, 1, 86-96, 1964.
  !  Accuracy: About 14 significant digits except for small regions
  !            in the vicinity of 1 and 2.
  !  Programmer: Alan Miller
  !              1 Creswick Street, Brighton, Vic. 3187, Australia
  !  e-mail: amiller @ bigpond.net.au
  !  Latest revision - 14 October 1996
  USE mo_kind, only: dp
  USE mo_nrutil, ONLY : assert
  IMPLICIT NONE
  REAL(DP), DIMENSION(:), INTENT(IN) :: z
  REAL(DP), DIMENSION(size(z)) :: dgammln_v
  ! Local variables
  REAL(DP)  :: a(9) = (/ 0.9999999999995183_dp, 676.5203681218835_dp, &
       -1259.139216722289_dp, 771.3234287757674_dp, &
       -176.6150291498386_dp, 12.50734324009056_dp, &
       -0.1385710331296526_dp, 0.9934937113930748D-05, &
       0.1659470187408462D-06 /)
  REAL(DP)  :: zero = 0.0_dp,   &
       one = 1.0_dp, &
       lnsqrt2pi =  0.9189385332046727_dp, &
       half = 0.5_dp, &
       sixpt5 = 6.5_dp, &
       seven = 7.0_dp
  REAL(DP), DIMENSION(size(z)) :: tmp
  INTEGER   :: j
  if (size(z) == 0) RETURN
  call assert(all(z > 0.0_dp), 'dgammln_v arg')
  dgammln_v = zero
  tmp = z + seven
  DO j = 9, 2, -1
     dgammln_v = dgammln_v + a(j)/tmp
     tmp = tmp - one
  END DO
  dgammln_v = dgammln_v  + a(1)
  dgammln_v = LOG(dgammln_v) + lnsqrt2pi - (z + sixpt5) + (z - half)*LOG(z + sixpt5)
END FUNCTION dgammln_v


text:
FUNCTION almost_equal(x, gold, tol) RESULT(b)
  implicit none
  REAL,  intent(in) :: x
  REAL,  intent(in) :: gold
  REAL,  intent(in) :: tol
  LOGICAL              :: b
  b = ( gold * (1 - tol)  <= x ).AND.( x <= gold * (1+tol) )
END FUNCTION almost_equal

! This example attempts to map over a slice of an array, operate on it, then map it back.
program map_array_slices
   implicit none

   integer i, j, num_slices, num_values, a

   LOGICAL :: almost_equal
   real, pointer, contiguous, dimension(:,:) :: x, y
   real :: expected_value;
   a = 5
   num_slices = 10
   num_values = 1024

   allocate(x(num_values,num_slices))
   allocate(y(num_values,num_slices))

   ! initialize arrays
   do j = 1, num_slices
      x(:,j) = DBLE(j)
      y(:,j) = DBLE(j)
   end do

   !$omp parallel do
   do j = 1, num_slices
   ! Map and work on slice
      !$omp target teams distribute parallel do private(i) shared(a, x, y, j, num_values) map(tofrom:x(:,j)) map(to:y(:,j)) default(none)
      do i = 1, num_values
         x(i,j) = a*x(i,j) + y(i,j)
      end do
      !$omp end target teams distribute parallel do
   end do
   !$omp end parallel do

   do j = 1, num_slices
     expected_value = a * DBLE(j) + DBLE(j)
     do i = 1, num_values
        IF ( .NOT.almost_equal( x(i,j),expected_value, 0.1) ) THEN
            print*,  "x(", i, ",", j, ")",  " | cur", x(i,j),  "ref", expected_value
            STOP 120
        ENDIF
     end do
   end do

end program map_array_slices

max_stars_repo_name: TApplencourt/FGPU

max_issues_repo_name: TApplencourt/FGPU

max_forks_repo_name: TApplencourt/FGPU


!
! Distributed under the OSI-approved Apache License, Version 2.0.  See
!  accompanying file Copyright.txt for details.
!
!  adios2_mod.f90 : ADIOS2 Fortran bindings central module
!
!   Created on: Mar 13, 2017
!       Author: William F Godoy godoywf@ornl.gov
!

module adios2

    use adios2_parameters
    use adios2_adios
    use adios2_io
    use adios2_variable
    use adios2_engine

end module


text:
module mod_monolis_solver_PipeBiCGSTAB
  use mod_monolis_prm
  use mod_monolis_com
  use mod_monolis_mat
  use mod_monolis_precond
  use mod_monolis_matvec
  use mod_monolis_linalg
  use mod_monolis_linalg_util
  use mod_monolis_converge

  implicit none

contains

  subroutine monolis_solver_PipeBiCGSTAB(monoPRM, monoCOM, monoMAT)
    implicit none
    type(monolis_prm) :: monoPRM
    type(monolis_com) :: monoCOM
    type(monolis_mat) :: monoMAT
    integer(kind=kint) :: N, NP, NDOF, NNDOF
    integer(kind=kint) :: i, iter, iter_RR
    real(kind=kdouble) :: CG(5), RR, RW, RR1, RS, RZ, R2, QY, YY
    real(kind=kdouble) :: alpha, beta, omega, omega1, B2
    real(kind=kdouble), allocatable :: R(:), RT(:), R0(:), W0(:), WT(:), T(:), PT(:), S(:), ST(:)
    real(kind=kdouble), allocatable :: Z(:), ZT(:), Q(:), QT(:), Y(:), V(:)
    real(kind=kdouble), pointer :: B(:), X(:)
    logical :: is_converge

    N     = monoMAT%N
    NP    = monoMAT%NP
    NDOF  = monoMAT%NDOF
    NNDOF = N*NDOF
    X => monoMAT%X
    B => monoMAT%B
    iter_RR = 50

    if(monoPRM%is_init_x) X = 0.0d0

    allocate(R (NDOF*NP), source = 0.0d0)
    allocate(RT(NDOF*NP), source = 0.0d0)
    allocate(R0(NDOF*NP), source = 0.0d0)
    allocate(W0(NDOF*NP), source = 0.0d0)
    allocate(WT(NDOF*NP), source = 0.0d0)
    allocate(T (NDOF*NP), source = 0.0d0)
    allocate(PT(NDOF*NP), source = 0.0d0)
    allocate(S (NDOF*NP), source = 0.0d0)
    allocate(ST(NDOF*NP), source = 0.0d0)
    allocate(Z (NDOF*NP), source = 0.0d0)
    allocate(ZT(NDOF*NP), source = 0.0d0)
    allocate(Q (NDOF*NP), source = 0.0d0)
    allocate(QT(NDOF*NP), source = 0.0d0)
    allocate(Y (NDOF*NP), source = 0.0d0)
    allocate(V (NDOF*NP), source = 0.0d0)

    call monolis_set_converge(monoPRM, monoCOM, monoMAT, B, B2, is_converge, monoPRM%tdotp, monoPRM%tcomm_dotp)
    if(is_converge) return
    call monolis_residual(monoCOM, monoMAT, X, B, R, monoPRM%tspmv, monoPRM%tcomm_spmv)

    call monolis_vec_copy_R(N, NDOF, R, R0)

    call monolis_precond_apply(monoPRM, monoCOM, monoMAT, R, RT)
    call monolis_matvec(monoCOM, monoMAT, RT, W0, monoPRM%tspmv, monoPRM%tcomm_spmv)
    call monolis_precond_apply(monoPRM, monoCOM, monoMAT, W0, WT)
    call monolis_matvec(monoCOM, monoMAT, WT, T, monoPRM%tspmv, monoPRM%tcomm_spmv)
    call monolis_inner_product_R(monoCOM, N, NDOF, R, R , RR, monoPRM%tdotp, monoPRM%tcomm_dotp)
    call monolis_inner_product_R(monoCOM, N, NDOF, R, W0, RW, monoPRM%tdotp, monoPRM%tcomm_dotp)

    alpha = RR / RW
    beta  = 0.0d0
    omega = 0.0d0

    do iter = 1, monoPRM%maxiter
      do i = 1, NNDOF
        PT(i) = RT(i) + beta*(PT(i) - omega*ST(i))
        S (i) = W0(i) + beta*(S (i) - omega*Z (i))
        ST(i) = WT(i) + beta*(ST(i) - omega*ZT(i))
        Z (i) = T (i) + beta*(Z (i) - omega*V (i))
        Q (i) = R (i) - alpha*S (i)
        QT(i) = RT(i) - alpha*ST(i)
        Y (i) = W0(i) - alpha*Z (i)
      enddo

      call monolis_inner_product_R(monoCOM, N, NDOF, Q, Y, CG(1), monoPRM%tdotp, monoPRM%tcomm_dotp)
      call monolis_inner_product_R(monoCOM, N, NDOF, Y, Y, CG(2), monoPRM%tdotp, monoPRM%tcomm_dotp)

      call monolis_precond_apply(monoPRM, monoCOM, monoMAT, Z, ZT)
      call monolis_matvec(monoCOM, monoMAT, ZT, V, monoPRM%tspmv, monoPRM%tcomm_spmv)

      QY = CG(1)
      YY = CG(2)
      if(YY /= 0.0d0)then
        omega1 = QY / YY
      else
        omega1 = 0.0d0
      endif

      if(mod(iter, iter_RR) == 0)then
        do i = 1, NNDOF
          X(i) = X(i) + alpha*PT(i) + omega1*T(i)
        enddo
        call monolis_residual(monoCOM, monoMAT, X, B, R, monoPRM%tspmv, monoPRM%tcomm_dotp)
        call monolis_precond_apply(monoPRM, monoCOM, monoMAT, R, RT)
        call monolis_matvec(monoCOM, monoMAT, RT, W0, monoPRM%tspmv, monoPRM%tcomm_spmv)
      else
        do i = 1, NNDOF
          X (i) = X(i)  + alpha * PT(i) + omega1*QT(i)
          R (i) = Q (i) - omega1* Y (i)
          RT(i) = QT(i) - omega1*(WT(i) - alpha *ZT(i))
          W0(i) = Y (i) - omega1*(T (i) - alpha *V (i))
        enddo
      endif

      call monolis_inner_product_R(monoCOM, N, NDOF, R0, R , CG(1), monoPRM%tdotp, monoPRM%tcomm_dotp)
      call monolis_inner_product_R(monoCOM, N, NDOF, R0, W0, CG(2), monoPRM%tdotp, monoPRM%tcomm_dotp)
      call monolis_inner_product_R(monoCOM, N, NDOF, R0, S , CG(3), monoPRM%tdotp, monoPRM%tcomm_dotp)
      call monolis_inner_product_R(monoCOM, N, NDOF, R0, Z , CG(4), monoPRM%tdotp, monoPRM%tcomm_dotp)
      call monolis_inner_product_R(monoCOM, N, NDOF, R , R , CG(5), monoPRM%tdotp, monoPRM%tcomm_dotp)

      call monolis_precond_apply(monoPRM, monoCOM, monoMAT, W0, WT)
      call monolis_matvec(monoCOM, monoMAT, WT, T, monoPRM%tspmv, monoPRM%tcomm_spmv)

      RR1= CG(1)
      RW = CG(2)
      RS = CG(3)
      RZ = CG(4)
      R2 = CG(5)

      call monolis_check_converge_2(monoPRM, monoCOM, monoMAT, R2, B2, iter, is_converge)
      if(is_converge) exit

      beta  = (alpha*RR1) / (RR*omega1)
      alpha = RR1 / (RW + beta*RS - beta*omega1*RZ)
      omega = omega1
      RR    = RR1
    enddo

    call monolis_update_R(monoCOM, NDOF, X, monoPRM%tcomm_spmv)

    deallocate(R )
    deallocate(RT)
    deallocate(R0)
    deallocate(W0)
    deallocate(WT)
    deallocate(T )
    deallocate(PT)
    deallocate(S )
    deallocate(ST)
    deallocate(Z )
    deallocate(ZT)
    deallocate(Q )
    deallocate(QT)
    deallocate(Y )
    deallocate(V )
  end subroutine monolis_solver_PipeBiCGSTAB
end module mod_monolis_solver_PipeBiCGSTAB

max_stars_repo_name: nqomorita/monolis

max_issues_repo_name: nqomorita/monolis

max_forks_repo_name: nqomorita/monolis