PROGRAM TIDEPREDICTION

  ! -----------------------------------------------------
  ! PURPOSE: TIDE PREDICTION
  ! -----------------------------------------------------
  ! CALL WITH: gfortran tidePrediction.f90 -o tidePrediction
  ! -----------------------------------------------------
  IMPLICIT NONE

  INTEGER           :: cnt, i,n, m, noFiles, sizeDoodson, maxDegree
  CHARACTER*130     :: fileName, path, readit, file1, file2, modelName
  CHARACTER*130     :: fileDoodson, fileAdmittance, fileList, outputFile
  DOUBLE PRECISION  :: mjd, doodsonArguments(6), r, GM
  DOUBLE PRECISION, ALLOCATABLE :: doodson(:,:), admittance(:,:), coeffCosine(:,:), coeffSine(:,:)
  DOUBLE PRECISION, ALLOCATABLE :: factorSin(:), factorCos(:), thetaf(:), anm(:,:)

  path = '../../models/EOT20_OCN/'
  fileName = 'EOT20_OCN'
  modelName = 'exampleTidesAtMJD51544.5'
  outputFile = 'tides_at_mjd_51544.5.gfc'

  fileDoodson = TRIM(path)//TRIM(fileName)//'_002doodson.txt'
  fileAdmittance = TRIM(path)//TRIM(fileName)//'_003admittance.txt'
  fileList =  TRIM(path)//TRIM(fileName)//'_001fileList.txt'

  OPEN(21,file=fileDoodson,STATUS='OLD',ACTION='READ')
  OPEN(22,file=fileAdmittance,STATUS='OLD',ACTION='READ')
  OPEN(23,file=fileList,STATUS='OLD',ACTION='READ')

  sizeDoodson = 0
  DO
     READ(21,'(A70)',end=1000) readit
     sizeDoodson = sizeDoodson + 1
  END DO
1000 CONTINUE
  REWIND(21)

  noFiles=0
  DO
     READ(23,'(A70)',end=2000) readit
     noFiles = noFiles + 1
  END DO
2000 CONTINUE
  REWIND(23)

  ALLOCATE(doodson(sizeDoodson, 6), thetaf(sizeDoodson), admittance(noFiles, sizeDoodson))
  ALLOCATE(factorSin(noFiles), factorCos(noFiles))

  DO cnt=1,sizeDoodson
     READ(21,*) doodson(cnt,:)
  END DO
  CLOSE(21)

  DO cnt=1,noFiles
     READ(22,*) admittance(cnt,:)
  END DO
  CLOSE(22)

  mjd = 51544.5D0
  CALL DOODSON_ARGUMENTS(mjd, doodsonArguments)

  thetaf = MATMUL(doodson, doodsonArguments)
  factorCos = MATMUL(admittance, COS(thetaf(:)))
  factorSin = MATMUL(admittance, SIN(thetaf(:)))

  DEALLOCATE(doodson, admittance)

  DO cnt=1,noFiles
     READ(23, *) file1, file2
     file1 =  TRIM(path)//TRIM(file1)
     file2 =  TRIM(path)//TRIM(file2)

     IF (cnt .EQ. 1) THEN
        OPEN(24,file=file1,STATUS='OLD',ACTION='READ')
        DO
           READ(24,'(A130)',end=3000) readit
           IF (readit(1:10) .EQ. 'max_degree') THEN
              READ(readit(24:26),'(I3)') maxDegree
              CLOSE(24)
              ALLOCATE(coeffCosine(maxDegree+1, maxDegree+1))
              ALLOCATE(coeffSine(maxDegree+1, maxDegree+1))
              ALLOCATE(anm(maxDegree+1, maxDegree+1))
           END IF
        END DO
3000    CONTINUE
     END IF

     CALL READ_ICGEM(file1, maxDegree, coeffCosine, r, GM)
     CALL READ_ICGEM(file2, maxDegree, coeffSine, r, GM)

     anm = anm + factorCos(cnt) * coeffCosine + factorSin(cnt) * coeffSine
  END DO

  DEALLOCATE(coeffCosine, coeffSine, thetaf, factorSin, factorCos)

  CALL WRITE_ICGEM(outputFile, modelName, maxDegree, anm, r, GM)

