C&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
C
C  VERSION 220124
C
C  THIS PROGRAM ASSEMBLES LOADING DISPLACEMENTS INTO A GEODYN "ATMOSPHERE"
C  LOADING FILE (WHICH MAY CONTAIN LOADING COMPONENTS COMING FROM SOURCES
C  OTHER THAN THE ATOSPHERE). LOADING DISPLACEMENTS MAY BE ASSEMBLED FROM
C  VARIOUS SOURCES THAT ARE  AVAILABLE AT THE WEBSITE: 
C  https://massloading.gsfc.nasa.gov  
C  THAT WEBSITE WAS CREATED BY LEONID PETROV AND CONTAINS THE FOLLOWING
C  LOADING MODELS:
C            ATMOSPHERE        GEOS-FPIT
C            ATMOSPHERE        MERRA2
C            HYDROLOGY         GEOS-FPIT
C            HYDROLOGY         MERRA2
C            NON-TIDAL OCEAN   MPIOM06
C            TIDAL OCEAN       GOT-48
C            TIDAL OCEAN       GOT-410.C
C            TIDAL OCEAN       FES-2012
C            TIDAL OCEAN       FES-2014B
C            TIDAL OCEAN       EQUIL01
C            TIDAL OCEAN       EQUIL02
C            
C
C   THIS PROGRAM WILL ASSEMBLE ATMOSPHERE, HYDROLOGY AND NON-TIDAL OCEAN 
C   COMPONENTS UPON USER REQUEST. TIDAL OCEAN LOADING IS NOT USED HERE
C   BECAUSE GEODYN ALREADY MODELS TIDAL OCEAN LOADING. 
C
C   NOTE THAT FOR EACH OF LEONID'S LOADING FILES, THERE ARE TWO OPTOPNS
C   FOR THE COORDINATE SYSTEM ORIGIN: (1) CENTER OF MASS OR (2)  "DEGREE
C   1 TERMS ONLY". THIS PROGTAM DEFAULTS TO CENTER OF MASS. THE USER CAN
C   SWITCH TO "DEGREE 1 TERMS ONLY" BY INSERTING A "DEG1" CARD IN THE
C   OPTIONS FILE.  
C
C
C   THE OUTPUT FILE FOR GEODYN IS AT 3 HOURS. THE GEOS-FPIT ATMOSPHERE MODEL
C   AS WELL AS THE MPIOM06 NON-TIDAL OCEAN MODEL AND BOTH HYDROLOGY MODELS
C   ARE ALSO AT 3 HOURS. THE CURRENT VERSION CAN NOT USE THE MERRA2 ATMOSPHERE
C   MODEL BECAUSE IT IS AT 6 HOURS.
C
C   REQUIRED INPUTS:
C
C     STAPOS    - THIS CAN BE A GEODYN SETUP DECK THAT CONTAINS STAPOS CARDS
C                 FOR THE STATIONS THAT WILL HAVE LOADING DISPLACEMENTS ON
C                 THE OUTPUT LOADING FILE. ALL THAT IS REALLY REQUIRED (AND
C                 USED) IS THE STAPOS SECTION.
C
C     OPTIONS   - INSTRUCTION CARD INPUT. HERE IS AN EXAMPLE OPTIONS 
C                 SETUP OF SIX CARDS:
C
C            START 20190421
C            END   20190423
C            ATM:GEOS-FPIT
C            LWS:GEOS-FPIT 
C            NTO:MPIOM06
C            PRINT
C            DEG1   <--- NOT THE DEFAULT
C
C
C                 THE OPTION CARDS CAN START IN ANY COLUMN.
C
C                 START AND END CARDS ARE REQUIRED. THESE CARDS MUST CONTAIN AN 
C                 8 DIGIT STRING OF THE FORM: YYYYMMDD. THE OUTPUT FILE WILL
C                 START THREE HOURS EARLIER THAN 0 HOURS OF THE START DAY AND END
C                 THREE HOURS LATER THAN 0 HOURS OF THE END DAY              
C
C                 IN THE ABOVE EXAMPLE ATMOSPHERE, HYDROLODY AND NON-TIDAL
C                 OCEAN WILL CONTRIBUTE  DISPLACEMENTS TO THE OUTPUT FILE.
C                 AT LEAST ONE OF THE THREE DISPLACEMENT TYPES (ATM, LWS OR NTO) 
C                 MUST BE REQUESTED. LWS:MERRA2 COULD BE REQUESTED INSTEAD OF 
C                 LWS:GEOSFPIT. IN FUTURE VERSIONS ATM:MERRA2 WILL BE AVAILABLE.  
C
C                 PRINT IS OPTIONAL. IF A PRINT CARD IS INCLUDED, THEN DISPLACEMENTS
C                 WILL BE PRINTED OUT FOR EVERY LOADING SITEAT EVERY THREE HOUR 
C                 INTERVAL
C
C     COMPILE:  ifort -i8 or gfortran -fdefault-integer-8
C
C&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
      IMPLICIT REAL*8(A-H,O-Z),LOGICAL(L)
      CHARACTER*120 COMND
      CHARACTER*120 COMND2
      CHARACTER*80 CARD,CARDD,CARDS
      CHARACTER*13 NAMATM,NAMLWS,NAMNTO,EPOCH
      CHARACTER*8 STANAM
      CHARACTER*2 IC
      REAL*8, ALLOCATABLE :: DISPL(:,:,:,:)
      PARAMETER(MAXSTA=500)
      PARAMETER(MAXD=5000)
      DIMENSION STACOR(3,MAXSTA)
      DIMENSION ITCOR(MAXSTA)
      DIMENSION ILODN(MAXSTA)
      DIMENSION STANAM(MAXSTA)
      DIMENSION IHOLD(MAXSTA)
      DIMENSION NMISS(MAXSTA)
      DIMENSION CARDD(MAXD)
      DIMENSION CARDS(MAXD)
      DIMENSION IYMD(1),IHM(1),SEC(1)
      DIMENSION FSEC(1)
      DIMENSION NMEP(3)
      DATA FSEC/0.D0/
      DATA NSTA/0/
      DATA IHOLD/MAXSTA*1/
      DATA NMISS/MAXSTA*0/
      DATA LATM/.FALSE./
      DATA LLWS/.FALSE./
      DATA LNTO/.FALSE./
      DATA LPRNT/.FALSE./
      DATA LMISST/.FALSE./
      DATA LD1/.FALSE./
      DATA NMEP/0,0,0/
      DATA STATLM/10000.D0/
      DATA STATLD/0.09D0/
C
      EPOCH='12345678_1234'
C
      INQUIRE(EXIST=LEXIST,FILE='STAPOS')
      IF(.NOT.LEXIST) THEN
        WRITE(6,6000)
        WRITE(6,6001)
        STOP
      ENDIF
      OPEN(UNIT=10,FILE='STAPOS',FORM='FORMATTED',STATUS='OLD')
      CALL GETSTA(10,MAXSTA,ILODN,STACOR,ITCOR,NSTA)
      CLOSE(10)
C
      INQUIRE(EXIST=LEXIST,FILE='OPTIONS')
      IF(.NOT.LEXIST) THEN
        WRITE(6,6000)
        WRITE(6,6002)
        STOP
      ENDIF
      OPEN(UNIT=10,FILE='OPTIONS',FORM='FORMATTED',STATUS='OLD')
      CALL GETOPT(10,MJDS1,MJDS2,LATM,LLWS,LNTO,NAMATM,NAMLWS,NAMNTO,
     &            LPRNT,LD1)
      CLOSE(10)
