C ----------------------------------------------------------------------!-------
C     gnwnl.for
C ----------------------------------------------------------------------!-------
! 31 Jul 2017 - Moved factor of 2 from gnwnl to Coriolis2
! 19 Jul 2017 - CORRECTIONS TO HAMILTONIAN terms
!  7 Jul 2017 - Added calculation of f1 to Coriolis6
!  7 Jul 2017 - Added f1 to TYPE :: value_and_derivs_r_and_k
! 29 Dec 2016 - Tested for a negative value before taking a square root
! 29 Dec 2016 - Changed DSQRT to SQRT in two places
! 29 Dec 2016 - Changed some REAL (KIND=8) TO DOUBLE PRECISION
! 14 Sep 2016 - Finished adding term6 to all of the appropriate equations
! 13 Sep 2016 - Finished adding values of components for term6
! 12 Sep 2016 - Started adding values of components for term6
! 12 Sep 2016 - Added values of components for Gamma in coriolis6
!  5 Aug 2016 - Added common block /rink3/ to help with reflection of wave
!  4 Aug 2016 - Generalized calculation of kz2a
!  4 Aug 2016 - Generalized calculation of kay2
!  4 Aug 2016 - Calculated values of those components for term1 through term5
!  4 Aug 2016 - Added new components for type for term1 through term6
! 29 Jul 2016 - Added correct calculations of of derivatives of Coriolis terms
!               with respect to omega, kr, kth, and kph
! 29 Jul 2016 - Added new TYPE for term1 through term6
! 28 Jul 2016 - Corrected formulas for coriolis term additions to kay2 and kz2a
! 22 Jul 2016 - Added the 4 more Coriolis terms to kay2 and kz2a
! 22 Jul 2016 - Added call to coriolis2, coriolis3, coriolis4, coriolis5 
!               to add four more Coriolis terms to Hamiltonian
! 22 Jul 2016 - Added subroutines coriolis2, coriolis3, coriolis4, coriolis5, 
!               and coriolis6,  to calculate the rest of the Coriolis terms
!  1 Jul 2016 - Added factor605=w605 as a factor of the main Coriolis term
! 30 Jun 2016 - Added the main Coriolis term to kay2 and kz2a
! 29 Jun 2016 - Added call to coriolis1 to add main Coriolis term to Hamiltonian
! 29 Jun 2016 - Added subroutine coriolis1 to calculate the main Coriolis term
! 27 Jun 2016 - Added terms to the /GV/ common block for Coriolis terms
!  8 Apr 2011 - Corrected calculation of kz2a for case where kr2, kth2, or kph2 
!               might be negative
!  8 Apr 2011 - Added common block /rink2/
!  9 Dec 2010 - Added calculation of Hrel
!  9 Dec 2010 - Added calculation of errRel
!  9 Dec 2010 - Added common block /rink/
! 30 Apr 2010 - Slight change in setting model labeling
! 14 Apr 2010 - Moved some statemnts - should have no effect, added comments
C 25 Mar 2010 - Changed CALL density to CALL DENSITYg
C 25 Mar 2010 - Added ENTRY DISPER
C 25 Mar 2010 - Changed the subroutine name to GNWNL
C 13 Mar 2008 - Corrected length in common block /PCONST/
C 13 Mar 2008 - Corrected length in common block /RIN/
C 19 Jul 2006 - Added g to end of names with arguments in entry idisper
C 16 Mar 2006 - Added kzb to /RIN/
C 16 Mar 2006 - Changed kz2 to kz2a
C     14 Mar 2006 - Added kz2 to /RIN/
C     28 Feb 2006 - corrections
C     27 Feb 2006 - Converted those entry points to subroutines
C     27 Feb 2006 - Added entry points SETWND, SETABS, IWINDR, WINDR
C     27 Feb 2006 - correction
C     24 Feb 2006 - Calculate acoustic cutoff frequency in separate subroutine
C     24 Feb 2006 - Calculate Brunt-Vaisala frequency in separate subroutine
C      2 Feb 2006 - Added call to igrav
C     31 Jan 2006 - Added /B1/, /B2/, /B17/, and /B18/
C     26 Jan 2006 - Added call to setmw
C     25 Jan 2006 - Added call to setgrav
C     24 Jan 2006 - calling different subroutines for acoustics and 
C                   gravity
C     23 Jan 2006
C ----------------------------------------------------------------------!-------
C     DISPERSION RELATION FOR ACOUSTIC GRAVITY WAVES: NO WIND, NO LOSSES
! Barotropic case, assumes the gradients of pressure and density are vertical

      SUBROUTINE GNWNL ! (NO WINDS NO LOSSES)                     ! 25 Mar 2010

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

      TYPE :: value_and_derivs
         REAL(KIND=8) :: value
         REAL(KIND=8) :: pt,pr,pth,pph
      END TYPE 

      TYPE (value_and_derivs) :: kA2
      TYPE (value_and_derivs) :: Ns

      TYPE :: value_and_derivs_r_and_k                             ! 29 Jul 2016
         REAL(KIND=8) :: value
         REAL(KIND=8) :: pt,pr,pth,pph
         REAL(KIND=8) :: pom,pkr,pkth,pkph
         REAL(KIND=8) :: f0,f1,f2                                  !  7 Jul 2017
         REAL(KIND=8) :: fr,fth,fph                                !  4 Aug 2016
         REAL(KIND=8) :: frr,frth,frph,fthth,fthph,fphph           !  4 Aug 2016
      END TYPE 

      TYPE (value_and_derivs_r_and_k) :: term1       ! 29 Jun 2016 ! 29 Jul 2016
      TYPE (value_and_derivs_r_and_k) :: term2,term3,term4,term5   ! 29 Jul 2016
      TYPE (value_and_derivs_r_and_k) :: term6                     ! 13 Sep 2016

      TYPE :: values_and_derivs_1_2
         REAL(KIND=8) :: value
         REAL(KIND=8) :: ptm,pr,pth,pph
         REAL(KIND=8) :: ptmtm,prtm,pthtm,pphtm
         REAL(KIND=8) :: ptmr,prr,pthr,pphr
         REAL(KIND=8) :: ptmth,prth,pthth,pphth
         REAL(KIND=8) :: ptmph,prph,pthph,pphph
      END TYPE 

      TYPE (values_and_derivs_1_2) :: pressure
      TYPE (values_and_derivs_1_2) :: rho
      TYPE (values_and_derivs_1_2) :: potentialrho

      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    

c      REAL(KIND=8) :: kr,kth,kph
      double precision kr,kth,kph                                  ! 29 Dec 2016
      COMMON//R,TH,PH,KR,KTH,KPH,RKVARS(14),TPULSE,CSTEP,DRDT(20)

C     COMMON DECK "CC" INSERTED HERE
      DOUBLE PRECISION MODC
      COMMON/CC/MODC(4),CS,PCST,PCSR,PCSTH,PCSPH

C     COMMON DECK "CONST" INSERTED HERE
      COMMON/PCONST/CREF,RGAS,GAMMA 
     1 ,MO2,MN2,MO,MN,MH2O,MO3                        ! 9 Mar 2007 ! 13 Mar 2008
      DOUBLE PRECISION MO2,MN2,MO,MN,MH2O,MO3                      ! 13 Mar 2008
      COMMON/MCONST/PI,PIT2,PID2,DEGS,RAD,ALN10  

C     COMMON DECK "RK" INSERTED HERE
C     DEFINE SIZE REQUIRED FOR RAY STATE SAVE ARRAY
      PARAMETER (LRKAMS=87+2*100,NXRKMS=12+LRKAMS,MXEQPT=21)
      PARAMETER (NRKSAV=NXRKMS+MXEQPT-1)
      COMMON/RK/NEQS,STEP,MODE,E1MAX,E1MIN,E2MAX,E2MIN,FACT,RSTART

C     COMMON DECK "RINREAL" INSERTED HERE
      REAL(KIND=8) :: LPOLAR,LPOLRI,KPHK,KPHKI,KAY2,KAY2I 
     & , kz2a,kz2ai,kzb,kzbi  ! 16 Mar 2006
      CHARACTER(LEN=64) :: modrin
      LOGICAL :: space
      COMMON/RIN/modrin,rayname(2,3),type(3),space
      COMMON/RIN/OMEGMIN,OMEGMAX,KAY2,KAY2I 
     & , kz2a,kz2ai,kzb,kzbi                ! 16 Mar 2006
     1 ,H,HI,PHT,PHTI,PHR,PHRI,PHTH,PHTHI,PHPH,PHPHI
     2 ,PHOW,PHOWI,PHKR,PHKRI,PHKTH,PHKTI,PHKPH,PHKPI 
     3 ,KPHK,KPHKI,POLAR,POLARI,LPOLAR,LPOLRI,SGN 

      DOUBLE PRECISION kr2,kth2,kph2                               !  8 Apr 2011
      COMMON /rink2/ kr2,kth2,kph2                                 !  8 Apr 2011

      DOUBLE PRECISION azz,bzz,czz                                 !  5 Aug 2016
      COMMON /rink3/ azz,bzz,czz                                   !  5 Aug 2016
     &,czz0a,czz0b,czzf0,czzfth,czzfph,czzfthth,czzfphph,czzfthph  !  5 Aug 2016

      common /rink/ kaysq(4),kaysqb(4),kaysqa(4)                   !  9 Dec 2010
     & ,kaysqi(4),kaysqbi(4),kaysqai(4)                           
     & ,sol(4),soli(4),solb(4),solbi(4),sola(4),solai(4)          
     & ,alpha,dd                                                  
     & ,relerr(4),relerra(4),relerrb(4)                           
     & ,errRell(4),errRela(4),errRelb(4)                          
     & ,n3,n4,nbig,nsmall                                         
     & ,Hrel,errRel                                               
      complex*16 relerr,relerra,relerrb                            !  9 Dec 2010
      complex*16 errRell,errRela,errRelb                          
      double precision kaysq,kaysqb,kaysqa                        
     & ,kaysqi,kaysqbi,kaysqai                                    
      dimension dd(0:4)
