C ---------------------------------------------------------------------
C     gconst.for
C ---------------------------------------------------------------------
!  1 Jul 2016 - Made the Earth's angular velocity input data
! 28 Jun 2016 - Added calculations of Earth's angular velocity
! 27 Jun 2016 - Changed calling sequence from (gravity) to (gravity,coriolis)
! 27 Jun 2016 - Added terms to the /GV/ common block for Coriolis terms
! 30 Apr 2010 - Correction
C 23 Mar 2010 - Changed subroutine name to gconst
C 23 Mar 2010 - Added ENTRY GRAV
C 23 Mar 2010 - Added ENTRY GRAVg
C 23 Mar 2010 - Added ENTRY GRAVa
C 23 Mar 2010 - Corrected some errors
C 23 Mar 2010 - 
C 23 Mar 2010 - 
C 10 Mar 2008 - corrected length errors in blank common
C     19 Jul 2006 - Separated igrav (no argument) and igravg (w. argument)
C     29 Apr 2006 - correction
C     26 January 2006 - 
C     25 January 2006 - Added perturbations.
C     23 January 2006
C     CONSTANT GRAVITY MODEL
C ---------------------------------------------------------------------
      SUBROUTINE gconst
C ---------------------------------------------------------------------

      IMPLICIT double precision (a-h,o-z)

C For gravity, the following variables represent first derivatives of the
C the three components (r, theta, and phi), NOT SECOND DERIVATIVES.
C Example: pthtm is the partial derivative of g_sub_th with respect to time.

      TYPE :: gravity_and_derivs_1
         REAL(KIND=8) :: gs,gsr,gsth,gsph
         REAL(KIND=8) :: ptm,pr,pth,pph
         REAL(KIND=8) :: prtm,prr,prth,prph
         REAL(KIND=8) :: pthtm,pthr,pthth,pthph
         REAL(KIND=8) :: pphtm,pphr,pphth,pphph
      END TYPE 
      TYPE (gravity_and_derivs_1) :: gravity
      TYPE (gravity_and_derivs_1) :: coriolis                      ! 27 Jun 2016

      double precision kr,kth,kph                                  ! 10 Mar 2008
      COMMON//r,th,ph,kr,kth,kph,rkvars(14),tpulse,cstep,drdt(20)  ! 10 Mar 2008

      PARAMETER (NWARSZ=1000)
      COMMON/WW/dum(10),MAXW,W(NWARSZ)
      DOUBLE PRECISION gvmodel,gvform,gvid

      EQUIVALENCE (W(1),earthr)

C     GRAVITY  600-624
      EQUIVALENCE (W(600),gvmodel),(W(601),gvform),(W(602),gvid)

C     COMMON DECK "GV" INSERTED HERE
      DOUBLE PRECISION modgv
      COMMON/GV/MODGV(4)
     1  ,GV  ,PGVTM  ,PGVR  ,PGVTH  ,PGVPH
     2  ,GVR ,PGVRTM ,PGVRR ,PGVRTH ,PGVRPH
     3  ,GVTH,PGVTHTM,PGVTHR,PGVTHTH,PGVTHPH
     4  ,GVPH,PGVPHTM,PGVPHR,PGVPHTH,PGVPHPH
     5  ,OMV  ,POMVTM  ,POMVR  ,POMVTH  ,POMVPH                    ! 27 Jun 2016
     6  ,OMVR ,POMVRTM ,POMVRR ,POMVRTH ,POMVRPH                   ! 27 Jun 2016
     7  ,OMVTH,POMVTHTM,POMVTHR,POMVTHTH,POMVTHPH                  ! 27 Jun 2016
     8  ,OMVPH,POMVPHTM,POMVPHR,POMVPHTH,POMVPHPH                  ! 27 Jun 2016

      REAL(KIND=8) :: gcnst
      EQUIVALENCE (W(603),gcnst),(W(604),omegaE)                   !  1 Jul 2016

C ---------------------------------------------------------------------
      ENTRY grav(gravity,coriolis)                   ! 23 Mar 2010 ! 27 Jun 2016
      ENTRY gravg(gravity,coriolis)                  ! 23 Mar 2010 ! 27 Jun 2016

!      earthr2 = earthr**2                          ! Do not need. ! 28 Jun 2016

      gravity%gs = gcnst
      gravity%gsr = -gcnst
c      gravity%gs = gcnst*earthr2/(r**2)
c      gravity%gsr = -gcnst*earthr2/(r**2)
      gravity%gsth = 0.d0
      gravity%gsph = 0.d0

      gravity%ptm = 0.d0