C
C
C
      COMND(1:44)='curl -k -O https://massloading.gsfc.nasa.gov'
      COMND(45:59)= '/imls/load_list'
      IF(LD1) COMND(45:62)='/imls/load_d1_list'
      COMND2(1:3)='rm '
C
C
C
C COMPUTE NUMBER OF EPOCHS AND ALLOCATE THE DISPLACEMENT ARRAY
      NEP=1+(MJDS2-MJDS1)/10800
      ALLOCATE(DISPL(3,NEP,4,NSTA))
      NZ=12*NSTA*NEP
      DO I=1,NZ
       DISPL(I,1,1,1)=-9999.D0
      ENDDO
      DO ISTA=1,NSTA
        DO IEP=1,NEP 
          DISPL(1,IEP,4,ISTA)=0.D0
        ENDDO
      ENDDO
C
C
C
C ATMOSPHERE SECTION THROUGH 100
C
C
C
C  IF ATMOSPHERE IS BEING USED GATHER UP FILES AND COLLECT DISPLACEMENTS   
      IF(.NOT.LATM) GO TO 100
      IF(LD1) THEN
        COMND(63:66)='/atm'
      ELSE
        COMND(60:63)='/atm'
      ENDIF
C ATM:GEOS-FPIT   atm/geosfpit/atm_geosfpit_20190421_0000.eph
C ATM:MERRA2      atm/merra2/atm_merra2_20190421_0000.eph
      IF(NAMATM.EQ.'ATM:GEOS-FPIT') THEN
        IF(.NOT.LD1) THEN
          COMND(64:103)='/geosfpit/atm_geosfpit_20190421_0000.eph'
          IEND=103
          IYS=87
          IHS=96
          IFS=74
        ELSE
          COMND(67:109)='/geosfpit/atm_geosfpit_d1_20190421_0000.eph'
          IEND=109
          IYS=93
          IHS=102
          IFS=77
        ENDIF
      ENDIF
      IF(NAMATM.EQ.'ATM:MERRA2   ') THEN
        IF(.NOT.LD1) THEN
          COMND(64:99)='/merra2/atm_merra2_20190421_0000.eph'
          IEND=99
          IYS=83
          IHS=92
          IFS=72
        ELSE
          COMND(67:105)='/merra2/atm_merra2_d1_20190421_0000.eph'
          IEND=105
          IYS=89
          IHS=98
          IFS=75
        ENDIF
      ENDIF
C
      MJDS=MJDS1
      DO 50 IEP=1,NEP
      CALL YMDHMS(MJDS,FSEC,IYMD,IHM,SEC,1)
      IYMD(1)=MOD(IYMD(1),1000000)
      IC='20'
      IY=IYMD(1)/10000
      IF(IY.GT.50) IC='19'
      COMND(IYS:IYS+1)=IC
      WRITE(COMND(IYS+2:IYS+7),7000) IYMD(1)
      JS=IYS+2     
      JF=IYS+7     
      DO J=JS,JF
       IF(COMND(J:J).EQ.' ') COMND(J:J)='0'
      ENDDO
      WRITE(COMND(IHS:IHS+3),7001) IHM(1)
      JS=IHS     
      JF=IHS+3     
      DO J=JS,JF
       IF(COMND(J:J).EQ.' ') COMND(J:J)='0'
      ENDDO
C
C
      CALL SYSTEM(COMND(1:IEND))
      OPEN(UNIT=10,FILE=COMND(IFS:IEND),FORM='FORMATTED',STATUS='OLD')
      CALL GETD(10,MAXD,NLODD,NLODS,CARDD,CARDS)
      CLOSE(10)
      ILEN=IEND-IFS+1
      COMND2(4:3+ILEN)=COMND(IFS:IEND)
      CALL SYSTEM(COMND2(1:3+ILEN))
C
C GET CORRESPONDENCE BETWEEN GEODYN STATIONS AND LOADING SITES
      DO J=1,NSTA
        DO JJ=1,8
          STANAM(JJ:JJ)=' '
        ENDDO
      ENDDO
C
      CALL GETNAM(NSTA,NLODS,CARDS,STACOR,ITCOR,STANAM,STATLM,STATLD)
C
C
C USING THE STATION NAMES GET THE STATION DISPLACEMENTS AT THIS EPOCH 
C
      CALL GETDSP(NSTA,NLODD,CARDD,STANAM,IEP,NEP,1,DISPL)
      MJDS=MJDS+10800
   50 CONTINUE
C
C CHECK FOR MISSING EPOCHS
      DO ISTA=1,NSTA
        DO IEP=1,NEP 
          IF(DISPL(1,IEP,1,ISTA).LT.-99.D0) NMEP(1)=NMEP(1)+1
        ENDDO
      ENDDO
C
C
C  END ATMOSPHERE SECTION
C
  100 CONTINUE
C
C
C
C HYDROLOGY SECTION THROUGH 200
C
C
      IF(.NOT.LLWS) GO TO 200
C
      IF(LD1) THEN
        COMND(63:66)='/lws'
      ELSE
        COMND(60:63)='/lws'
      ENDIF
C LWS:GEOS-FPIT   lws/geosfpit/lws_geosfpit_20190421_0030.eph
C LWS:MERRA2      lws/merra2/lws_merra2_20190421_0000.eph
      IF(NAMLWS.EQ.'LWS:GEOS-FPIT') THEN
        IF(.NOT.LD1) THEN
          COMND(64:103)='/geosfpit/lws_geosfpit_20190421_0030.eph'
          IEND=103
          IYS=87
          IHS=96
          IFS=74
        ELSE
          COMND(67:109)='/geosfpit/lws_geosfpit_d1_20190421_0030.eph'
          IEND=109
          IYS=93
          IHS=102
          IFS=77
        ENDIF
      ENDIF
      IF(NAMLWS.EQ.'LWS:MERRA2   ') THEN
        IF(.NOT.LD1) THEN
          COMND(64:99)='/merra2/lws_merra2_20190421_0000.eph'
          IEND=99
          IYS=83
          IHS=92
          IFS=72
        ELSE
          COMND(67:105)='/merra2/lws_merra2_d1_20190421_0000.eph'
          IEND=105
          IYS=89
          IHS=95
          IFS=75
        ENDIF
      ENDIF
C
      MJDS=MJDS1
      DO 150 IEP=1,NEP
      CALL YMDHMS(MJDS,FSEC,IYMD,IHM,SEC,1)
      IYMD(1)=MOD(IYMD(1),1000000)
      IC='20'
      IY=IYMD(1)/10000
      IF(IY.GT.50) IC='19'
      COMND(IYS:IYS+1)=IC
      WRITE(COMND(IYS+2:IYS+7),7000) IYMD(1)
      JS=IYS+2
      JF=IYS+7
      DO J=JS,JF
       IF(COMND(J:J).EQ.' ') COMND(J:J)='0'
      ENDDO
      WRITE(COMND(IHS:IHS+3),7001) IHM(1)
      JS=IHS
      JF=IHS+3
      DO J=JS,JF
       IF(COMND(J:J).EQ.' ') COMND(J:J)='0'
      ENDDO
      IF(NAMLWS.EQ.'LWS:GEOS-FPIT') COMND(IHS+2:IHS+2)='3'
C
C
      CALL SYSTEM(COMND(1:IEND))
      OPEN(UNIT=10,FILE=COMND(IFS:IEND),FORM='FORMATTED',STATUS='OLD')
      CALL GETD(10,MAXD,NLODD,NLODS,CARDD,CARDS)
      CLOSE(10)
      ILEN=IEND-IFS+1
      COMND2(4:3+ILEN)=COMND(IFS:IEND)
      CALL SYSTEM(COMND2(1:3+ILEN))