!     & ,sol(4),soli(4),solb(4),solbi(4),sola(4),solai(4)          

C     COMMON DECK "UU" INSERTED HERE
      DOUBLE PRECISION MODU
      COMMON/UU/MODU(4)
     1  ,V  ,PVT  ,PVR  ,PVTH  ,PVPH
     2  ,VR ,PVRT ,PVRR ,PVRTH ,PVRPH
     3  ,VTH,PVTHT,PVTHR,PVTHTH,PVTHPH
     4  ,VPH,PVPHT,PVPHR,PVPHTH,PVPHPH

      COMMON/WW/dum(10),MAXW,W(1000)

      EQUIVALENCE (w(1),earthr)
      REAL(KIND=8) :: maxstp,maxerr,intyp,llat,llong
c      REAL(KIND=8) :: ks
      double precision ks                                          ! 29 Dec 2016
      REAL(KIND=8),DIMENSION(24) :: kvect
      REAL(KIND=8),DIMENSION(20) :: vset
c      REAL(KIND=8) :: ow,ows                                    ! 27 Feb 2006
      double precision ow,ows                                      ! 29 Dec 2016
      EQUIVALENCE(kvect,kay2),(vset,v),(ow,w(6))
     & ,(W(606),factor606)                                         ! 29 Jul 2016

      ENTRY disper                                                ! 25 Mar 2010

      space = .false.
      ks = kr**2 + kth**2 + kph**2

      ows = ow**2

c      print *, ' before speed '
      CALL speed
c      print *, ' after speed '
c      print *, ' before cutoff'
      CALL cutoff(kA2)                       ! 24 Feb 2006
c      print *, ' before Brunt '
      CALL Brunt(Ns)                         ! 24 Feb 2006

      call coriolis1(term1)                                        ! 29 Jun 2016
      call coriolis2(term2)                                        ! 22 Jul 2016
      call coriolis3(term3)                                        ! 22 Jul 2016
      call coriolis4(term4)                                        ! 22 Jul 2016
      call coriolis5(term5)                                        ! 22 Jul 2016
      call coriolis6(term6)                                        ! 13 Sep 2016

c      print *, ' ow = ', ow
c      print *, ' ows = ', ows
c      print *, ' cs in gnwnl = ', cs
c      print *, ' kr = ', kr
c      print *, ' kth = ', kth
c      print *, ' kph = ', kph
c      print *, ' ks = ', ks
c      print *, ' kA2%value = ', kA2%value
c      print *, ' Ns%value = ', Ns%value

C ----------------------------------------------------------------------!-------

      cee = (ows/cs)-kA2%value
     & +(term1%f0+term2%f0+term3%f0+term4%f0+term5%f0+term6%f0)/ows 

      aaa = (ows-Ns%value)*(kth**2+kph**2)/(ows*ks) + (kr**2/ks)
     & -(term1%f2+term2%f2+term3%f2+term4%f2+term5%f2+term6%f2)/ows

      bee = (term1%fr*kr+term1%fth*kth+term1%fph*kph
     & + term2%fr*kr+term2%fth*kth+term2%fph*kph
     & + term3%fr*kr+term3%fth*kth+term3%fph*kph
     & + term4%fr*kr+term4%fth*kth+term4%fph*kph
     & + term5%fr*kr+term5%fth*kth+term5%fph*kph
     & + term6%fr*kr+term6%fth*kth+term6%fph*kph)/ows

      rad = 0.d0
      rad2 = 0.d0                                                  ! 29 Dec 2016
c      if(ks.ne.0.d0) rad = dsqrt(cee*aaa/ks+bee**2/(4.d0*ks*ks))
      if(ks.ne.0.d0) rad2 = (cee*aaa/ks+bee**2/(4.d0*ks*ks))       ! 29 Dec 2016
      if(rad2.gt.0.d0) rad = sqrt(rad2)                            ! 29 Dec 2016
      dee=0.d0
      if(ks.ne.0.d0) dee = bee**2/(2.d0*ks)

!      kay2 = cee/aaa       
      kay2 = (cee*aaa+dee+bee*rad)/aaa**2                          !  4 Aug 2016

c      print *, 'aaa disper= ', aaa
c      print *, 'bee disper= ', bee
c      print *, 'cee disper= ', cee
c      print *, 'dee disper= ', dee
c      print *, 'rad disper= ', rad
      if(isnan(rad)) print *, 'aaa disper= ', aaa
      if(isnan(rad)) print *, 'bee disper= ', bee
      if(isnan(rad)) print *, 'cee disper= ', cee
      if(isnan(rad)) print *, ' ks = ', ks
      if(isnan(rad)) print *, ' rad2 = ', rad2
c      print *, 'kay2 disper= ', kay2
      if(kay2.lt.0.d0) print *, 'kay2 disper= ', kay2,  ' ks = ', ks
c     & ,'Coriolis term 5 = ',term5%value,term5%pom,term5%f0,term5%f2
     & , 'aaa disper= ', aaa, 'cee disper= ', cee

C ----------------------------------------------------------------------!-------

      CALL GRAVg(gravity,coriolis)  

      czz0a = (ows/cs)-kA2%value 
      czz0b =  -(kth**2+kph**2)*(ows-Ns%value)/ows 
      czzf0 =
     & (term1%f0+term2%f0+term3%f0+term4%f0+term5%f0+term6%f0
     & )/ows                                  
      czzfth =
     & (term1%fth+term2%fth+term3%fth+term4%fth+term5%fth+term6%fth)*kth
     & /ows                                  
      czzfph =
     & (term1%fph+term2%fph+term3%fph+term4%fph+term5%fph+term6%fph)*kph
     & /ows                                  
      czzfthth =
     & (term1%fthth+term2%fthth+term3%fthth+term4%fthth+term5%fthth
     & +term6%fthth)
     &  *kth*kth
     & /ows                                  
      czzfphph =
     & (term1%fphph+term2%fphph+term3%fphph+term4%fphph+term5%fphph
     & +term5%fphph)
     &  *kph*kph
     & /ows                                  
      czzfthph =
     & (term1%fthph+term2%fthph+term3%fthph+term4%fthph+term5%fthph
     & +term5%fthph)
     &  *kth*kph                                                   ! 31 Jul 2017
     & /ows                                  

      czz = (ows/cs)-kA2%value -(kth**2+kph**2)*(ows-Ns%value)/ows  
     & +(term1%f0+term2%f0+term3%f0+term4%f0+term5%f0+term6%f0
     & +(term1%fth+term2%fth+term3%fth+term4%fth+term5%fth
     &+term6%fth)*kth
     & +(term1%fph+term2%fph+term3%fph+term4%fph+term5%fph
     & +term6%fph)*kph
     & +(term1%fthth+term2%fthth+term3%fthth+term4%fthth+term5%fthth
     &+term6%fthth)
     &  *kth*kth
     & +(term1%fphph+term2%fphph+term3%fphph+term4%fphph+term5%fphph
     & +term6%fphph)
     &  *kph*kph
     & +(term1%fthph+term2%fthph+term3%fthph+term4%fthph+term5%fthph
     & +term6%fthph)
     &  *kth*kph                                                   ! 31 Jul 2017
     & )/ows                                  

      azz=1.d0-(term1%frr+term2%frr+term3%frr+term4%frr+term5%frr
     & +term6%frr)/ows   

      bzz=(term1%fr+term2%fr+term3%fr+term4%fr+term5%fr
     & + (term1%frth+term2%frth+term3%frth+term4%frth+term5%frth   ! 31 Jul 2017
     & +term6%frth)
     & *kth
     & + (term1%frph+term2%frph+term3%frph+term4%frph+term5%frph   ! 31 Jul 2017
     & +term6%frph)
     & *kph
     & )/(2.d0*ows)

c      rad = dsqrt(azz*czz+bzz*bzz)
      rad2 = azz*czz+bzz*bzz                                       ! 29 Dec 2016
      rad = 0.d0                                                   ! 29 Dec 2016
      if(rad2.gt.0.d0) rad = dsqrt(rad2)                           ! 29 Dec 2016

!      kz2a=czz/azz
      kz2a=(azz*czz+2.d0*bzz*bzz+2.d0*bzz*rad)/azz**2              !  4 Aug 2016

! When only the ratios of kr, kth, and kph are correct, the following is the 
! correct formula.
!  7 Apr 2011      kz2a = kay2*kr**2/ks                            !  7 Apr 2011
!      kz2a = kay2*kr**2/ks                                         ! 29 Jul 2016
      kz2ai = 0.d0 ! 16 Mar 2006
      kzb = 0.d0 ! 16 Mar 2006
      kzbi = 0.d0 ! 16 Mar 2006

      H=(kth**2+kph**2)*(Ns%value-ows) - ows*(kr**2+kA2%value-(ows/cs))
     & + term1%value                                               ! 29 Jun 2016
     & + term2%value + term3%value + term4%value + term5%value     ! 22 Jul 2016
     & + term6%value                                               ! 13 Sep 2016
      H1=(kth**2+kph**2)*(Ns%value)                                !  9 Dec 2010
      H2=(kth**2+kph**2)*(-ows)                                    !  9 Dec 2010
      H3=-ows*(kr**2)                                !  9 Dec 2010 ! 19 Jul 2017
      H4=-ows*(kA2%value)                            !  9 Dec 2010 ! 19 Jul 2017
      H5=-ows*(-(ows/cs))                            !  9 Dec 2010 ! 19 Jul 2017
      H6 = term1%value                                             ! 29 Jun 2016
      H7 = term2%value                                             ! 22 Jul 2016
      H8 = term3%value                                             ! 22 Jul 2016
      H9 = term4%value                                             ! 22 Jul 2016
      H10 = term5%value                                            ! 22 Jul 2016
      H11 = term6%value                                            ! 13 Sep 2016
      Hrel=H/max(abs(H1),abs(H2),abs(H3),abs(H4),abs(H5),abs(H6)   !  9 Dec 2010
     & ,abs(H7),abs(H8),abs(H9),abs(H10),abs(H11))  ! 22 Jul 2016  ! 13 Sep 2016
      height = r - earthr