c      gravity%pr = -2*gcnst*earthr2/(r**3)
      gravity%pr =  0.d0
      gravity%pth = 0.d0
      gravity%pph = 0.d0

      gravity%prtm = 0.d0
c      gravity%prr = 2*gcnst*earthr2/(r**3)
      gravity%prr = 0.d0
      gravity%prth = 0.d0
      gravity%prph = 0.d0

      gravity%pthtm = 0.d0
      gravity%pthr = 0.d0
      gravity%pthth = 0.d0
      gravity%pthph = 0.d0

      gravity%pphtm = 0.d0
      gravity%pphr = 0.d0
      gravity%pphth = 0.d0
      gravity%pphph = 0.d0
  
      costh=dcos(th)                                               ! 28 Jun 2016
      sinth=dsin(th)                                               ! 28 Jun 2016
!      omegaE=7.2921158D-5                           ! 28 Jun 2016 !  1 Jul 2016
      coriolis%gs = omegaE                                         ! 28 Jun 2016
      coriolis%gsr = omegaE*costh                                  ! 28 Jun 2016
      coriolis%gsth = -omegaE*sinth                                ! 28 Jun 2016
      coriolis%gsph = 0.d0                                         ! 28 Jun 2016
        
      coriolis%ptm = 0.d0                                          ! 28 Jun 2016
      coriolis%pr =  0.d0                                          ! 28 Jun 2016
      coriolis%pth = 0.d0                                          ! 28 Jun 2016
      coriolis%pph = 0.d0                                          ! 28 Jun 2016

      coriolis%prtm = 0.d0                                         ! 28 Jun 2016
      coriolis%prr = 0.d0                                          ! 28 Jun 2016
      coriolis%prth = -omegaE*sinth                                ! 28 Jun 2016
      coriolis%prph = 0.d0                                         ! 28 Jun 2016

      coriolis%pthtm = 0.d0                                        ! 28 Jun 2016
      coriolis%pthr = 0.d0                                         ! 28 Jun 2016
      coriolis%pthth = -omegaE*costh                               ! 28 Jun 2016
      coriolis%pthph = 0.d0                                        ! 28 Jun 2016

      coriolis%pphtm = 0.d0                                        ! 28 Jun 2016
      coriolis%pphr = 0.d0                                         ! 28 Jun 2016
      coriolis%pphth = 0.d0                                        ! 28 Jun 2016
      coriolis%pphph = 0.d0                                        ! 28 Jun 2016

!      CALL clear(gv,20) ! 29 Apr 2006
      CALL clear(gv,40)                                          ! 27 Jun 2016
      gv = gravity%gs                                            ! 23 Mar 2010
      gvr = gravity%gsr                                          ! 23 Mar 2010
      gvth = gravity%gsth                                        ! 23 Mar 2010
      gvph = gravity%gsph                                        ! 23 Mar 2010

      omv = coriolis%gs                                            ! 28 Jun 2016
      omvr = coriolis%gsr                                          ! 28 Jun 2016
      omvth = coriolis%gsth                                        ! 28 Jun 2016
      pomvrth = coriolis%prth                                      ! 28 Jun 2016
      pomvthth = coriolis%pthth                                    ! 28 Jun 2016

c      print *, ' before pgrav'
C      CALL pgrav(gravity) ! 25 January 2006
      CALL PGRAVg(gravity,coriolis)                  ! 23 Mar 2010 ! 27 Jun 2016

c      print *, ' after pgrav'

      RETURN
C ---------------------------------------------------------------------
      ENTRY GRAVa                                         ! 23 Mar 2010
      CALL clear(gv,40)                                   ! 27 Jun 2016
      gv = gcnst                                          ! 23 Mar 2010
      gvr = -gcnst                                        ! 23 Mar 2010
      gvth = 0.d0                                         ! 23 Mar 2010
      gvph = 0.d0                                         ! 23 Mar 2010

c      print *, ' before PGRAva'
      CALL PGRAVa
c      print *, ' after PGRAVa'

      RETURN
C ---------------------------------------------------------------------
      END SUBROUTINE gconst
C ---------------------------------------------------------------------


C ---------------------------------------------------------------------
      SUBROUTINE igrav ! 19 Jul 2006
C ---------------------------------------------------------------------
C     This subroutine is called once for each new W array.

      IMPLICIT double precision (a-h,o-z)

      TYPE :: model
         CHARACTER :: model_name
         REAL(KIND=8) :: data_id
         CHARACTER :: pert_name
         REAL(KIND=8) :: pert_id
      END TYPE model
      TYPE (model) :: modgg

      PARAMETER (NWARSZ=1000)
      COMMON/WW/dum(10),MAXW,W(NWARSZ)
      DOUBLE PRECISION gvmodel,gvform,gvid