C
C GET CORRESPONDENCE BETWEEN GEODYN STATIONS AND LOADING SITES
      DO J=1,NSTA
        DO JJ=1,8
          STANAM(JJ:JJ)=' '
        ENDDO
      ENDDO
C
      CALL GETNAM(NSTA,NLODS,CARDS,STACOR,ITCOR,STANAM,STATLM,STATLD)
C
C
C USING THE STATION NAMES GET THE STATION DISPLACEMENTS AT THIS EPOCH
C
      CALL GETDSP(NSTA,NLODD,CARDD,STANAM,IEP,NEP,2,DISPL)
      MJDS=MJDS+10800
  150 CONTINUE
C
C
C
C CHECK FOR MISSING EPOCHS
      DO ISTA=1,NSTA
        DO IEP=1,NEP 
          IF(DISPL(1,IEP,2,ISTA).LT.-99.D0) NMEP(2)=NMEP(2)+1
        ENDDO
      ENDDO
C  END HYDROLGY SECTION
C
  200 CONTINUE
C
C
C
C NON-TIDAL OCEAN LOADING SECTION THROUGH 300
C
C
      IF(.NOT.LNTO) GO TO 300
C
      IF(LD1) THEN
        COMND(63:66)='/nto'
      ELSE
        COMND(60:63)='/nto'
      ENDIF
C NTO:MPIOM06     nto/mpiom06/nto_mpiom06_20190421_0000.eph
      IF(NAMNTO.EQ.'NTO:MPIOM06  ') THEN
        IF(.NOT.LD1) THEN
          COMND(64:101)='/mpiom06/nto_mpiom06_20190421_0000.eph'
          IEND=101
          IYS=85
          IHS=94
          IFS=73
        ELSE
          COMND(67:107)='/mpiom06/nto_mpiom06_d1_20190421_0000.eph'
          IEND=107
          IYS=91
          IHS=100
          IFS=76
        ENDIF
      ENDIF
C
      MJDS=MJDS1
      DO 250 IEP=1,NEP
      CALL YMDHMS(MJDS,FSEC,IYMD,IHM,SEC,1)
      IYMD(1)=MOD(IYMD(1),1000000)
      IC='20'
      IY=IYMD(1)/10000
      IF(IY.GT.50) IC='19'
      COMND(IYS:IYS+1)=IC
      WRITE(COMND(IYS+2:IYS+7),7000) IYMD(1)
      JS=IYS+2
      JF=IYS+7
      DO J=JS,JF
       IF(COMND(J:J).EQ.' ') COMND(J:J)='0'
      ENDDO
      WRITE(COMND(IHS:IHS+3),7001) IHM(1)
      JS=IHS
      JF=IHS+3
      DO J=JS,JF
       IF(COMND(J:J).EQ.' ') COMND(J:J)='0'
      ENDDO
C
C
      CALL SYSTEM(COMND(1:IEND))
      OPEN(UNIT=10,FILE=COMND(IFS:IEND),FORM='FORMATTED',STATUS='OLD')
      CALL GETD(10,MAXD,NLODD,NLODS,CARDD,CARDS)
      CLOSE(10)
      ILEN=IEND-IFS+1
      COMND2(4:3+ILEN)=COMND(IFS:IEND)
      CALL SYSTEM(COMND2(1:3+ILEN))
C
C GET CORRESPONDENCE BETWEEN GEODYN STATIONS AND LOADING SITES
      DO J=1,NSTA
        DO JJ=1,8
          STANAM(JJ:JJ)=' '
        ENDDO
      ENDDO
C
      CALL GETNAM(NSTA,NLODS,CARDS,STACOR,ITCOR,STANAM,STATLM,STATLD)
C
C
C USING THE STATION NAMES GET THE STATION DISPLACEMENTS AT THIS EPOCH
C
      CALL GETDSP(NSTA,NLODD,CARDD,STANAM,IEP,NEP,3,DISPL)
      MJDS=MJDS+10800
  250 CONTINUE
C
C
C CHECK FOR MISSING EPOCHS
      DO ISTA=1,NSTA
        DO IEP=1,NEP 
          IF(DISPL(1,IEP,3,ISTA).LT.-99.D0) NMEP(3)=NMEP(3)+1
        ENDDO
      ENDDO
C
C  END NON-TIDAL OCEAN LOADING SECTION
C
  300 CONTINUE
C
C
      IF(NMEP(1).GT.0) WRITE(6,6011) NMEP(1)
      IF(NMEP(2).GT.0) WRITE(6,6012) NMEP(2)
      IF(NMEP(3).GT.0) WRITE(6,6014) NMEP(3)
C
C SUM UP CORRECTIONS
C
      DO 400 ISTA=1,NSTA
      DO 390 IEP=1,NEP
      DISPL(1,IEP,4,ISTA)=0.D0
      DISPL(2,IEP,4,ISTA)=0.D0
      DISPL(3,IEP,4,ISTA)=0.D0
      LMISS=.FALSE.
      IF(DISPL(1,IEP,1,ISTA).GT.-99.D0) THEN
        DISPL(1,IEP,4,ISTA)=DISPL(1,IEP,4,ISTA)+DISPL(1,IEP,1,ISTA)
        DISPL(2,IEP,4,ISTA)=DISPL(2,IEP,4,ISTA)+DISPL(2,IEP,1,ISTA)
        DISPL(3,IEP,4,ISTA)=DISPL(3,IEP,4,ISTA)+DISPL(3,IEP,1,ISTA)
      ELSE
        LMISS=.TRUE.
      ENDIF
      IF(DISPL(1,IEP,2,ISTA).GT.-99.D0) THEN
        DISPL(1,IEP,4,ISTA)=DISPL(1,IEP,4,ISTA)+DISPL(1,IEP,2,ISTA)
        DISPL(2,IEP,4,ISTA)=DISPL(2,IEP,4,ISTA)+DISPL(2,IEP,2,ISTA)
        DISPL(3,IEP,4,ISTA)=DISPL(3,IEP,4,ISTA)+DISPL(3,IEP,2,ISTA)
      ELSE
        LMISS=.TRUE.
      ENDIF
      IF(DISPL(1,IEP,3,ISTA).GT.-99.D0) THEN
        DISPL(1,IEP,4,ISTA)=DISPL(1,IEP,4,ISTA)+DISPL(1,IEP,3,ISTA)
        DISPL(2,IEP,4,ISTA)=DISPL(2,IEP,4,ISTA)+DISPL(2,IEP,3,ISTA)
        DISPL(3,IEP,4,ISTA)=DISPL(3,IEP,4,ISTA)+DISPL(3,IEP,3,ISTA)
      ELSE
        LMISS=.TRUE.
      ENDIF
      IF(LMISS) THEN
        NMISS(ISTA)=NMISS(ISTA)+1
        LMISST=.TRUE.
      ENDIF
  390 CONTINUE
  400 CONTINUE
      IF(LMISST) THEN
        WRITE(6,6020)
        WRITE(6,6020)
        WRITE(6,6015)
        DO ISTA=1,NSTA
         IF(NMISS(ISTA).GT.0) WRITE(6,6016) ILODN(ISTA),NMISS(ISTA)
        ENDDO
        WRITE(6,6020)
        WRITE(6,6020)
      ENDIF
C
      IF(.NOT.LPRNT) GO TO 500