c      print *,'height = ',height
c      print *,'Ns%value = ',Ns%value,'Ns%pr = ',Ns%pr
c      print *,'Ns%value = ',Ns%value
c      print *,'Ns%pt = ',Ns%pt
c      print *,'Ns%pr = ',Ns%pr
c      print *,'Ns%pth = ',Ns%pth
c      print *,'Ns%pph = ',Ns%pph

c      print *,'kA2%pt = ',kA2%pt
c      print *,'kA2%pr = ',kA2%pr
c      print *,'kA2%pth = ',kA2%pth
c      print *,'kA2%pph = ',kA2%pph

      pHow = -2.d0*(kth**2+kph**2)*ow
     &       -2.d0*(kr**2+kA2%value)*ow+4*ow*ows/cs
     &       +term1%pom+term2%pom+term3%pom+term4%pom+term5%pom    ! 29 Jul 2016
     &       +term6%pom                                            ! 13 Sep 2016
      pHkr = -2.d0*ows*kr
     &       +term1%pkr+term2%pkr+term3%pkr+term4%pkr+term5%pkr    ! 29 Jul 2016
     &       +term6%pkr                                            ! 13 Sep 2016
      pHkth = 2.d0*(Ns%value-ows)*kth 
     &     +term1%pkth+term2%pkth+term3%pkth+term4%pkth+term5%pkth ! 29 Jul 2016
     &     +term6%pkth                                             ! 13 Sep 2016
      pHkph = 2.d0*(Ns%value-ows)*kph 
     &     +term1%pkph+term2%pkph+term3%pkph+term4%pkph+term5%pkph ! 29 Jul 2016
     &     +term6%pkph                                             ! 13 Sep 2016
      kpHk = -2.d0*ows*ks + 2.d0*Ns%value*(kth**2+kph**2)
     &     +kr*(term1%pkr+term2%pkr+term3%pkr+term4%pkr+term5%pkr  ! 29 Jul 2016
     &         +term6%pkr)                                         ! 13 Sep 2016
     &+kth*(term1%pkth+term2%pkth+term3%pkth+term4%pkth+term5%pkth ! 29 Jul 2016
     &     +term6%pkth)                                            ! 13 Sep 2016
     &+kph*(term1%pkph+term2%pkph+term3%pkph+term4%pkph+term5%pkph ! 29 Jul 2016
     &     +term6%pkph)                                            ! 13 Sep 2016

c      print *,'kr = ',kr
c      print *,'kth = ',kth
c      print *,'kph = ',kph
c      print *,'pHkr = ',pHkr
cc      print *,'pHkth = ',pHkth = 0.00000000
cc      print *,'pHkph = ',pHkph = 0.00000000

      pHt = (kth**2+kph**2)*Ns%pt - ((ows/cs)**2)*pcst - ows*kA2%pt
     & +term1%pt                                                   ! 29 Jun 2016
     & +term2%pt+term3%pt+term4%pt+term5%pt+term6%pt ! 22 Jul 2016 ! 13 Sep 2016
      pHr = (kth**2+kph**2)*Ns%pr - ((ows/cs)**2)*pcsr - ows*kA2%pr
     & +term1%pr                                                   ! 29 Jun 2016
     & +term2%pr+term3%pr+term4%pr+term5%pr+term6%pr ! 22 Jul 2016 ! 13 Sep 2016
      pHth = (kth**2+kph**2)*Ns%pth - ((ows/cs)**2)*pcsth - ows*kA2%pth
     & +term1%pth                                                  ! 29 Jun 2016
     & +term2%pth+term3%pth+term4%pth+term5%pth+term6%pth!22Jul2016! 13 Sep 2016
      pHph = (kth**2+kph**2)*Ns%pph - ((ows/cs)**2)*pcsph - ows*kA2%pph
     & +term1%pph                                                  ! 29 Jun 2016
     & +term2%pph+term3%pph+term4%pph+term5%pph+term6%pph!22Jul 2016 13 Sep 2016

c      print *,'H = ',H
c      print *,'pHt = ',pHt
c      print *,'pHr = ',pHr
c      print *,'pHth = ',pHth
c      print *,'pHph = ',pHph
cc      print *,'drdt(3) = ',drdt(3) = 0.00000000

      errRel=1.d10                                                 !  9 Dec 2010
      IF (KAY2.NE.0.d0) errRel =ks/KAY2-1.d0                       !  9 Dec 2010

      ENTRY WINDR                     ! 27 Feb 2006
      RETURN
C ----------------------------------------------------------------------!-------
      END SUBROUTINE GNWNL 
C ----------------------------------------------------------------------!-------

C ----------------------------------------------------------------------!-------
      SUBROUTINE idisper 
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 
      TYPE (model) :: modpg
      TYPE (model) :: moddg
      TYPE (model) :: modgg ! 2 Feb 2006

C     COMMON DECK "PP" INSERTED HERE
      DOUBLE PRECISION MODP
      COMMON/PP/MODP(4),P,PPT,PPR,PPTH,PPPH

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

C     COMMON DECK "DS" INSERTED HERE                            ! 27 Mar 2007
      DOUBLE PRECISION ds,modds
      COMMON/DS/MODDS(4),ds,pdst,pdsr,pdsth,pdsph

      CALL ispeed
      CALL irecvr   
      CALL itopog   
      CALL isurfac
      CALL ipresg(modpg) ! 24 January 2006
!      modp(1) = modpg%model_name                                  ! 30 Apr 2010
!      modp(2) = modpg%data_id                                     ! 30 Apr 2010
!      MODP(3) = modpg%pert_name                                   ! 30 Apr 2010
!      MODP(4) = modpg%pert_id                                     ! 30 Apr 2010

      CALL igravg(modgg) ! 2 Feb 2006 ! 19 Jul 2006
!      modgv(1) = modgg%model_name                                 ! 30 Apr 2010
!      modgv(2) = modgg%data_id                                    ! 30 Apr 2010
!      MODGV(3) = modgg%pert_name                                  ! 30 Apr 2010
!      MODGV(4) = modgg%pert_id                                    ! 30 Apr 2010

      CALL idensityg(moddg) ! 19 Jul 2006
!      modds(1) = moddg%model_name                                 ! 30 Apr 2010
!      modds(2) = moddg%data_id                                    ! 30 Apr 2010
!      MODDS(3) = moddg%pert_name                                  ! 30 Apr 2010
!      MODDS(4) = moddg%pert_id                                    ! 30 Apr 2010

      ENTRY IWINDR                     ! 27 Feb 2006
      RETURN
C ----------------------------------------------------------------------!-------
      END SUBROUTINE idisper
C ----------------------------------------------------------------------!-------

C ----------------------------------------------------------------------!-------
      SUBROUTINE setdsp
C ----------------------------------------------------------------------!-------
C     This subroutine is called once.

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

      CHARACTER(LEN=6) :: dispm
      CHARACTER(LEN=64) :: modrin
      CHARACTER(LEN=5) :: dlspm
      REAL(KIND=8),DIMENSION(3) :: tbpe
      LOGICAL :: space

      COMMON/RINPL/dispm
      COMMON/RIN/modrin,rayname(2,3),type(3),space,stuff(33)       ! 13 Mar 2008

      DATA dlspm /'GNWNL'/
      DATA tbpe /3*'0'/ 

      dispm = dlspm
      CALL rmove(type,tbpe,3)
      modrin = 'ACOUSTIC GRAVITY WAVE --- NO WIND --- NO LOSSES'

      CALL setspd 
      CALL setrcv 
      CALL settop 
      CALL setsur
      CALL setpres
      CALL setdensity
      CALL setgrav
      CALL setmw
      CALL setwnd                                               ! 28 Feb 2006
      CALL setabs

      RETURN
C ----------------------------------------------------------------------!-------
      END SUBROUTINE setdsp
C ----------------------------------------------------------------------!-------

C ----------------------------------------------------------------------!-------
      SUBROUTINE SETWND                     ! 27 Feb 2006
C ----------------------------------------------------------------------!-------

C     This subroutine is called once.

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

C     COMMON DECK "B1" INSERTED HERE ! 31 Jan 2006
      INTEGER UMX,UNTBL,UITBL,UFRMTBL,IDSU(10) ! 31 Jan 2006
      COMMON/B1/UMX,UNTBL(10),UITBL(10),UFRMTBL(10),UGP(10)      ! 31 Jan 2006
      EQUIVALENCE (UGP,IDSU) ! 31 Jan 2006
      INTEGER UPX,UQTBL(10),ULTBL(10),UIRMTBL(10)                ! 31 Jan 2006
C
C     COMMON DECK "B2" INSERTED HERE ! 31 Jan 2006
C     DELTA WIND        125-149                                  ! 31 Jan 2006
      INTEGER DUMX,DUNTBL,DUITBL,DUFRMTB,IDSDU(20)               ! 31 Jan 2006
      COMMON/B2/DUMX,DUNTBL(10),DUITBL(10),DUFRMTB(10),DUGP(10)  ! 31 Jan 2006
      EQUIVALENCE (DUGP,IDSDU)                                   ! 31 Jan 2006
      INTEGER DXMX    ,DXNTBL(10),DXITBL(10),DXFRMTB(10)         ! 31 Jan 2006