CONTAINS
  SUBROUTINE READ_ICGEM(inputFile, maxDegree, gravityData, radius, gravityConstant)
    IMPLICIT NONE
    INTEGER                        :: n, m, stat
    INTEGER, INTENT(IN)            :: maxDegree
    CHARACTER*130, INTENT(IN)      :: inputFile
    CHARACTER*130                  :: header
    DOUBLE PRECISION               :: cnm, snm
    DOUBLE PRECISION, INTENT(OUT)  :: radius, gravityConstant
    DOUBLE PRECISION, INTENT(OUT)  :: gravityData(maxDegree+1,maxDegree+1)

    OPEN(900,file=inputFile,STATUS='OLD',ACTION='READ',IOSTAT=stat)

    IF (stat .EQ. 0) THEN
       DO
          READ(900,'(A130)',end=4000) header

          IF (header(1:3) .EQ. 'gfc' ) THEN
             READ(header(6:9),'(I4)') n
             READ(header(11:14),'(I4)') m
             READ(header(16:34),'(E19.12)') cnm
             READ(header(36:54),'(E19.12)') snm

             gravityData(n+1, m+1) = cnm
             IF (m .NE. 0) THEN
                gravityData(m, n+1) = snm
             END IF

          ELSEIF (header(1:6) .EQ. 'radius') THEN
             READ(header(24:39),'(E16.20)') radius

          ELSEIF (header(1:22) .EQ. 'earth_gravity_constant') THEN
             READ(header(24:39),'(E16.10)') gravityConstant
          END IF
       END DO
    ELSE
       WRITE(*,*) 'Can not find input file: ',inputfile
    END IF
4000 CONTINUE
    CLOSE(900)

  END SUBROUTINE READ_ICGEM

END PROGRAM TIDEPREDICTION

SUBROUTINE DOODSON_ARGUMENTS(mjd, doodson)
  IMPLICIT NONE
  DOUBLE PRECISION              :: t, gmst
  DOUBLE PRECISION              :: era
  DOUBLE PRECISION, INTENT(IN)  :: mjd
  DOUBLE PRECISION, INTENT(OUT) :: doodson(6)
  DOUBLE PRECISION, PARAMETER   :: PI = 3.141592653589793238462643383279502884D0

  t = (mjd - 51544.5D0) / 365250.D0
  era = 2D0 * PI * (MODULO(mjd-51544.5D0, 1D0) + 0.7790572732640D0 + 0.00273781191135448D0 * (mjd-51544.5D0))
  gmst=era+(0.014506D0+(4612.156534D0+(1.391581D0+(-0.00000044D0+(-0.000029956D0+(-0.0000000368D0*t))*t)*t)*t)*t)*PI/(180D0*3600D0)

  doodson(2) = (218.31664562999D0 + (4812678.81195750D0 + (-0.14663889D0 + ( 0.00185140D0 + (-0.00015355D0*t))*t)*t)*t) * PI / 180D0
  doodson(3) = (280.46645016002D0 + ( 360007.69748806D0 + ( 0.03032222D0 + ( 0.00002000D0 + (-0.00006532D0*t))*t)*t)*t) * PI / 180D0
  doodson(4) = ( 83.35324311998D0 + (  40690.13635250D0 + (-1.03217222D0 + (-0.01249168D0 + ( 0.00052655D0*t))*t)*t)*t) * PI / 180D0
  doodson(5) = (234.95544499000D0 + (  19341.36261972D0 + (-0.20756111D0 + (-0.00213942D0 + ( 0.00016501D0*t))*t)*t)*t) * PI / 180D0
  doodson(6) = (282.93734098001D0 + (     17.19457667D0 + ( 0.04568889D0 + (-0.00001776D0 + (-0.00003323D0*t))*t)*t)*t) * PI / 180D0
  doodson(1) = gmst + PI - doodson(2)

END SUBROUTINE DOODSON_ARGUMENTS

SUBROUTINE WRITE_ICGEM(outputFile, modelName, maxDegree, anm, r, GM)
  IMPLICIT NONE

  INTEGER                       :: n, m
  INTEGER, INTENT(IN)           :: maxDegree
  CHARACTER*130, INTENT(IN)     :: outputFile, modelName
  DOUBLE PRECISION, INTENT(IN)  :: anm(maxdegree+1, maxdegree+1), r, GM

  WRITE(*,*) outputFile
  OPEN(31,file=outputFile,STATUS='REPLACE')

  WRITE(31,'(A54)') 'begin_of_head ========================================'
  WRITE(31,'(A22,X,A29)') 'modelname             ',modelName
  WRITE(31,'(A36)') 'product_type           gravity_field'
  WRITE(31,'(A22,X,ES18.12)') 'earth_gravity_constant',GM
  WRITE(31,'(A22,X,ES18.12)') 'radius                ',r
  WRITE(31,'(A22,X,I3)') 'max_degree            ',maxDegree
  WRITE(31,'(A39)') 'norm                   fully_normalized'
  WRITE(31,'(A25)') 'errors                 no'
  WRITE(31,*) ''
  WRITE(31,'(A53)') 'key   L   M       C                   S              '
  WRITE(31,'(A54)') 'end_of_head =========================================='

  131 FORMAT(A4,I4,X,I4,X,eS19.12,X,eS19.12)
  DO n = 0, maxDegree
     DO m = 0, n
        IF (m .EQ. 0) THEN
           WRITE(31,131) 'gfc ',n, m,anm(n+1,m+1), 0.D0
        ELSE
           WRITE(31,131) 'gfc ',n, m,anm(n+1,m+1), anm(m, n+1)
        END IF
     END DO
  END DO
  CLOSE(31)

END SUBROUTINE WRITE_ICGEM