C PRINT OUT LOADING IF REQUESTED
      DO 490 ISTA=1,NSTA
      WRITE(6,6020)
      WRITE(6,6020)
      WRITE(6,6021) STANAM(ISTA),ILODN(ISTA)
C
      MJDS=MJDS1
      DO 450 IEP=1,NEP
      CALL YMDHMS(MJDS,FSEC,IYMD,IHM,SEC,1)
      IYMD(1)=MOD(IYMD(1),1000000)
      IC='20'
      IY=IYMD(1)/10000
      IF(IY.GT.50) IC='19'
      EPOCH(1:2)=IC
      WRITE(EPOCH(3:8),7000) IYMD
      WRITE(EPOCH(10:13),7001) IHM
      DO J=1,13
        IF(EPOCH(J:J).EQ.' ') EPOCH(J:J)='0'
      ENDDO
      WRITE(6,6022) EPOCH
      IF(LATM) WRITE(6,6023) (DISPL(JQ,IEP,1,ISTA),JQ=1,3)
      IF(LLWS) WRITE(6,6024) (DISPL(JQ,IEP,2,ISTA),JQ=1,3)
      IF(LNTO) WRITE(6,6025) (DISPL(JQ,IEP,3,ISTA),JQ=1,3)
      WRITE(6,6026) (DISPL(JQ,IEP,4,ISTA),JQ=1,3)
      MJDS=MJDS+10800
  450 CONTINUE
  490 CONTINUE
C
  500 CONTINUE
C
C  OUTPUT GEODYN LOADING FILE
C
      INQUIRE(EXIST=LEXIST,FILE='LOADFILE')
      IF(LEXIST) CALL SYSTEM('rm LOADFILE')
      OPEN(UNIT=30,FILE='LOADFILE',STATUS='NEW',FORM='UNFORMATTED')
C
      INDX=1
      INTRVL=10800
      WRITE(30) INDX,NSTA,(IHOLD(I),I=1,NSTA),MJDS1,MJDS2,INTRVL
C
      INDX=2
      WRITE(30) INDX,NSTA,(ILODN(I),I=1,NSTA)
C
      DO 600 ISTA=1,NSTA
      INDX=3
      IDUM=1
      WRITE(30) INDX,IDUM,ILODN(ISTA),NEP
      DO 590 IEP=1,NEP
      INDX=4
      I1=(DISPL(2,IEP,4,ISTA)*1.D9)
      I2=(DISPL(3,IEP,4,ISTA)*1.D9)
      I3=(DISPL(1,IEP,4,ISTA)*1.D9)
      WRITE(30) INDX,I1,I2,I3
 590  CONTINUE
 600  CONTINUE
      CLOSE(30)
      DEALLOCATE(DISPL)
      STOP
 6000 FORMAT(' EXECUTION TERMINATING')
 6001 FORMAT(' FILE STAPOS NOT FOUND')
 6002 FORMAT(' FILE OPTIONS NOT FOUND')
 6011 FORMAT(' ',I8,' MISSING ATMOSPHERE EPOCHS')
 6012 FORMAT(' ',I8,' MISSING HYDROLOGY EPOCHS')
 6014 FORMAT(' ',I8,' MISSING NON TIDAL OCEAN EPOCHS')
 6015 FORMAT(' THE FOLLOWING SITES HAVE MISSING DISPLACEMENT EPOCHS:') 
 6016 FORMAT(' SITE # ',I6,' HAS ',I6,' MISSING EPOCHS ')
 6020 FORMAT(' ')
 6021 FORMAT(' LOADING NAME: ',A8,' GEODYN SITE # ',I4)
 6022 FORMAT(' ',A13,'          Up            East          North')
 6023 FORMAT('  ATMOSPHERE  ',3F15.6)
 6024 FORMAT('  HYDROLOGY   ',3F15.6)
 6025 FORMAT('  NON-TIDAL O ',3F15.6)
 6026 FORMAT('  TOTAL       ',3F15.6)
 7000 FORMAT(I6)
 7001 FORMAT(I4)
      END
      SUBROUTINE GETDSP(NSTA,NLODD,CARDD,STANAM,IEP,NEP,ITYP,DISPL)
      IMPLICIT REAL*8(A-H,O-Z),LOGICAL(L)
      CHARACTER*80 CARD,CARDD(NLODD)
      CHARACTER*8 STANAM(NSTA)
      DIMENSION DISPL(3,NEP,4,NSTA)
C
      DO 100 ISTA=1,NSTA
      DO 90 JD=1,NLODD
      CARD=CARDD(JD)
      IF(STANAM(ISTA).NE.CARD(46:53)) GO TO 90
      READ(CARD(54:80),*) D1,D2,D3
      DISPL(1,IEP,ITYP,ISTA)=D1
      DISPL(2,IEP,ITYP,ISTA)=D2
      DISPL(3,IEP,ITYP,ISTA)=D3
      GO TO 100
  90  CONTINUE
 100  CONTINUE
      RETURN
      END
      SUBROUTINE GETNAM(NSTA,NLODS,CARDS,STACOR,ITCOR,STANAM,STATLM,
     &                  STATLD)
      IMPLICIT REAL*8(A-H,O-Z),LOGICAL(L)
      CHARACTER*80 CARD,CARDS(NLODS)
      CHARACTER*8 STANAM(NSTA)
      DIMENSION STACOR(3,NSTA),ITCOR(NSTA)
C
      DO 100 ISTA=1,NSTA
      TOL=STATLM
      IF(ITCOR(ISTA).EQ.1) THEN
       XLATR=STACOR(1,ISTA)*DACOS(-1.D0)/180.D0
       TOL=STATLD
       TOL2=TOL/DCOS(XLATR)
      ENDIF
      DO 50 J=1,NLODS
      CARD=CARDS(J)
      READ(CARD(13:54),*) C1,C2,C3
      IF(ITCOR(ISTA).EQ.1) THEN
        READ(CARD(56:80),*) C1,C2,C3
        TEST1=DABS(C1-STACOR(1,ISTA))
        TEST2=DABS(C2-STACOR(2,ISTA))
        IF(TEST1.LT.TOL.AND.TEST2.LT.TOL2) THEN
         TOL=TEST1
         TOL=TEST2
         STANAM(ISTA)=CARD(4:11)
        ENDIF
      ELSE
        READ(CARD(13:54),*) C1,C2,C3
        CC1=STACOR(1,ISTA)-C1
        CC2=STACOR(2,ISTA)-C2
        CC3=STACOR(3,ISTA)-C3
        TEST=DSQRT(CC1*CC1+CC2*CC2+CC3*CC3)
        IF(TEST.LT.TOL) THEN
         TOL=TEST
         STANAM(ISTA)=CARD(4:11)
        ENDIF
      ENDIF
  50  CONTINUE
 100  CONTINUE
      RETURN
      END
      SUBROUTINE GETD(IUN,MAXD,NLODD,NLODS,CARDD,CARDS)
      IMPLICIT REAL*8(A-H,O-Z),LOGICAL(L)
      CHARACTER*80 CARD,CARDD(MAXD),CARDS(MAXD)