C                            
      DATA UPX/1/ ! 31 Jan 2006
      DATA UQTBL/1,11,8*0/ ! 31 Jan 2006
      DATA ULTBL/1,9*0/ ! 31 Jan 2006
      DATA UIRMTBL/1,9*0/ ! 31 Jan 2006

      DATA DXMX/1/            ! 31 Jan 2006
      DATA DXNTBL/1,11,8*0/   ! 31 Jan 2006
      DATA DXITBL/1,9*0/      ! 31 Jan 2006
      DATA DXFRMTB/1,9*0/     ! 31 Jan 2006

      UMX=UPX ! 31 Jan 2006
      CALL IMOVE(UNTBL,UQTBL,10)     ! 31 Jan 2006
      CALL IMOVE(UITBL,ULTBL,10)     ! 31 Jan 2006
      CALL IMOVE(UFRMTBL,UIRMTBL,10) ! 31 Jan 2006

      DUMX=DXMX                      ! 31 Jan 2006
      CALL IMOVE(DUNTBL,DXNTBL,10)   ! 31 Jan 2006
      CALL IMOVE(DUITBL,DXITBL,10)   ! 31 Jan 2006
      CALL IMOVE(DUFRMTB,DXFRMTB,10) ! 31 Jan 2006

      RETURN
C ----------------------------------------------------------------------!-------
      END SUBROUTINE setwnd ! 28 Feb 2006
C ----------------------------------------------------------------------!-------

C ----------------------------------------------------------------------!-------
      SUBROUTINE SETABS                     ! 27 Feb 2006
C ----------------------------------------------------------------------!-------

C     This subroutine is called once.

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

C     COMMON DECK "CB17" INSERTED HERE                           ! 31 Jan 2006
      INTEGER VMX,VNTBL,VITBL,VFRMTBL,IDSV(10)                   ! 31 Jan 2006
      COMMON/B17/VMX,VNTBL(10),VITBL(10),VFRMTBL(10),VGP(53)     ! 31 Jan 2006
      EQUIVALENCE (VGP,IDSV),(ANV,VGP(11))                       ! 31 Jan 2006
      INTEGER VPX    ,VQTBL(10),VLTBL(10),VIRMTBL(10)            ! 31 Jan 2006
C                         
C     COMMON DECK "CB18" INSERTED HERE ! 31 Jan 2006
      INTEGER DVMX,DVNTBL,DVITBL,DVFRMTB,IDSDV(10)               ! 31 Jan 2006
      COMMON/B18/DVMX,DVNTBL(10),DVITBL(10),DVFRMTB(10),DVGP(11) ! 31 Jan 2006
      EQUIVALENCE (DVGP,IDSDV),(ANDV,DVGP(11))                   ! 31 Jan 2006
      INTEGER DYMX    ,DYNTBL(10),DYITBL(10),DYFRMTB(10)         ! 31 Jan 2006
C                           
      DATA VPX/1/          ! 31 Jan 2006
      DATA VQTBL/1,11,8*0/ ! 31 Jan 2006
      DATA VLTBL/1,9*0/    ! 31 Jan 2006
      DATA VIRMTBL/1,9*0/  ! 31 Jan 2006

      DATA DYMX/1/           ! 31 Jan 2006
      DATA DYNTBL/1,11,8*0/  ! 31 Jan 2006
      DATA DYITBL/1,9*0/     ! 31 Jan 2006
      DATA DYFRMTB/1,9*0/    ! 31 Jan 2006

      VMX=VPX                        ! 31 Jan 2006
      CALL IMOVE(VNTBL,VQTBL,10)     ! 31 Jan 2006
      CALL IMOVE(VITBL,VLTBL,10)     ! 31 Jan 2006
      CALL IMOVE(VFRMTBL,VIRMTBL,10) ! 31 Jan 2006

      DVMX=DYMX                      ! 31 Jan 2006
      CALL IMOVE(DVNTBL,DYNTBL,10)   ! 31 Jan 2006
      CALL IMOVE(DVITBL,DYITBL,10)   ! 31 Jan 2006
      CALL IMOVE(DVFRMTB,DYFRMTB,10) ! 31 Jan 2006

      RETURN
C ----------------------------------------------------------------------!-------
      END SUBROUTINE setabs
C ----------------------------------------------------------------------!-------

C ----------------------------------------------------------------------!-------
      SUBROUTINE cutoff(kA2)
C ----------------------------------------------------------------------!-------
C The acoustic cutoff frequency equals C times kA, where C is sound speed
! The formulas use only the vertical component of kA

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

      TYPE :: value_and_derivs
         REAL(KIND=8) :: value
         REAL(KIND=8) :: pt,pr,pth,pph
      END TYPE 

      TYPE (value_and_derivs) :: kA2

      TYPE :: values_and_derivs_1_2
         REAL(KIND=8) :: value
         REAL(KIND=8) :: ptm,pr,pth,pph
         REAL(KIND=8) :: ptmtm,prtm,pthtm,pphtm
         REAL(KIND=8) :: ptmr,prr,pthr,pphr
         REAL(KIND=8) :: ptmth,prth,pthth,pphth
         REAL(KIND=8) :: ptmph,prph,pthph,pphph
      END TYPE 

      TYPE (values_and_derivs_1_2) :: rho
      TYPE (values_and_derivs_1_2) :: potentialrho

c      print *, ' before density'
      CALL DENSITYg(rho,potentialrho)                              ! 25 Mar 2010
c      print *, ' after density'

      kA2%value = (rho%pr/(2.d0*rho%value))**2

      kA2%pt = (rho%pr*rho%prtm)/(2.d0*rho%value**2) -
     & ((rho%pr**2)*rho%ptm)/(2.d0*rho%value**3)
      kA2%pr = (rho%pr*rho%prr)/(2.d0*rho%value**2) -
     & ((rho%pr**3))/(2.d0*rho%value**3)
      kA2%pth = (rho%pr*rho%prth)/(2.d0*rho%value**2) -
     & ((rho%pr**2)*rho%pth)/(2.d0*rho%value**3)
      kA2%pph = (rho%pr*rho%prph)/(2.d0*rho%value**2) - 
     & ((rho%pr**2)*rho%pph)/(2.d0*rho%value**3)
c      print*,'ka2 cutoff = ',ka2%value
c      print*,'rho = ',rho%value
c      print*,'rho%pr = ',rho%pr

      RETURN
C ----------------------------------------------------------------------!-------
      END SUBROUTINE cutoff
C ----------------------------------------------------------------------!-------

C ----------------------------------------------------------------------!-------
      SUBROUTINE Brunt(Ns)
C ----------------------------------------------------------------------!-------

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

      REAL(KIND=8) :: kr,kth,kph
      COMMON//R,TH,PH,KR,KTH,KPH,RKVARS(14),TPULSE,CSTEP,DRDT(20)

      TYPE :: value_and_derivs
         REAL(KIND=8) :: value
         REAL(KIND=8) :: pt,pr,pth,pph
      END TYPE 

      TYPE (value_and_derivs) :: Ns

      TYPE :: values_and_derivs_1_2
         REAL(KIND=8) :: value
         REAL(KIND=8) :: ptm,pr,pth,pph
         REAL(KIND=8) :: ptmtm,prtm,pthtm,pphtm
         REAL(KIND=8) :: ptmr,prr,pthr,pphr
         REAL(KIND=8) :: ptmth,prth,pthth,pphth
         REAL(KIND=8) :: ptmph,prph,pthph,pphph
      END TYPE 

      TYPE (values_and_derivs_1_2) :: pressure
      TYPE (values_and_derivs_1_2) :: rho
      TYPE (values_and_derivs_1_2) :: potentialrho

c      print*,'Beginning of Brunt'

      CALL presg(pressure) 
      CALL DENSITYg(rho,potentialrho)                             ! 25 Mar 2010

c      print *,'potentialrho%pr = ',potentialrho%pr
c      print *,'potentialrho%pth = ',potentialrho%pth
c      print *,'potentialrho%pph = ',potentialrho%pph
c      print *,'pressure%pr = ',pressure%pr
c      print *,'pressure%pth = ',pressure%pth
c      print *,'pressure%pph = ',pressure%pph
c      print *,'rho%value = ',rho%value
c      print *,'th = ',th

      Ns%value = (potentialrho%pr*pressure%pr +
     & (1.d0/(r**2))*(potentialrho%pth*pressure%pth) +
     & ((potentialrho%pph*pressure%pph)/((r*dsin(th))**2)))/
     & (rho%value**2)
c      print *, 'Ns%value = ', Ns%value

      Ns%pt = -2.d0*Ns%value*rho%ptm/rho%value +
     & ((potentialrho%prtm*pressure%pr) +
     & (potentialrho%pr*pressure%prtm) +
     & ((potentialrho%ptmth*pressure%pth)/(r**2)) +
     & ((potentialrho%pth*pressure%pthtm)/(r**2)) +
     & ((potentialrho%pphtm*pressure%pph +
     & potentialrho%pph*pressure%pphtm)/((r*dsin(th))**2)))/
     & (rho%value**2)

      Ns%pr = -2.d0*Ns%value*rho%pr/rho%value +
     & (((potentialrho%prr*pressure%pr) +
     & (potentialrho%pr*pressure%prr) -
     & ((2.d0*potentialrho%pth*pressure%pth)/(r**3)) +
     & ((potentialrho%prth*pressure%pth)/(r**2)) +
     & ((potentialrho%pth*pressure%pthr)/(r**2)) -
     & ((2.d0*potentialrho%pph*pressure%pph)/(r**3*(dsin(th))**2)) +
     & ((potentialrho%pphr*pressure%pph +
     & potentialrho%pph*pressure%pphr)/((r*dsin(th))**2))))/
     & (rho%value**2)

      Ns%pth = -2.d0*Ns%value*rho%pth/rho%value +
     & ((potentialrho%prth*pressure%pr) +
     & (potentialrho%pr*pressure%prth) +
     & ((potentialrho%pthth*pressure%pth)/(r**2)) +
     & ((potentialrho%pth*pressure%pthth)/(r**2)) +
     & ((potentialrho%pphth*pressure%pph +
     & potentialrho%pph*pressure%pphth)/((r*dsin(th))**2)) -
     & (2.d0*dcos(th)*(potentialrho%pphth*pressure%pph +
     & potentialrho%pph*pressure%pphth))/
     & (r**2*(dsin(th))**3))/(rho%value**2)

      Ns%pph = -2.d0*Ns%value*rho%pph/rho%value +
     & ((potentialrho%prph*pressure%pr) +
     & (potentialrho%pr*pressure%prph) +
     & ((potentialrho%pthph*pressure%pth)/(r**2)) +
     & ((potentialrho%pth*pressure%pthph)/(r**2)) +
     & ((potentialrho%pphph*pressure%pph +
     & potentialrho%pph*pressure%pphph)/((r*dsin(th))**2)))/
     & (rho%value**2)

      RETURN