C     GRAVITY     600-624
      EQUIVALENCE (W(600),gvmodel),(W(601),gvform),(W(602),gvid)

C     COMMON DECK "GV" INSERTED HERE
      DOUBLE PRECISION modgv
      COMMON/GV/MODGV(4)
     1  ,GV  ,PGVTM  ,PGVR  ,PGVTH  ,PGVPH
     2  ,GVR ,PGVRTM ,PGVRR ,PGVRTH ,PGVRPH
     3  ,GVTH,PGVTHTM,PGVTHR,PGVTHTH,PGVTHPH
     4  ,GVPH,PGVPHTM,PGVPHR,PGVPHTH,PGVPHPH
     5  ,OMV  ,POMVTM  ,POMVR  ,POMVTH  ,POMVPH                    ! 27 Jun 2016
     6  ,OMVR ,POMVRTM ,POMVRR ,POMVRTH ,POMVRPH                   ! 27 Jun 2016
     7  ,OMVTH,POMVTHTM,POMVTHR,POMVTHTH,POMVTHPH                  ! 27 Jun 2016
     8  ,OMVPH,POMVPHTM,POMVPHR,POMVPHTH,POMVPHPH                  ! 27 Jun 2016

      DOUBLE PRECISION gcnst
      EQUIVALENCE (W(603),gcnst)

C     COMMON DECK "B13" INSERTED HERE
c      INTEGER gvmx,gvntbl,gvitbl,gvfrmtbl,idsgv(10)
c      DOUBLE PRECISION gvgp
c      COMMON/B13/gvmx,gvntbl(10),gvitbl(10),gvfrmtbl(10),gvgp(10)
c      EQUIVALENCE (gvgp,idsgv)
c      DATA gvmx/1/
c      DATA gvntbl/1,11,8*0/
c      DATA gvitbl/1,9*0/
c      DATA gvfrmtbl/1,9*0/

      IF (gvmodel /= 1.d0) THEN
           PRINT *,' Model check number mismatch: '
           STOP ' W(600) should be 1 if you are using model GCONST.'
      ENDIF

      modgv(1) = 'GCONST'
!      modgv(1) = 6hGCONST
      modgv(2) = gvid

      CALL ipgrav ! 19 Jul 2006

      RETURN

      ENTRY igravg(modgg) ! 19 Jul 2006
      IF (gvmodel /= 1.d0) THEN
           PRINT *,' Model check number mismatch: '
           STOP ' W(600) should be 1 if you are using model GCONST.'
      ENDIF

      modgg%model_name = 'GCONST'
      modgg%data_id = gvid
      modgg%pert_name = ''
      modgg%pert_id = 0.d0

!      modgv(1) = modgg%model_name                                ! 30 Apr 2010
!      modgv(2) = modgg%data_id                                   ! 30 Apr 2010

      modgv(1) = 'GCONST'                                         ! 30 Apr 2010
      modgv(2) = gvid                                             ! 30 Apr 2010

      CALL ipgravg(modgg) ! 25 January 2006

      RETURN
C ---------------------------------------------------------------------
      END SUBROUTINE igrav
C ---------------------------------------------------------------------


C ---------------------------------------------------------------------
      SUBROUTINE SETGRAV
C ---------------------------------------------------------------------
C
C     This subroutine is called once.
C
      implicit double precision (a-h,o-z)

C     COMMON DECK "B13" INSERTED HERE
      INTEGER DTMX,DTNTBL,DTITBL,DTFRMTB,IDSDT(10)
      COMMON/B13/DTMX,DTNTBL(10),DTITBL(10),DTFRMTB(10),DTGP(10)
      EQUIVALENCE (DTGP,IDSDT)
C
      INTEGER DWMX ,DWNTBL(10),DWITBL(10),DWFRMTB(10)
C
      DATA DWMX/1/
      DATA DWNTBL/1,11,8*0/
      DATA DWITBL/1,9*0/
      DATA DWFRMTB/1,9*0/
C
      DTMX=DWMX
      CALL IMOVE(DTNTBL,DWNTBL,10)
      CALL IMOVE(DTITBL,DWITBL,10)
      CALL IMOVE(DTFRMTB,DWFRMTB,10)
C
      CALL setpgrav

      RETURN
C
C ---------------------------------------------------------------------
      END SUBROUTINE SETGRAV
C ---------------------------------------------------------------------