C
      NLODD=0
      NLODS=0
  10  CONTINUE
      READ(IUN,5000,END=100) CARD
      IF(CARD(1:1).EQ.'S') THEN
        NLODS=NLODS+1
        IF(NLODS.GT.MAXD) THEN
          WRITE(6,6000)
          WRITE(6,6001) MAXD
        ENDIF
        CARDS(NLODS)=CARD
      ENDIF
      IF(CARD(1:1).EQ.'D') THEN
        NLODD=NLODD+1
        IF(NLODD.GT.MAXD) THEN
          WRITE(6,6000)
          WRITE(6,6001) MAXD
        ENDIF
        CARDD(NLODD)=CARD
      ENDIF
      GO TO 10
 100  CONTINUE
      RETURN
 5000 FORMAT(A80)
 6000 FORMAT(' EXECUTION TERMINATING IN GETD')
 6001 FORMAT(' MAXD= ',I8,' EXCEEDED')
      END
      SUBROUTINE GETOPT(IUN,MJDS1,MJDS2,LATM,LLWS,LNTO,NAMATM,NAMLWS,   &                 
     &                  NAMNTO,LPRNT,LD1)
      IMPLICIT REAL*8(A-H,O-Z),LOGICAL(L)
      CHARACTER*80 CARD
      CHARACTER*13 NAMATM,NAMLWS,NAMNTO
      DIMENSION IYMD(1),IHM(1),SEC(1)
      DIMENSION FSEC(1)
      DATA FSEC/0.D0/
C
      LSTART=.FALSE.
      LEND=.FALSE.
C
  10  CONTINUE
      READ(10,5000,END=100) CARD
C
      IF(.NOT.LSTART) THEN
        DO I=1,70
          IF(CARD(I:I+4).EQ.'START') THEN
            READ(CARD(I+6:80),*) ISTART
            IF(ISTART.LT.19700101.OR.ISTART.GT.20301231)  THEN
               WRITE(6,6000)
               WRITE(6,6006) ISTART
               STOP
            ENDIF
            LSTART=.TRUE.
            GO TO 10
          ENDIF
        ENDDO
      ENDIF
C
      IF(.NOT.LEND) THEN
        DO I=1,70
          IF(CARD(I:I+2).EQ.'END') THEN
            READ(CARD(I+6:80),*) IEND
            IF(IEND.LT.19700101.OR.IEND.GT.20301231)  THEN
               WRITE(6,6000)
               WRITE(6,6007) IEND
               STOP
            ENDIF
            LEND=.TRUE.
            GO TO 10
          ENDIF
        ENDDO
      ENDIF
C
C ATM:GEOS-FPIT   atm/geosfpit/atm_geosfpit_20190421_0000.eph
C ATM:MERRA2      atm/merra2/atm_merra2_20190421_0000.eph
C
C LWS:GEOS-FPIT   lws/geosfpit/lws_geosfpit_20190421_0030.eph
C LWS:MERRA2      lws/merra2/lws_merra2_20190421_0000.eph
C
C NTO:MPIOM06     nto/mpiom06/nto_mpiom06_20190421_0000.eph
C
      DO I=1,68
        IF(CARD(I:I+4).EQ.'PRINT') THEN
          LPRNT=.TRUE.
          GO TO 10
        ENDIF
        IF(CARD(I:I+3).EQ.'DEG1') THEN
          LD1=.TRUE.
          GO TO 10
        ENDIF
        IF(CARD(I:I+12).EQ.'ATM:GEOS-FPIT') THEN
          LATM=.TRUE.
          NAMATM='ATM:GEOS-FPIT'
          GO TO 10
        ENDIF
        IF(CARD(I:I+9).EQ.'ATM:MERRA2') THEN
          LATM=.TRUE.
          NAMATM='ATM:MERRA2   '
          WRITE(6,6000)
          WRITE(6,6020)
          WRITE(6,6021)
          WRITE(6,6022)
          STOP
          GO TO 10
        ENDIF
        IF(CARD(I:I+12).EQ.'LWS:GEOS-FPIT') THEN
          LLWS=.TRUE.
          NAMLWS='LWS:GEOS-FPIT'
          GO TO 10
        ENDIF
        IF(CARD(I:I+9).EQ.'LWS:MERRA2') THEN
          LLWS=.TRUE.
          NAMLWS='LWS:MERRA2   '
          GO TO 10
        ENDIF
        IF(CARD(I:I+10).EQ.'NTO:MPIOM06') THEN
          LNTO=.TRUE.
          NAMNTO='NTO:MPIOM06  '
          GO TO 10
        ENDIF
      ENDDO 
 100  CONTINUE
      IF(.NOT.LSTART) THEN
        WRITE(6,6000)
        WRITE(6,6001)
        STOP
      ENDIF
      IF(.NOT.LSTART) THEN
        WRITE(6,6000)
        WRITE(6,6002)
        STOP
      ENDIF
      IF(.NOT.LATM.AND..NOT.LLWS.AND..NOT.LNTO) THEN
        WRITE(6,6000)
        WRITE(6,6003)
        STOP
      ENDIF
C
C AT THIS POINT THERE IS A VALID SETUP DECK
      WRITE(6,6004)
      IF(LATM) WRITE(6,6005) NAMATM
      IF(LLWS) WRITE(6,6005) NAMLWS
      IF(LNTO) WRITE(6,6005) NAMNTO
      WRITE(6,6008)
      WRITE(6,6009)
C
C
      IYIMID=MOD(ISTART,1000000)
      IHIMIS=0
      CALL YMDTIS(IYIMID,IHIMIS,ISEC)
      MJDS1=ISEC-10800
      IYIMID=MOD(IEND,1000000)
      IHIMIS=0
      CALL YMDTIS(IYIMID,IHIMIS,ISEC)
      MJDS2=ISEC+10800
C
      CALL YMDHMS(MJDS1,FSEC,IYMD,IHM,SEC,1)
      IYMD(1)=MOD(IYMD(1),1000000)
      WRITE(6,6010) IYMD,IHM
      CALL YMDHMS(MJDS2,FSEC,IYMD,IHM,SEC,1)
      IYMD(1)=MOD(IYMD(1),1000000)
      WRITE(6,6011) IYMD,IHM
C
      NEP=1+(MJDS2-MJDS1)/10800
      WRITE(6,6012) NEP
C
      WRITE(6,6014)
      WRITE(6,6014)
C
      RETURN
 5000 FORMAT(A80)
 6000 FORMAT(' EXECUTION TERMINATING IN SUBROUTINE GETOPT')
 6001 FORMAT(' NO START CARD PROVIDED IN OPTIONS DATA SET') 
 6002 FORMAT(' NO START CARD PROVIDED IN OPTIONS DATA SET') 
 6003 FORMAT(' NO LOADING MODELS REQUESTED IN OPTIONS DATA SET') 
 6004 FORMAT(' THE FOLLOWING LOADING MODELS HAVE BEEN REQUESTED:')
 6005 FORMAT(' ',A13)
 6006 FORMAT(' INVALID START TIME: ',I14)
 6007 FORMAT(' INVALID END TIME:   ',I14)
 6008 FORMAT(' A FILE OF DISPLACEMENTS WILL BE CONSTRUCTED FROM ')
 6009 FORMAT(' THESE MODELS AT 3 HOUR INTERVALS')
 6010 FORMAT(' START TIME OF FILE ',I6,2X,I4)
 6011 FORMAT('   END TIME OF FILE ',I6,2X,I4)
 6012 FORMAT(' EACH SITE WILL HAVE ',I6,' DISPLACEMENT EPOCHS')
 6014 FORMAT(' ')
 6020 FORMAT(' MERRA2 ATMOSPHERE SELECTED')
 6021 FORMAT(' THIS PROGRAM WILL NEED A FEW UPGRADES TO ACCOMMODATE')
 6022 FORMAT(' THE MERRA2 ATMOSPHERE 6 HOUR INTERCAL')
      END
      SUBROUTINE GETSTA(IUN,MAXSTA,ILODN,STACOR,ITCOR,NSTA)
      IMPLICIT REAL*8(A-H,O-Z),LOGICAL(L)
      CHARACTER*80 CARD
      DIMENSION STACOR(3,MAXSTA)
      DIMENSION ITCOR(MAXSTA)
      DIMENSION ILODN(MAXSTA)
      DATA TOL/2000.D0/