C ----------------------------------------------------------------------!-------
      END SUBROUTINE Brunt
C ----------------------------------------------------------------------!-------

C ----------------------------------------------------------------------!-------
      SUBROUTINE coriolis1(term1)                                  ! 29 Jun 2016
C ----------------------------------------------------------------------!-------
!  1 Jul 2016 - Added factor605=w605 as a factor of the main Coriolis term
C Term1 equals 4 times the square of the dot product of kA and OmegaE, where
!     OmegaE is the Earth's angular velocity in radians per second.

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

      TYPE :: value_and_derivs_r_and_k                             ! 29 Jul 2016
         REAL(KIND=8) :: value
         REAL(KIND=8) :: pt,pr,pth,pph
         REAL(KIND=8) :: pom,pkr,pkth,pkph
         REAL(KIND=8) :: f0,f1,f2                                  !  7 Jul 2017
         REAL(KIND=8) :: fr,fth,fph                                !  4 Aug 2016
         REAL(KIND=8) :: frr,frth,frph,fthth,fthph,fphph           !  4 Aug 2016
      END TYPE 

      TYPE (value_and_derivs_r_and_k) :: term1       ! 29 Jun 2016 ! 29 Jul 2016

      TYPE :: values_and_derivs_1_2
         REAL(KIND=8) :: value
         REAL(KIND=8) :: ptm,pr,pth,pph
         REAL(KIND=8) :: ptmtm,prtm,pthtm,pphtm
         REAL(KIND=8) :: ptmr,prr,pthr,pphr
         REAL(KIND=8) :: ptmth,prth,pthth,pphth
         REAL(KIND=8) :: ptmph,prph,pthph,pphph
      END TYPE 

      TYPE (values_and_derivs_1_2) :: rho
      TYPE (values_and_derivs_1_2) :: potentialrho

      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    
      TYPE (gravity_and_derivs_1) :: kA   

      double precision kAdotOmega

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

      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)
     & ,(W(605),factor605)                                         !  1 Jul 2016

c      print *, ' before density'
      CALL DENSITYg(rho,potentialrho) 
c      print *, ' after density'

      sinth = dsin(th)
      costh = dcos(th)

      kA%gsr = (rho%pr/(2.d0*rho%value))  

      kA%prtm = (rho%prtm)/(2.d0*rho%value) -
     & ((rho%pr)*rho%ptm)/(2.d0*rho%value**2)
      kA%prr = (rho%prr)/(2.d0*rho%value) -
     & ((rho%pr**2))/(2.d0*rho%value**2)
      kA%prth = (rho%prth)/(2.d0*rho%value) -
     & ((rho%pr)*rho%pth)/(2.d0*rho%value**2)
      kA%prph = (rho%prph)/(2.d0*rho%value) - 
     & ((rho%pr)*rho%pph)/(2.d0*rho%value**2)

      kA%gsth = (rho%pth/(2.d0*r*rho%value)) 

      kA%pthtm = (rho%pthtm)/(2.d0*r*rho%value) -
     & ((rho%pth)*rho%ptm)/(2.d0*r*rho%value**2)
      kA%pthr = (rho%pthr)/(2.d0*r*rho%value) -
     & ((rho%pth*rho%pr))/(2.d0*r*rho%value**2) -
     & ((rho%pth))/(2.d0*r**2*rho%value)
      kA%pthth = (rho%pthth)/(2.d0*r*rho%value) -
     & ((rho%pth)*rho%pth)/(2.d0*r*rho%value**2)
      kA%pthph = (rho%pthph)/(2.d0*r*rho%value) - 
     & ((rho%pth)*rho%pph)/(2.d0*r*rho%value**2)

      kA%gsph = (rho%pph/(2.d0*r*sinth*rho%value))

      kA%pphtm = (rho%pphtm)/(2.d0*r*sinth*rho%value) -
     & ((rho%pph)*rho%ptm)/(2.d0*r*sinth*rho%value**2)
      kA%pphr = (rho%pphr)/(2.d0*r*sinth*rho%value) -
     & ((rho%pph*rho%pr))/(2.d0*r*sinth*rho%value**2) -
     & ((rho%pph))/(2.d0*r**2*sinth*rho%value)
      kA%pphth = (rho%pphth)/(2.d0*r*sinth*rho%value) -
     & ((rho%pph)*rho%pth)/(2.d0*r*sinth*rho%value**2) -
     & ((rho%pph)*costh)/(2.d0*r*sinth**2*rho%value)
      kA%pphph = (rho%pphph)/(2.d0*r*sinth*rho%value) - 
     & ((rho%pph)*rho%pph)/(2.d0*r*sinth*rho%value**2)

! The next 5 statements use only the vertical component of kA, 
! but that is OK, because they are not actually needed for Coriolis term 1

      kA%gs = kA%gsr 
      kA%ptm = kA%prtm 
      kA%pr = kA%prr 
      kA%pth = kA%prth 
      kA%pph = kA%prph 

      CALL GRAVg(gravity,coriolis)  

      kAdotOmega = kA%gsr*coriolis%gsr + kA%gsth*coriolis%gsth
     &  + kA%gsph*coriolis%gsph 

      term1%f0 = factor605*4.d0*kAdotOmega**2                      !  4 Aug 2016
      term1%f1=0.d0                                                !  7 Jul 2017
      term1%f2=0.d0
      term1%fr=0.d0
      term1%fth=0.d0
      term1%fph=0.d0
      term1%frr=0.d0
      term1%frth=0.d0
      term1%frph=0.d0
      term1%fthth=0.d0
      term1%fthph=0.d0
      term1%fphph=0.d0

      term1%value = factor605*4.d0*kAdotOmega**2                   !  1 Jul 2016
      term1%pt = factor605*8.d0*kAdotOmega*(                       !  1 Jul 2016
     &   kA%prtm*coriolis%gsr + kA%gsr*coriolis%prtm 
     & + kA%pthtm*coriolis%gsth + kA%gsth*coriolis%pthtm 
     & + kA%pphtm*coriolis%gsph + kA%gsph*coriolis%pphtm)
      term1%pr = factor605*8.d0*kAdotOmega*(                       !  1 Jul 2016
     &   kA%prr*coriolis%gsr + kA%gsr*coriolis%prr 
     & + kA%pthr*coriolis%gsth + kA%gsth*coriolis%pthr 
     & + kA%pphr*coriolis%gsph + kA%gsph*coriolis%pphr)
      term1%pth = factor605*8.d0*kAdotOmega*(                      !  1 Jul 2016
     &   kA%prth*coriolis%gsr + kA%gsr*coriolis%prth 
     & + kA%pthth*coriolis%gsth + kA%gsth*coriolis%pthth 
     & + kA%pphth*coriolis%gsph + kA%gsph*coriolis%pphth)
      term1%pph = factor605*8.d0*kAdotOmega*(                      !  1 Jul 2016
     &   kA%prph*coriolis%gsr + kA%gsr*coriolis%prph 
     & + kA%pthph*coriolis%gsth + kA%gsth*coriolis%pthph 
     & + kA%pphph*coriolis%gsph + kA%gsph*coriolis%pphph)

      term1%pom=0.d0                                               ! 29 Jul 2016
      term1%pkr=0.d0                                               ! 29 Jul 2016
      term1%pkth=0.d0                                              ! 29 Jul 2016
      term1%pkph=0.d0                                              ! 29 Jul 2016

c      print*,'ka = ',kA%gsr
c      print*,'Coriolis term 1 = ',term1%value

      RETURN
C ----------------------------------------------------------------------!-------
      END SUBROUTINE coriolis1
C ----------------------------------------------------------------------!-------


C ----------------------------------------------------------------------!-------
      SUBROUTINE coriolis2(term2)                                  ! 22 Jul 2016
C ----------------------------------------------------------------------!-------
C Term2 equals 4 times the square of the dot product of k and OmegaE, where
!     OmegaE is the Earth's angular velocity in radians per second.

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

      TYPE :: value_and_derivs_r_and_k                             ! 29 Jul 2016
         REAL(KIND=8) :: value
         REAL(KIND=8) :: pt,pr,pth,pph
         REAL(KIND=8) :: pom,pkr,pkth,pkph
         REAL(KIND=8) :: f0,f1,f2                                  !  7 Jul 2017
         REAL(KIND=8) :: fr,fth,fph                                !  4 Aug 2016
         REAL(KIND=8) :: frr,frth,frph,fthth,fthph,fphph           !  4 Aug 2016
      END TYPE 

      TYPE (value_and_derivs_r_and_k) :: term2                     ! 29 Jul 2016

      TYPE :: values_and_derivs_1_2
         REAL(KIND=8) :: value
         REAL(KIND=8) :: ptm,pr,pth,pph
         REAL(KIND=8) :: ptmtm,prtm,pthtm,pphtm
         REAL(KIND=8) :: ptmr,prr,pthr,pphr
         REAL(KIND=8) :: ptmth,prth,pthth,pphth
         REAL(KIND=8) :: ptmph,prph,pthph,pphph
      END TYPE 

      TYPE (values_and_derivs_1_2) :: rho
      TYPE (values_and_derivs_1_2) :: potentialrho

      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    
      TYPE (gravity_and_derivs_1) :: kA   

      double precision kdotOmega

      double precision kr,kth,kph,ks    
      COMMON//r,th,ph,kr,kth,kph,rkvars(14),tpulse,cstep,drdt(20)

      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)
     & ,(W(606),factor606)                                         

      sinth = dsin(th)
      costh = dcos(th)
      ks = kr**2 + kth**2 + kph**2

      CALL GRAVg(gravity,coriolis)  

      kdotOmega = kr*coriolis%gsr 
     &  + kth*coriolis%gsth 
     &  + kph*coriolis%gsph 

      term2%f0=0.d0
      term2%f1=0.d0                                                !  7 Jul 2017
      term2%f2=0.d0
      if(ks.ne.0.d0) term2%f2 = factor606*4.d0*kdotOmega**2/ks     !  4 Aug 2016
!      term2%f2=(term2%frr*kr*kr+term2%fthth*kth*kth+term2%fphph*kph*kph
!     & +2.d0*(term2%frth*kr*kth+term2%frph*kr*kph+term2%fthph*kth*kph))
!     & /ks
      term2%fr=0.d0
      term2%fth=0.d0
      term2%fph=0.d0
      term2%frr = factor606*4.d0*coriolis%gsr*coriolis%gsr
      term2%frth = factor606*8.d0*coriolis%gsr*coriolis%gsth       ! 31 Jul 2017
      term2%frph = factor606*8.d0*coriolis%gsr*coriolis%gsth       ! 31 Jul 2017
      term2%fthth = factor606*4.d0*coriolis%gsth*coriolis%gsth
      term2%fthph = factor606*8.d0*coriolis%gsth*coriolis%gsph     ! 31 Jul 2017
      term2%fphph = factor606*4.d0*coriolis%gsph*coriolis%gsph

      term2%value = factor606*4.d0*kdotOmega**2      
      term2%pt = factor606*8.d0*kdotOmega*(     
     &    kr*coriolis%prtm 
     & +  kth*coriolis%pthtm 
     & +  kph*coriolis%pphtm
     & )
      term2%pr = factor606*8.d0*kdotOmega*(       
     &    kr*coriolis%prr 
     & +  kth*coriolis%pthr 
     & +  kph*coriolis%pphr
     & )
      term2%pth = factor606*8.d0*kdotOmega*(   
     &    kr*coriolis%prth 
     & +  kth*coriolis%pthth 
     & +  kph*coriolis%pphth
     & )
      term2%pph = factor606*8.d0*kdotOmega*(  
     &    kr*coriolis%prph 
     & +  kth*coriolis%pthph 
     & +  kph*coriolis%pphph
     & )

      term2%pom=0.d0                                               ! 29 Jul 2016
      term2%pkr=factor606*8.d0*kdotOmega*coriolis%gsr              ! 29 Jul 2016
      term2%pkth=0.d0
      term2%pkth=factor606*8.d0*kdotOmega*coriolis%gsth            ! 29 Jul 2016
      term2%pkph=0.d0
      term2%pkph=factor606*8.d0*kdotOmega*coriolis%gsph            ! 29 Jul 2016

c      print*,'Coriolis term 2 = ',term2%value

      RETURN
C ----------------------------------------------------------------------!-------
      END SUBROUTINE coriolis2
C ----------------------------------------------------------------------!-------


C ----------------------------------------------------------------------!-------
      SUBROUTINE coriolis3(term3)                                  ! 22 Jul 2016
C ----------------------------------------------------------------------!-------
C Term3 equals 4 times the square of (k sub z times OmegaE sub z), where
!     OmegaE is the Earth's angular velocity in radians per second.

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

      TYPE :: value_and_derivs_r_and_k                             ! 29 Jul 2016
         REAL(KIND=8) :: value
         REAL(KIND=8) :: pt,pr,pth,pph
         REAL(KIND=8) :: pom,pkr,pkth,pkph
         REAL(KIND=8) :: f0,f1,f2                                  !  7 Jul 2017
         REAL(KIND=8) :: fr,fth,fph                                !  4 Aug 2016
         REAL(KIND=8) :: frr,frth,frph,fthth,fthph,fphph           !  4 Aug 2016
      END TYPE 

      TYPE (value_and_derivs_r_and_k) :: term3                     ! 29 Jul 2016

      TYPE :: values_and_derivs_1_2
         REAL(KIND=8) :: value
         REAL(KIND=8) :: ptm,pr,pth,pph
         REAL(KIND=8) :: ptmtm,prtm,pthtm,pphtm
         REAL(KIND=8) :: ptmr,prr,pthr,pphr
         REAL(KIND=8) :: ptmth,prth,pthth,pphth
         REAL(KIND=8) :: ptmph,prph,pthph,pphph
      END TYPE 

      TYPE (values_and_derivs_1_2) :: rho
      TYPE (values_and_derivs_1_2) :: potentialrho

      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    
      TYPE (gravity_and_derivs_1) :: kA   

      double precision kdotOmega,ks

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

      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)
     & ,(W(607),factor607)          

      sinth = dsin(th)
      costh = dcos(th)
      ks = kr**2 + kth**2 + kph**2

      CALL GRAVg(gravity,coriolis)  

      kdotOmega = kr*coriolis%gsr 

      term3%f0=0.d0
      term3%f1=0.d0                                                !  7 Jul 2017
      term3%f2=0.d0
      if(ks.ne.0.d0) term3%f2 = factor607*4.d0*kdotOmega**2/ks     !  4 Aug 2016
      term3%fr=0.d0
      term3%fth=0.d0
      term3%fph=0.d0
      term3%frr = factor607*4.d0*coriolis%gsr**2                   !  4 Aug 2016
      term3%frth=0.d0
      term3%frph=0.d0
      term3%fthth=0.d0
      term3%fthph=0.d0
      term3%fphph=0.d0

      term3%value = factor607*4.d0*kdotOmega**2   
      term3%pt = factor607*8.d0*kdotOmega*(     
     &    kr*coriolis%prtm )
      term3%pr = factor607*8.d0*kdotOmega*(  
     &    kr*coriolis%prr )
      term3%pth = factor607*8.d0*kdotOmega*(  
     &    kr*coriolis%prth )
      term3%pph = factor607*8.d0*kdotOmega*(   
     &    kr*coriolis%prph )

      term3%pom=0.d0                                               ! 29 Jul 2016
      term3%pkr=2.d0*factor607*4.d0*kdotOmega*coriolis%gsr         ! 29 Jul 2016
      term3%pkth=0.d0                                              ! 29 Jul 2016
      term3%pkph=0.d0                                              ! 29 Jul 2016

c      print*,'Coriolis term 3 = ',term3%value

      RETURN
C ----------------------------------------------------------------------!-------
      END SUBROUTINE coriolis3
C ----------------------------------------------------------------------!-------


C ----------------------------------------------------------------------!-------
      SUBROUTINE coriolis4(term4)                                  ! 22 Jul 2016
C ----------------------------------------------------------------------!-------
! Term4 equals -4 times the square of (omega times OmegaE/C), where
!     OmegaE is the Earth's angular velocity in radians per second.

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

      TYPE :: value_and_derivs_r_and_k                             ! 29 Jul 2016
         REAL(KIND=8) :: value
         REAL(KIND=8) :: pt,pr,pth,pph
         REAL(KIND=8) :: pom,pkr,pkth,pkph
         REAL(KIND=8) :: f0,f1,f2                                  !  7 Jul 2017
         REAL(KIND=8) :: fr,fth,fph                                !  4 Aug 2016
         REAL(KIND=8) :: frr,frth,frph,fthth,fthph,fphph           !  4 Aug 2016
      END TYPE 

      TYPE (value_and_derivs_r_and_k) :: term4                     ! 29 Jul 2016

      TYPE :: values_and_derivs_1_2
         REAL(KIND=8) :: value
         REAL(KIND=8) :: ptm,pr,pth,pph
         REAL(KIND=8) :: ptmtm,prtm,pthtm,pphtm
         REAL(KIND=8) :: ptmr,prr,pthr,pphr
         REAL(KIND=8) :: ptmth,prth,pthth,pphth
         REAL(KIND=8) :: ptmph,prph,pthph,pphph
      END TYPE 

      TYPE (values_and_derivs_1_2) :: rho
      TYPE (values_and_derivs_1_2) :: potentialrho

      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    
      TYPE (gravity_and_derivs_1) :: kA   

!      double precision kAdotOmega

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

C     COMMON DECK "CC" INSERTED HERE
      DOUBLE PRECISION MODC
      COMMON/CC/MODC(4),CS,PCST,PCSR,PCSTH,PCSPH

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

      EQUIVALENCE (W(1),earthr),(w(6),ow)

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

      sinth = dsin(th)
      costh = dcos(th)

      ows = ow**2

      CALL GRAVg(gravity,coriolis)  

      Omegasquare = coriolis%gsr*coriolis%gsr 
     &  + coriolis%gsth*coriolis%gsth + coriolis%gsph*coriolis%gsph 

      term4%f0 = -factor608*4.d0*ows*Omegasquare/cs                !  4 Aug 2016
      term4%f1=0.d0                                                !  7 Jul 2017
      term4%f2=0.d0
      term4%fr=0.d0
      term4%fth=0.d0
      term4%fph=0.d0
      term4%frr=0.d0
      term4%frth=0.d0
      term4%frph=0.d0
      term4%fthth=0.d0
      term4%fthph=0.d0
      term4%fphph=0.d0

      term4%value = -factor608*4.d0*ows*Omegasquare/cs 
      term4%pt = -factor608*8.d0*ows*(    
     &    coriolis%gsr*coriolis%prtm 
     & +  coriolis%gsth*coriolis%pthtm 
     & +  coriolis%gsph*coriolis%pphtm
     & )/cs 
     & + factor608*4.d0*ows*Omegasquare/cs**2*pcst
      term4%pr = -factor608*8.d0*ows*(  
     &    coriolis%gsr*coriolis%prr 
     & +  coriolis%gsth*coriolis%pthr 
     & +  coriolis%gsph*coriolis%pphr
     & )/cs
     & + factor608*4.d0*ows*Omegasquare/cs**2*pcsr
      term4%pth = -factor608*8.d0*ows*(  
     &    coriolis%gsr*coriolis%prth 
     & +  coriolis%gsth*coriolis%pthth 
     & +  coriolis%gsph*coriolis%pphth
     & )/cs
     & + factor608*4.d0*ows*Omegasquare/cs**2*pcsth
      term4%pph = -factor608*8.d0*ows*(   
     &    coriolis%gsr*coriolis%prph 
     & +  coriolis%gsth*coriolis%pthph 
     & +  coriolis%gsph*coriolis%pphph
     & )/cs
     & + factor608*4.d0*ows*Omegasquare/cs**2*pcsph

      term4%pom=2.d0*term4%value/ow                                ! 29 Jul 2016
      term4%pkr=0.d0                                               ! 29 Jul 2016
      term4%pkth=0.d0                                              ! 29 Jul 2016
      term4%pkph=0.d0                                              ! 29 Jul 2016