C
      LSTAP=.FALSE.
      SEC=1.D0/3600.D0
      XMIN=1.D0/60.D0
C
  10  CONTINUE
C
      READ(IUN,5000,END=100) CARD
      IF(CARD(1:6).EQ.'STAPOS') THEN
        LSTAP=.TRUE.
        GO TO 10
      ENDIF
      IF(.NOT.LSTAP) GO TO 10
      IF(CARD(1:8).EQ.'ELCUTOFF') GO TO 10
      IF(CARD(1:8).EQ.'INSTRMNT') GO TO 10
      IF(CARD(1:8).EQ.'CONSTADJ') GO TO 10
      IF(CARD(1:8).EQ.'ADJUSTED') GO TO 10
      IF(CARD(1:8).EQ.'CONSTEND') GO TO 10
      IF(CARD(1:8).EQ.'EXTRAGEO') GO TO 10
      IF(CARD(1:8).EQ.'GEODETIC') GO TO 10
      IF(CARD(1:5).EQ.'FIXED') GO TO 10
      IF(CARD(1:6).EQ.'STAVEL') GO TO 10
      IF(CARD(1:6).EQ.'SIGVEL') GO TO 10
      IF(CARD(1:6).EQ.'TIMVEL') GO TO 10
      IF(CARD(1:6).EQ.'STATH2') GO TO 10
      IF(CARD(1:6).EQ.'STATL2') GO TO 10
      IF(CARD(1:6).EQ.'CORREL') GO TO 10
      IF(CARD(1:6).EQ.'ENDSTA') THEN
         LSTAP=.FALSE.
         GO TO 10
      ENDIF
C
C CHECK COLUMNS 21-70
      DO I=21,70
        CALL NUMCHK(CARD(I:I),LOK,LDIG) 
        IF(.NOT.LOK) GO TO 10
      ENDDO
C
C CHECK 67-70 FOR INTEGER
      LDIGT=.FALSE.
      DO I=67,70
        CALL NUMCHK(CARD(I:I),LOK,LDIG) 
        IF(LDIG) LDIGT=.TRUE.
        IF(.NOT.LDIG.AND.CARD(I:I).NE.' ') GO TO 10
      ENDDO
      IF(.NOT.LDIGT) GO TO 10
      READ(CARD(67:70),*) JLODN
C MAKE SURE THERE IS A VALID LOADING SITE NUMBER
      IF(JLODN.LE.0) GO TO 10
C CHECK WHETHER ANOTHER STATION ALREADY HAS THIS SITE NUMBER
C IF SO SKIP THIS STATION
      IF(NSTA.GT.0) THEN
        DO ISTA=1,NSTA
          IF(ILODN(ISTA).EQ.JLODN) GO TO 10
        ENDDO
      ENDIF
C CHECK STATION COORDINATES
      READ(CARD(21:65),5001) C1,C2,C3
      TEST=DSQRT(C1*C1+C2*C2+C3*C3)
      IF(TEST.GT.6382000.D0) GO TO 10
      NSTA=NSTA+1
      IF(NSTA.GT.MAXSTA) THEN
        WRITE(6,6000)
        WRITE(6,6001) MAXSTA
        STOP
      ENDIF
      ILODN(NSTA)=JLODN
      IF(TEST.LT.6355000.D0) THEN
C COORDINATES MUST BE GEODETIC
C CHANGE MMSS TO DECIMAL IN LAT AND LONG
        ITCOR(NSTA)=1
C
        C=C1
        IC=C
        CP=(C-DFLOAT(IC))*1000.D0
        IMMSS=CP
        IMM=IMMSS/100
        ISS=MOD(IMMSS,100)
        STACOR(1,NSTA)=DFLOAT(IC)+DFLOAT(IMM)*XMIN+DFLOAT(ISS)*SEC
C
        C=C2
        IC=C
        CP=(C-DFLOAT(IC))*1000.D0
        IMMSS=CP
        IMM=IMMSS/100
        ISS=MOD(IMMSS,100)
        STACOR(2,NSTA)=DFLOAT(IC)+DFLOAT(IMM)*XMIN+DFLOAT(ISS)*SEC
C
        C=C3
        IC=C
        CP=(C-DFLOAT(IC))*1000.D0
        IMMSS=CP
        IMM=IMMSS/100
        ISS=MOD(IMMSS,100)
        STACOR(3,NSTA)=DFLOAT(IC)+DFLOAT(IMM)*XMIN+DFLOAT(ISS)*SEC
      ELSE
C COORDINATES MUST BE CARTESIAN
        ITCOR(NSTA)=0
        STACOR(1,NSTA)=C1
        STACOR(2,NSTA)=C2
        STACOR(3,NSTA)=C3
      ENDIF
      GO TO 10
 100  CONTINUE
      IF(NSTA.LE.0) THEN
         WRITE(6,6000)
         WRITE(6,6002)
         STOP
      ENDIF
      RETURN
 5000 FORMAT(A80)
 5001 FORMAT(3D15.6)
 6000 FORMAT(' EXECITION TERMINATING IN SUBROUTINE GETSTA')
 6001 FORMAT(' MAXSTA= ',I8,' EXCEEDED')
 6002 FORMAT(' NO STATIONS FOUND')
      END
      SUBROUTINE NUMCHK(C,LOK,LDIG) 
      IMPLICIT REAL*8(A-H,O-Z),LOGICAL(L)
      CHARACTER*10 DIG
      CHARACTER*1 C
C
      DIG='0123456789'
      LOK=.FALSE.
      LDIG=.FALSE.
C
      DO I=1,10
       IF(DIG(I:I).EQ.C) LDIG=.TRUE.
      ENDDO
      IF(LDIG) THEN
       LOK=.TRUE.
       RETURN
      ENDIF
      IF(C.EQ.' ') LOK=.TRUE.
      IF(C.EQ.'-') LOK=.TRUE.
      IF(C.EQ.'.') LOK=.TRUE.
      RETURN
      END
C$ADDYMD
      SUBROUTINE ADDYMD(IYMD,IDAY)
C*********1*********2*********3*********4*********5*********6*********7**
C* NAME              ADDYMD
C*
C* FUNCTION          TO ADD OR SUBTRACT DAYS FROM A DATE IN THE FORM
C*                   YYMMDD AND TO PROVIDE THE USER WITH THE NEW DATE
C*
C* INPUT PARAMETERS:
C*
C*
C*                  IYMD    - SIX DIGIT DATE IN THE FORM YYMMDD
C*
C*                  IDAY    - NUMBER OF DAYS TO BE ADDED OR SUBTRACTED
C*                            FROM INPUT DATE
C*
C* OUTPUT PARAMETERS:
C*                  IYMD    - NEW DATE AFTER IDAY ADDED TO INPUT DATE
C*
C*
C* RESTRICTIONS      ....old and new dates must lie within the years:
C*                          March 1, 1900 to Feb. 28, 2100
C*
C*********1*********2*********3*********4*********5*********6*********7**
      IMPLICIT DOUBLE PRECISION (A-H,O-Z), LOGICAL(L)
      SAVE
C*
C***********************************************************************
C* START OF EXECUTABLE CODE ********************************************
C***********************************************************************
C*
      izero = 0
      call mjdymd( modjd, iymd, izero, 2 )
      newjd = modjd +IDAY
      call mjdymd( newjd, iymd, izero, 1 )
      if(iymd.lt.491231)iymd=iymd+1000000
C*
      RETURN
      END
C$MJDYMD
      SUBROUTINE MJDYMD(MODJD,IYMD,IHMS,IFLAG)
C*********1*********2*********3*********4*********5*********6*********7**
C* MJDYMD           CREATION DATE 07/18/83     PGMR - D. SEVITSKI
C*                       UPDATED  10/30/84     PGMR - B EDDY
C*                       UPDATED  12/13/90     PGMR - J. MCCARTHY
C*
C* FUNCTION          CONVERT TIMES FROM MJD AND MJDS TO YMD AND
C*                   VICE VERSA
C*
C* INPUT PARAMETERS:
C*                   MODJD - MODIFIED JULIAN DATE        (IFLAG =1)
C*                   MODJD - MODIFIED JULIAN DAY SECONDS (IFLAG =3 OR 4)
C*                   IYMD  - YEAR,MONTH,DAY(YYMMDD)      (IFLAG =2)
C*                   IFLAG - FLAG INDICATING INPUT AND OUTPUT
C*                           TIME CONVERSION AS FOLLOWS:
C*                              1 = MJD    TO YYMMDD
C*                              2 = YYMMDD TO MJD
C*                              3 = MJDS   TO YYMMDD
C*                              4 = MJDS   TO YYMMDD HHMMSS
C*
C* OUTPUT PARAMETERS:
C*                   MODJD - MODIFIED JULIAN DATE        (IFLAG =2)
C*                   IYMD  - YEAR,MONTH,DAY(YYMMDD)      (IFLAG =1,3,4)
C*                   IHMS  - HOUR,MINUTE, SECOND(HHMMSS) (IFLAG =4)
C*
C*
C* RESTRICTIONS      ....all dates must lie within the years:
C*                          March 1, 1900 to Feb. 28, 2100
C*
C*
C* UPDATE HISTORY:
C*
C* 10/30/84  B EDDY     - OPTION TO OUTPUT HOUR MINUTES SECONDS ADDED
C* 12/13/90  J McCarthy - capability to handle dates in 21th century added
C*
C***********1*********2*********3*********4*********5*********6*********7
C*
      IMPLICIT DOUBLE PRECISION (A-H,O-Z),LOGICAL(L)
      SAVE


C*      parameter ( half =  0.5d0 ) , ( one = 1.d0 ),
C*     1          ( seven = 7.d0 )  , ( ten = 10.d0 ),
C*     2          ( twelve = 12.d0 ), ( d14 = 14.d0 )
      parameter ( half =  0.5d0 )
      parameter ( one = 1.d0 )
      parameter ( seven = 7.d0 )
      parameter ( ten = 10.d0 )
      parameter ( twelve = 12.d0 )
      parameter ( d14 = 14.d0 )


      parameter ( xjd0 = 2400000.5d0 )
      parameter ( d1537 =  1537.d0 )
      parameter ( d122  =   122.1d0 )
      parameter ( d36525 =  365.25d0 )
      parameter ( d30600 =  30.6001d0 )
      parameter ( d4715 = 4715.d0 )
      parameter ( ib = -15 )
      parameter ( d17209 = 1720996.5d0 )
      parameter (  d3600 = 3600.d0 )
      parameter (  d60 =  60.d0 )


      DATA TGTYMD/410106.D0/
      DATA TMGDN2/30000.D0/

C*      DATA MJD00/15019/,MD1901/15385/,MD1940/29629/
C*      DATA IYMD01/010101/,IYMD40/400101/

C***********************************************************************
C* START OF EXECUTABLE CODE ********************************************
C***********************************************************************

      IF(IFLAG.LT.1 .OR. IFLAG.GT.4) RETURN

      GOTO (100,200,300,300),IFLAG

C*     MJD TO YMD

100   CONTINUE

      xjd = modjd + xjd0
C     write(6,*)'mjdymd: xjd ', xjd
      c = INT( xjd + half ) + 1537
C     write(6,*)'mjdymd: c ',c
      nd = INT( (c - d122) / d36525 )
C*     write(6,*)'mjdymd: nd ', nd
      e = INT( d36525 * nd )
C*     write(6,*)'mjdymd: e ', e
      nf = INT( ( c - e ) / d30600 )

      id = c - e - INT( d30600 * nf ) + frac( xjd + half )
C*     write(6,*)'mjdymd: id ',id
      im = nf - 1 - 12 * INT( nf / 14 )
C*     write(6,*)'mjdymd: im ', im
      iy = nd - 4715 - INT(  ( 7 + im ) / 10 )
C*     write(6,*)'mjdymd: iy ', iy


C*     ....need this mod to get 2-digit output years

      iy = MOD( iy, 100 )
C*     write(6,*)'mjdymd: iy ', iy

      iymd = iy * 10000 + im * 100 + id

C*     write(6,*) ' iy, im, id ', iy, im, id, '  mjd ', mjd

      RETURN

C*     YMD TO MJD

200   CONTINUE
      IY=IYMD/10000
      IM=(IYMD-IY*10000)/100
      ID=IYMD-IY*10000-IM*100


      if( iy .lt. 1900 ) then
         if( iy .gt. 50 ) then
            y = 1900 + iy - 1
C*     write (6,*)'in 200 block: y ', y
         else
            y = 2000 + iy - 1
C*     write (6,*)'in 200 block: y ', y
         endif
      endif

      if( im .gt. 2 ) then
         m = im
C*   write (6,*)'in 200 block: m ', m
         y = y + 1
      else
         m = im + 12
C*   write (6,*)'in 200 block: m ', m
      endif

      xjd = INT( d36525 * y ) + INT( d30600 * (m + 1) ) + ib
     &    + d17209 + id
C*     write(6,*) ' xjd ', xjd

      modjd = xjd - xjd0
C*     write(6,*) ' modjd ', modjd

C*     ....below line not needed for iflag=2
Cccc      fsec = ih * d3600 + imin * d60 + sec
C*
C*     write(6,*) ' iy, im, id ', iy, im, id, '  mjd ', mjd

      RETURN
C*
C*     MJDS TO YYMMDD (AND HHMMSS IF IFLAG=4)
C*
300   CONTINUE
C*     CONVERT FROM MJDS TO MJD
C*     ....mjds = 86400 * ( jd - geodyn reference time in jd )
C*     ....     = 86400 * ( mjd - geodyn reference time in mjd )
      IDY=(MODJD/86400)+INT(TMGDN2+0.1)

C*     CONVERT FROM MJD TO YMD (SAME AS ABOVE)
C*
      xjd = IDY + xjd0
      c = INT( xjd + half ) + 1537
      nd = INT( (c - d122) / d36525 )
      e = INT( d36525 * nd )
      nf = INT( ( c - e ) / d30600 )

      id = c - e - INT( d30600 * nf ) + frac( xjd + half )
      im = nf - 1 - 12 * INT( nf / 14 )
      iy = nd - 4715 - INT(  ( 7 + im ) / 10 )

C*     ....need this mod to get 2-digit output years

      iy = MOD( iy, 100 )

      iymd = iy * 10000 + im * 100 + id