c      print*,'Coriolis term 4 = ',term4%value

      RETURN
C ----------------------------------------------------------------------!-------
      END SUBROUTINE coriolis4
C ----------------------------------------------------------------------!-------


C ----------------------------------------------------------------------!-------
      SUBROUTINE coriolis5(term5)                                  ! 22 Jul 2016
C ----------------------------------------------------------------------!-------
! Term5 equals -4 times the square of (omega times OmegaE sub z/C), where
!     OmegaE is the Earth's angular velocity in radians per second.

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

      TYPE :: value_and_derivs_r_and_k                             ! 29 Jul 2016
         REAL(KIND=8) :: value
         REAL(KIND=8) :: pt,pr,pth,pph
         REAL(KIND=8) :: pom,pkr,pkth,pkph
         REAL(KIND=8) :: f0,f1,f2                                  !  7 Jul 2017
         REAL(KIND=8) :: fr,fth,fph                                !  4 Aug 2016
         REAL(KIND=8) :: frr,frth,frph,fthth,fthph,fphph           !  4 Aug 2016
      END TYPE 

      TYPE (value_and_derivs_r_and_k) :: term5                     ! 29 Jul 2016

      TYPE :: values_and_derivs_1_2
         REAL(KIND=8) :: value
         REAL(KIND=8) :: ptm,pr,pth,pph
         REAL(KIND=8) :: ptmtm,prtm,pthtm,pphtm
         REAL(KIND=8) :: ptmr,prr,pthr,pphr
         REAL(KIND=8) :: ptmth,prth,pthth,pphth
         REAL(KIND=8) :: ptmph,prph,pthph,pphph
      END TYPE 

      TYPE (values_and_derivs_1_2) :: rho
      TYPE (values_and_derivs_1_2) :: potentialrho

      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    
      TYPE (gravity_and_derivs_1) :: kA   

!      double precision kAdotOmega

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

C     COMMON DECK "CC" INSERTED HERE
      DOUBLE PRECISION MODC
      COMMON/CC/MODC(4),CS,PCST,PCSR,PCSTH,PCSPH

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

      EQUIVALENCE (W(1),earthr),(w(6),ow)

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

      sinth = dsin(th)
      costh = dcos(th)

      ows = ow**2

      CALL GRAVg(gravity,coriolis)  

      Omegasquare = coriolis%gsr*coriolis%gsr 
!     &  + coriolis%gsth*coriolis%gsth + coriolis%gsph*coriolis%gsph 

      term5%f0 = -factor609*4.d0*ows*Omegasquare/cs                !  4 Aug 2016
      term5%f1=0.d0                                                !  7 Jul 2017
      term5%f2=0.d0
      term5%fr=0.d0
      term5%fth=0.d0
      term5%fph=0.d0
      term5%frr=0.d0
      term5%frth=0.d0
      term5%frph=0.d0
      term5%fthth=0.d0
      term5%fthph=0.d0
      term5%fphph=0.d0

      term5%value = -factor609*4.d0*ows*Omegasquare/cs 
      term5%pt = -factor609*8.d0*ows*(   
     &    coriolis%gsr*coriolis%prtm 
!     & +  coriolis%gsth*coriolis%pthtm 
!     & +  coriolis%gsph*coriolis%pphtm
     & )/cs
     & + factor609*4.d0*ows*Omegasquare/cs**2*pcst
      term5%pr = -factor609*8.d0*ows*( 
     &    coriolis%gsr*coriolis%prr 
!     & +  coriolis%gsth*coriolis%pthr 
!     & +  coriolis%gsph*coriolis%pphr
     & )/cs
     & + factor609*4.d0*ows*Omegasquare/cs**2*pcsr
      term5%pth = -factor609*8.d0*ows*(  
     &    coriolis%gsr*coriolis%prth 
!     & +  coriolis%gsth*coriolis%pthth 
!     & +  coriolis%gsph*coriolis%pphth
     & )/cs
     & + factor609*4.d0*ows*Omegasquare/cs**2*pcsth
      term5%pph = -factor609*8.d0*ows*(  
     &    coriolis%gsr*coriolis%prph 
!     & +  coriolis%gsth*coriolis%pthph 
!     & +  coriolis%gsph*coriolis%pphph
     & )/cs
     & + factor609*4.d0*ows*Omegasquare/cs**2*pcsph
      term5%pom=2.d0*term5%value/ow                                ! 29 Jul 2016
      term5%pkr=0.d0                                               ! 29 Jul 2016
      term5%pkth=0.d0                                              ! 29 Jul 2016
      term5%pkph=0.d0                                              ! 29 Jul 2016

c      print*,'Coriolis term 5 = ',term5%value

      RETURN
C ----------------------------------------------------------------------!-------
      END SUBROUTINE coriolis5
C ----------------------------------------------------------------------!-------


C ----------------------------------------------------------------------!-------
      SUBROUTINE coriolis6(term6)                                  ! 13 Sep 2016
C ----------------------------------------------------------------------!-------
C Term6 equals 4 times omega OmegaE cross Gamma dot k, where
!     OmegaE is the Earth's angular velocity in radians per second.

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

      TYPE :: value_and_derivs_r_and_k                            
         REAL(KIND=8) :: value
         REAL(KIND=8) :: pt,pr,pth,pph
         REAL(KIND=8) :: pom,pkr,pkth,pkph
         REAL(KIND=8) :: f0,f1,f2                                  !  7 Jul 2017
         REAL(KIND=8) :: fr,fth,fph                               
         REAL(KIND=8) :: frr,frth,frph,fthth,fthph,fphph          
      END TYPE 

      TYPE (value_and_derivs_r_and_k) :: term6                    

      TYPE :: values_and_derivs_1_2
         REAL(KIND=8) :: value
         REAL(KIND=8) :: ptm,pr,pth,pph
         REAL(KIND=8) :: ptmtm,prtm,pthtm,pphtm
         REAL(KIND=8) :: ptmr,prr,pthr,pphr
         REAL(KIND=8) :: ptmth,prth,pthth,pphth
         REAL(KIND=8) :: ptmph,prph,pthph,pphph
      END TYPE 

      TYPE (values_and_derivs_1_2) :: rho
      TYPE (values_and_derivs_1_2) :: potentialrho
      TYPE (values_and_derivs_1_2) :: pressure                     ! 12 Sep 2016

      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    
      TYPE (gravity_and_derivs_1) :: kA   
      TYPE (gravity_and_derivs_1) :: Gamma                         ! 12 Sep 2016

      double precision kAdotOmega

      double precision kr,kth,kph,ks                               !  7 Jul 2017
      COMMON//r,th,ph,kr,kth,kph,rkvars(14),tpulse,cstep,drdt(20)

C     COMMON DECK "CC" INSERTED HERE
      DOUBLE PRECISION MODC
      COMMON/CC/MODC(4),CS,PCST,PCSR,PCSTH,PCSPH

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

      EQUIVALENCE (W(1),earthr),(W(6),ow)

C     GRAVITY  600-624
      EQUIVALENCE (W(600),gvmodel),(W(601),gvform),(W(602),gvid)
     & ,(W(610),factor610)                                         !  1 Jul 2016

c      print *, ' before density'
      CALL DENSITYg(rho,potentialrho) 
c      print *, ' after density'

      CALL PRESg(pressure)                                        ! 12 Sep 2016

      sinth = dsin(th)
      costh = dcos(th)

      kA%gsr = (rho%pr/(2.d0*rho%value))  

      kA%prtm = (rho%prtm)/(2.d0*rho%value) -
     & ((rho%pr)*rho%ptm)/(2.d0*rho%value**2)
      kA%prr = (rho%prr)/(2.d0*rho%value) -
     & ((rho%pr**2))/(2.d0*rho%value**2)
      kA%prth = (rho%prth)/(2.d0*rho%value) -
     & ((rho%pr)*rho%pth)/(2.d0*rho%value**2)
      kA%prph = (rho%prph)/(2.d0*rho%value) - 
     & ((rho%pr)*rho%pph)/(2.d0*rho%value**2)

      kA%gsth = (rho%pth/(2.d0*r*rho%value)) 

      kA%pthtm = (rho%pthtm)/(2.d0*r*rho%value) -
     & ((rho%pth)*rho%ptm)/(2.d0*r*rho%value**2)
      kA%pthr = (rho%pthr)/(2.d0*r*rho%value) -
     & ((rho%pth*rho%pr))/(2.d0*r*rho%value**2) -
     & ((rho%pth))/(2.d0*r**2*rho%value)
      kA%pthth = (rho%pthth)/(2.d0*r*rho%value) -
     & ((rho%pth)*rho%pth)/(2.d0*r*rho%value**2)
      kA%pthph = (rho%pthph)/(2.d0*r*rho%value) - 
     & ((rho%pth)*rho%pph)/(2.d0*r*rho%value**2)

      kA%gsph = (rho%pph/(2.d0*r*sinth*rho%value))

      kA%pphtm = (rho%pphtm)/(2.d0*r*sinth*rho%value) -
     & ((rho%pph)*rho%ptm)/(2.d0*r*sinth*rho%value**2)
      kA%pphr = (rho%pphr)/(2.d0*r*sinth*rho%value) -
     & ((rho%pph*rho%pr))/(2.d0*r*sinth*rho%value**2) -
     & ((rho%pph))/(2.d0*r**2*sinth*rho%value)
      kA%pphth = (rho%pphth)/(2.d0*r*sinth*rho%value) -
     & ((rho%pph)*rho%pth)/(2.d0*r*sinth*rho%value**2) -
     & ((rho%pph)*costh)/(2.d0*r*sinth**2*rho%value)
      kA%pphph = (rho%pphph)/(2.d0*r*sinth*rho%value) - 
     & ((rho%pph)*rho%pph)/(2.d0*r*sinth*rho%value**2)