C*      write(6,*) ' iy, im, id ', iy, im, id, '  mjd ', mjd


      IF(IFLAG.EQ.3) RETURN

      ISEC  =MOD(MODJD,86400)
      IH    =ISEC/3600
      IM    =(ISEC-IH*3600)/60
      IS    =ISEC-IH*3600-IM*60
      IHMS  =IH*10000+IM*100+IS
      RETURN
      END
C$YMDHMS
      SUBROUTINE YMDHMS(MJDS,FSEC,IYMD,IHM,SEC,NPTS)
C*********1*********2*********3*********4*********5*********6*********7**
C* YMDHMS           00/00/00            0000.0    PGMR - ?
C*
C*
C* FUNCTION:
C*
C* I/O PARAMETERS:
C*
C*   NAME    I/O  A/S   DESCRIPTION OF PARAMETERS
C*   ------  ---  ---   ------------------------------------------------
C*   MJDS
C*   FSEC
C*   IYMD
C*   IHM
C*   SEC
C*   NPTS
C*
C*
C* COMMENTS:
C*
C*
C*********1*********2*********3*********4*********5*********6*********7**
      IMPLICIT DOUBLE PRECISION (A-H,O-Z),LOGICAL(L)
      SAVE
      DIMENSION FSEC(NPTS),IYMD(NPTS),IHM(NPTS),SEC(NPTS)
C*     DATA IYMDRF/410106/
      DATA SECMIN/6.0D1/
      DATA TGTYMD/410106.D0/
      DATA TMGDN2/30000.D0/
      DATA SECDAY/86400.D0/
C*
C************************************************************************
C* START OF EXECUTABLE CODE *********************************************
C************************************************************************
C*
C* COMPUTE MODIFIED JULIAN DAY = JULIAN DAY - GEODYN REFERENCE TIME IN JD
      MJD=(DBLE(MJDS)+FSEC(1))/SECDAY
C* ADD ELAPSED DAYS SINCE GEODYN TIME TO REFERENCE CALENDAR DATE
      IYMDRF=INT(TGTYMD+0.001)
      IYMD0=IYMDRF
      CALL ADDYMD(IYMD0,MJD)
C* COMPUTE START TIME ELAPSED SECONDS WITHIN DAY
      ISEC0=MJDS-MJD*86400
      SEC0=DBLE(ISEC0)
C* ADD TO ELAPSED SECONDS SINCE START TIME
      DO 1000 N=1,NPTS
      SEC(N)=SEC0+FSEC(N)
 1000 CONTINUE
C* COMPUTE INTEGRAL ELAPSED DAYS FROM START DATE
      DO 2000 N=1,NPTS
      IYMD(N)=SEC(N)/SECDAY
 2000 CONTINUE
C* SUBTRACT INTEGRAL DAYS FROM ELAPSED SECONDS
      DO 3000 N=1,NPTS
      SEC(N)=SEC(N)-DBLE(IYMD(N)*86400)
 3000 CONTINUE
C* COMPUTE INTEGRAL ELAPSED MINUTES WITHIN DAY
      DO 4000 N=1,NPTS
      IHM(N)=SEC(N)/SECMIN
 4000 CONTINUE
C* SUBTRACT INTEGRAL MINUTES FROM ELAPSED SECONDS
      DO 5000 N=1,NPTS
      SEC(N)=SEC(N)-DBLE(IHM(N)*60)
 5000 CONTINUE
C* CONVERT INTEGRAL MINUTES TO HOURS AND MINUTES
      DO 6000 N=1,NPTS
C*     IH=IHM(N)/60
C*     IHM(N)=IHM(N)-IH*60
C*     IHM(N)=IHM(N)+IH*100
      IHM(N)=(IHM(N)/60)*40+IHM(N)
 6000 CONTINUE
C* REPLACE INTEGRAL DAYS SINCE START DATE WITH CALENDAR DATE
      IDAYS=0
      IYMDAY=IYMD0
      DO 8000 N=1,NPTS
      IF(IYMD(N).EQ.IDAYS) GO TO 7000
      IDAYS=IYMD(N)
      IYMDAY=IYMD0
      CALL ADDYMD(IYMDAY,IDAYS)
 7000 CONTINUE
      IYMD(N)=IYMDAY
 8000 CONTINUE
      RETURN
      END
C$FRAC
      FUNCTION FRAC(x)
      IMPLICIT REAL*8 (A-H,O-Z),LOGICAL(L)
C
C
      ix = int(x)
      frac = x - ix
      return
      end
C$YMDTIS
      SUBROUTINE YMDTIS(IYIMID,IHIMIS,ISEC)
C********1*********2*********3*********4*********5*********6*********7**
C YMDTIS           82/11/30            0000.0    PGMR - D. ROWLANDS
C                  90/07/90                      PGMR - LUCIA TSAOUSSI
C FUNCTION:
C
C  OLD VERSION        11/30/82          D. ROWLANDS
C                    CALCULATE THE INTEGER NUMBER OF UT SECONDS
C                    SINCE JD 2430000.5D0 TO  TIME IYIMID,IHIMIS
C
C  NEW VERSION        07/31/90          LUCIA TSAOUSSI
C                    EXTEND CALCULATION TO 21ST CENTURY AND EXTEND
C                    INPUT DATE TO YYYY MM DD (INSTEAD OF YY MM DD)
C                    USE DIFFERENT ALGORITHM FOR LEAP YEARS!
C
C I/O PARAMETERS:
C
C   NAME    I/O  A/S   DESCRIPTION OF PARAMETERS
C   ------  ---  ---   ------------------------------------------------
C   IYIMID   I         REFERENCE DATE
C   IHIMIS   I         REFERENCE TIME
C   ISEC     O         INTEGER UT SECONDS FROM JD 2430000.5 TO
C                      TIME IYIMID IHIMIS
C
C   RESTRICTIONS:    DATE MUST BE IN 20TH OR 21ST CENTURY
C
C********1*********2*********3*********4*********5*********6*********7**
      IMPLICIT REAL (A-H,O-Z), LOGICAL (L)
      SAVE
      DIMENSION MONTAB(24),MYRAB(2)
      DATA MONTAB/0,31,59,90,120,151,181,212,243,273,304,334,
     .            0,31,60,91,121,152,182,213,244,274,305,335/
      DATA MYRAB/2415019,2451544/
C
C**********************************************************************
C START OF EXECUTABLE CODE ********************************************
C**********************************************************************
C
      IYR=IYIMID/10000
      IF (IYR.GT.99) THEN
        IYRH=MOD(IYR,100)
        IYRT=IYR/100
      ELSE
        IYRH=IYR
        IF (IYRH.GE.50) THEN
          IYRT=19
        ELSE
          IYRT=20
        END IF
      END IF
      IADD=0
      IYR4=(IYRH-1)/4
      IF(IYR.EQ.0) IYR4=-1
      IJD=365*IYRH+IYR4
      IMO=MOD(IYIMID,10000)
      IMO=IMO/100
      IF (IYRH.EQ.0) THEN
         IYEAR=IYRT*100+IYRH
         IF (MOD(IYEAR,400).EQ.0) IADD=12
      ELSE
         IF(MOD(IYRH,4).EQ.0) IADD=12
      END IF
      IJD=IJD+MONTAB(IMO+IADD)
      ID=MOD(IYIMID,100)
      IJD=IJD+ID+MYRAB(IYRT-18)
      IH=IHIMIS/10000
      IM=IHIMIS-IH*10000
      IM=IM/100
      IS=MOD(IHIMIS,100)
      ISEC=(IJD-2430000)*86400+IH*3600+IM*60+IS
      RETURN
      END