! The next 5 statements use only the vertical component of kA, 
! but that is OK, because they are not actually needed for Coriolis term 6

      kA%gs = kA%gsr 
      kA%ptm = kA%prtm 
      kA%pr = kA%prr 
      kA%pth = kA%prth 
      kA%pph = kA%prph 

      CALL GRAVg(gravity,coriolis)  

Calculating Gamma

      Gamma%gsr = kA%gsr - (pressure%pr/(2.d0*rho%value*cs))  

      Gamma%prtm = kA%prtm - (pressure%prtm)/(2.d0*rho%value*cs) 
     & +((pressure%pr)*rho%ptm)/(2.d0*rho%value**2*cs)
     & +((pressure%pr)*pcst)/(2.d0*rho%value*cs**2)
      Gamma%prr = kA%prr - (pressure%prr)/(2.d0*rho%value*cs) 
     & +((pressure%pr)*rho%pr)/(2.d0*rho%value**2*cs)
     & +((pressure%pr)*pcsr)/(2.d0*rho%value*cs**2)
      Gamma%prth = kA%prth - (pressure%prth)/(2.d0*rho%value*cs) 
     & +((pressure%pr)*rho%pth)/(2.d0*rho%value**2*cs)
     & +((pressure%pr)*pcsth)/(2.d0*rho%value*cs**2)
      Gamma%prph = kA%prph - (pressure%prph)/(2.d0*rho%value*cs) 
     & +((pressure%pr)*rho%pph)/(2.d0*rho%value**2*cs)
     & +((pressure%pr)*pcsph)/(2.d0*rho%value*cs**2)

      Gamma%gsth = kA%gsth - (pressure%pth/(2.d0*r*rho%value*cs)) 

      Gamma%pthtm = kA%pthtm - (pressure%pthtm)/(2.d0*r*rho%value*cs) 
     & +((pressure%pth)*rho%ptm)/(2.d0*r*rho%value**2*cs)
     & +((pressure%pth)*pcst)/(2.d0*r*rho%value*cs**2)
      Gamma%pthr = kA%pthr - (pressure%pthr)/(2.d0*r*rho%value*cs) 
     & +((pressure%pth*rho%pr))/(2.d0*r*rho%value**2*cs) 
     & +((pressure%pth))/(2.d0*r**2*rho%value*cs)
     & +((pressure%pth*rho%pr))/(2.d0*r*rho%value*cs**2) 
      Gamma%pthth = kA%pthth - (pressure%pthth)/(2.d0*r*rho%value*cs) 
     & +((pressure%pth)*rho%pth)/(2.d0*r*rho%value**2*cs)
     & +((pressure%pth)*pcsth)/(2.d0*r*rho%value*cs**2)
      Gamma%pthph = kA%pthph - (pressure%pthph)/(2.d0*r*rho%value*cs) 
     & +((pressure%pth)*rho%pph)/(2.d0*r*rho%value**2*cs)
     & +((pressure%pth)*pcsph)/(2.d0*r*rho%value*cs**2)

      Gamma%gsph = kA%gsph - (pressure%pph/(2.d0*r*sinth*rho%value*cs))

      Gamma%pphtm = kA%pphtm 
     & - (pressure%pphtm)/(2.d0*r*sinth*rho%value*cs) 
     & +((pressure%pph)*rho%ptm)/(2.d0*r*sinth*rho%value**2*cs)
     & +((pressure%pph)*pcst)/(2.d0*r*sinth*rho%value*cs**2)
      Gamma%pphr = kA%pphr 
     & - (pressure%pphr)/(2.d0*r*sinth*rho%value*cs) 
     & +((pressure%pph*rho%pr))/(2.d0*r*sinth*rho%value**2*cs) 
     & +((pressure%pph))/(2.d0*r**2*sinth*rho%value*cs)
     & +((pressure%pph*pcsr))/(2.d0*r*sinth*rho%value*cs**2) 
      Gamma%pphth = kA%pphth 
     & - (pressure%pphth)/(2.d0*r*sinth*rho%value*cs) 
     & +((pressure%pph)*rho%pth)/(2.d0*r*sinth*rho%value**2*cs) 
     & +((pressure%pph)*costh)/(2.d0*r*sinth**2*rho%value*cs)
     & +((pressure%pph)*pcsth)/(2.d0*r*sinth*rho%value*cs**2) 
      Gamma%pphph = kA%pphph 
     & - (pressure%pphph)/(2.d0*r*sinth*rho%value*cs)  
     & +((pressure%pph)*rho%pph)/(2.d0*r*sinth*rho%value**2*cs)
     & +((pressure%pph)*pcsph)/(2.d0*r*sinth*rho%value*cs**2)

! The next 5 statements use only the vertical component of Gamma, 
! but that is OK, because they are not actually needed for Coriolis term 6

      Gamma%gs = Gamma%gsr
      Gamma%ptm = Gamma%prtm 
      Gamma%pr = Gamma%prr 
      Gamma%pth = Gamma%prth 
      Gamma%pph = Gamma%prph 



      kAdotOmega = kA%gsr*coriolis%gsr + kA%gsth*coriolis%gsth
     &  + kA%gsph*coriolis%gsph 

      term6%f0=0.d0
      term6%f1=0.d0                                                !  7 Jul 2017
      term6%f2=0.d0
      term6%fr=factor610*4.d0*ow*(
     &  coriolis%gsth*Gamma%gsph
     & -coriolis%gsph*Gamma%gsth
     & )
      term6%fth=factor610*4.d0*ow*(
     &  coriolis%gsph*Gamma%gsr
     & -coriolis%gsr*Gamma%gsph
     & )
      term6%fph=factor610*4.d0*ow*(
     &  coriolis%gsr*Gamma%gsth
     & -coriolis%gsth*Gamma%gsr
     & )
      term6%frr=0.d0
      term6%frth=0.d0
      term6%frph=0.d0
      term6%fthth=0.d0
      term6%fthph=0.d0
      term6%fphph=0.d0

      term6%value = term6%fr*kr+term6%fth*kth+term6%fph*kph
      ks = kr**2 + kth**2 + kph**2                                 !  7 Jul 2017
      if(ks.ne.0.d0) term6%f1=term6%value/sqrt(ks)                 !  7 Jul 2017

      term6%pt = factor610*4.d0*ow*(
     & kr*(coriolis%gsth*Gamma%pphtm+coriolis%pthtm*Gamma%gsph
     & -coriolis%gsph*Gamma%pthtm-coriolis%pphtm*Gamma%gsth) +
     & kth*(coriolis%gsph*Gamma%prtm+coriolis%pphtm*Gamma%gsr
     & -coriolis%gsr*Gamma%pphtm-coriolis%prtm*Gamma%gsph) +
     & kph*(coriolis%gsr*Gamma%pthtm+coriolis%prtm*Gamma%gsth
     & -coriolis%gsth*Gamma%prtm-coriolis%pthtm*Gamma%gsr))
      term6%pr = factor610*4.d0*ow*(
     & kr*(coriolis%gsth*Gamma%pphr+coriolis%pthr*Gamma%gsph
     & -coriolis%gsph*Gamma%pthr-coriolis%pphr*Gamma%gsth) +
     & kth*(coriolis%gsph*Gamma%prr+coriolis%pphr*Gamma%gsr
     & -coriolis%gsr*Gamma%pphr-coriolis%prr*Gamma%gsph) +
     & kph*(coriolis%gsr*Gamma%pthr+coriolis%prr*Gamma%gsth
     & -coriolis%gsth*Gamma%prr-coriolis%pthr*Gamma%gsr))
      term6%pth = factor610*4.d0*ow*(
     & kr*(coriolis%gsth*Gamma%pphth+coriolis%pthth*Gamma%gsph
     & -coriolis%gsph*Gamma%pthth-coriolis%pphth*Gamma%gsth) +
     & kth*(coriolis%gsph*Gamma%prth+coriolis%pphth*Gamma%gsr
     & -coriolis%gsr*Gamma%pphth-coriolis%prth*Gamma%gsph) +
     & kph*(coriolis%gsr*Gamma%pthth+coriolis%prth*Gamma%gsth
     & -coriolis%gsth*Gamma%prth-coriolis%pthth*Gamma%gsr))
      term6%pph = factor610*4.d0*ow*(
     & kr*(coriolis%gsth*Gamma%pphph+coriolis%pthph*Gamma%gsph
     & -coriolis%gsph*Gamma%pthph-coriolis%pphph*Gamma%gsth) +
     & kth*(coriolis%gsph*Gamma%prph+coriolis%pphph*Gamma%gsr
     & -coriolis%gsr*Gamma%pphph-coriolis%prph*Gamma%gsph) +
     & kph*(coriolis%gsr*Gamma%pthph+coriolis%prph*Gamma%gsth
     & -coriolis%gsth*Gamma%prph-coriolis%pthph*Gamma%gsr))

      term6%pom=term6%value/ow
      term6%pkr=term6%fr
      term6%pkth=term6%fth
      term6%pkph=term6%fph

c      print*,'ka = ',kA%gsr
c      print*,'Coriolis term 6 = ',term6%value

      RETURN
C ----------------------------------------------------------------------!-------
      END SUBROUTINE coriolis6
C ----------------------------------------------------------------------!-------

