C ----------------------------------------------------------------------!-------
C     gwwnl.for
C ----------------------------------------------------------------------!-------
!  6 Jul 2018 - Correct calculation of kz2a and kz2ai
! 28 Jun 2018 - Added Coriolis terms to calculation of kz2a
! 27 Jun 2018 - Added common block /rink4/ to communicate integer propagate(i)
! 31 Jul 2017 - Added in czz0a, czz0b, and other terms.
! 31 Jul 2017 - Declared kperpvperp double precision
! 31 Jul 2017 - Corrected ddz(0) and ddz(1)
! 28 Jul 2017 - Corrected Coriolis2
! 28 Jul 2017 - Added Coriolis terms to the formulas for ground reflection 
!               (valid only for reflection from a horizontal surface)
! 27 Jul 2017 - Corrected error in Coriolis5
! 21 Jul 2017 - Corrected Coriolis additions to dd(0), dd(1), and dd(2)
! 17 Jul 2017 - CORRECTIONS TO HAMILTONIAN
! 14 Jul 2017 - Added TYPE :: gravity_and_derivs_1
! 11 Jul 2017 - Added Coriolis terms to the Hamiltonian and its derivatives
! 11 Jul 2017 - Made additions to the calculation of dd(0), dd(1), and dd(2)
!               for Coriolis terms
! 11 Jul 2017 - Added calculation of derivatives with respect to 
!               ow, kr, kth, and kph to Coriolis6
! 10 Jul 2017 - Added calculation of derivatives to Coriolis4 through Coriolis6
! 10 Jul 2017 - Added calculation of frr to fthph to Coriolis4 through Coriolis6
! 10 Jul 2017 - Added calculation of fr, fth, fph to Coriolis4 and Coriolis5
!  7 Jul 2017 - Added calculation of fr, fth, fph to Coriolis6
!  7 Jul 2017 - Added calculation of f0, f1, f2 to Coriolis4 through Coriolis6
!  7 Jul 2017 - Added f1 to TYPE :: value_and_derivs_r_and_k
!  6 Jul 2017 - Added wind velocity to Coriolis4 through Coriolis6
!  6 Jul 2017 - Changed some REAL (KIND=8) TO DOUBLE PRECISION
!  6 Jul 2017 - Added subroutines Coriolis1 through Coriolis6
!  6 Jul 2017 - Added some double precision declarations
!  5 Aug 2016 - Added common block /rink3/ to help with reflection of wave
! 27 Jun 2016 - Added terms to the /GV/ common block for Coriolis terms
!  3 Jun 2013 - The identification of roots of the quartic is still not perfect
! 31 May 2013 - Corrected the identification of roots
! 29 May 2013 - Switching from using identifyroot to using identifyroot2
! 28 May 2013 - Added subroutine identifyroot2
! 24 May 2013 - Added kv and propagate to the argument list for subroutine quart
! 24 May 2013 - Added Nsol, kv, nsmall2, and propagate to the 
!               argument list for subroutine identifyroot
!  9 Mar 2012 - Adjust for different definition of ray (one smaller)
!  9 Mar 2012 - Resets RAYNAME to the 3 types of gravity waves (O, X, lg X)
!  2 Mar 2012 - Added the possibility of choosing any root of the quartic
! 13 Apr 2011 - Corrected calculation of kz2a for case where kr2, kth2, or kph2 
!               might be negative
! 13 Apr 2011 - Added common block /rink2/
! 11 Mar 2011 - allowed calculation of error when kay2=0
! 29 Nov 2010 - Altered the algorithm for identifying roots of the quartic
! 22 Nov 2010 - Corrected the use of ray=W2
! 19 Nov 2010 - Added calculation of errRell,errRelb,errRela 
! 19 Nov 2010 - Added errRell,errRelb,errRela to /rink/
! 10 Aug 2010 - Altered the algorithm for identifying roots of the quartic
!  9 Aug 2010 - Altered the algorithm for identifying roots of the quartic
!  5 Aug 2010 - Sets RAYNAME to the 3 types of gravity waves
!  5 Aug 2010 - Added kaysqi,kaysqbi,kaysqai to /rink/
!  3 Aug 2010 - Altered the algorithm for identifying roots of the quartic
!  1 Aug 2010 - Altered the algorithm for identifying roots of the quartic
! 30 Jul 2010 - Altered the algorithm for identifying roots of the quartic
! 28 Jul 2010 - Added alpha,dd to /rink/
! 28 Jul 2010 - Added subroutine identifyroot
! 27 Jul 2010 - Added calculation of Hrel
! 27 Jul 2010 - Added calculation of errRel
! 27 Jul 2010 - Added Hrel,errRel to /rink/
! 27 Jul 2010 - Cross-wind calculations
! 26 Jul 2010 - Corrections
! 23 Jul 2010 - Used a different algorithm to identify roots for large wind
! 23 Jul 2010 - Added test on ray type to no-wind case
! 20 Jul 2010 - Added iteration subroutine
! 19 Jul 2010 - Test for convergence
! 16 Jul 2010 - Added some tests to avoid square roots of negative numbers
! 15 Jul 2010 - Added printing of n3, n4, nbig, nsmall
! 15 Jul 2010 - Added n3,n4,nbig,nsmall to common block /rink/
! 15 Jul 2010 - Improve identification of roots 
! 28 Jun 2010 - Choose ray type from w(2) in data
! 28 Jun 2010 - Added equivalence of ray with W(2) for type of ray
! 28 Jun 2010 - Added relerr,relerra,relerrb to common block /rink/
! 28 Jun 2010 - Added subroutines bigroot, smallroot, root3, and root4
! 24 Jun 2010 - better organization of testing
! 23 Jun 2010 - better organization of testing
! 23 Jun 2010 - More accurate calculation of smaller root
! 22 Jun 2010 - More accurate calculation of smaller root
! 18 Jun 2010 - refine the quartic rootfinding
! 16 Jun 2010 - refine the quartic rootfinding
! 15 Jun 2010 - refine the quartic rootfinding
! 14 Jun 2010 - refine the quartic rootfinding
! 11 Jun 2010 - refine the quartic rootfinding
! 28 May 2010 - Get error from drdt(19)
! 27 May 2010 - Added test printing depending on drdt(20), set in reflect in s2
! 30 Apr 2010 - Slight change in setting model labeling
! 29 Apr 2010 - Finished debugging
! 27 Apr 2010 - Used a solution of a quartic equation to calculate kay2
! 14 Apr 2010 - Deleted Subroutine Setwnd
! 14 Apr 2010 - Added wind
! 14 Apr 2010 - Copied from gnwnl.for
C ----------------------------------------------------------------------!-------
! 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: WITH WIND, NO LOSSES
! Barotropic case, assumes the gradients of pressure and density are vertical

      SUBROUTINE GWWNL ! (WITH WINDS NO LOSSES)                    ! 25 Mar 2010

      IMPLICIT double precision (a-h,o-z)
!      complex*16 p4,xxx                                           ! 27 Jul 2010
      double precision errorcalc,kperpvperp          ! 11 Mar 2011 ! 31 Jul 2017
      external errorcalc                                           ! 11 Mar 2011

      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                             !  6 Jul 2017
         REAL(KIND=8) :: value                                     !  6 Jul 2017
         REAL(KIND=8) :: pt,pr,pth,pph                             !  6 Jul 2017
         REAL(KIND=8) :: pom,pkr,pkth,pkph                         !  6 Jul 2017
         REAL(KIND=8) :: f0,f1,f2                                  !  7 Jul 2017
         REAL(KIND=8) :: fr,fth,fph                                !  6 Jul 2017
         REAL(KIND=8) :: frr,frth,frph,fthth,fthph,fphph           !  6 Jul 2017
      END TYPE 

      TYPE (value_and_derivs_r_and_k) :: term1                     !  6 Jul 2017
      TYPE (value_and_derivs_r_and_k) :: term2,term3,term4,term5   !  6 Jul 2017
      TYPE (value_and_derivs_r_and_k) :: term6                     !  6 Jul 2017

      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                                 ! 14 Jul 2017
         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    

      complex*16 test(4),testa(4),testb(4)                         ! 24 Jun 2010
      integer propagate(4)                                         ! 24 May 2013
      COMMON /rink4/ propagate                                     ! 27 Jun 2018

c      REAL(KIND=8) :: kr,kth,kph
      double precision kr,kth,kph                                  !  6 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

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,kay 
     & , 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                               ! 13 Apr 2011
      COMMON /rink2/ kr2,kth2,kph2                                 ! 13 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)                   ! 14 Jun 2010
     & ,kaysqi(4),kaysqbi(4),kaysqai(4)                            !  5 Aug 2010
     & ,sol(4),soli(4),solb(4),solbi(4),sola(4),solai(4)           ! 15 Jun 2010
     & ,alpha,dd                                                   ! 28 Jul 2010
     & ,relerr(4),relerra(4),relerrb(4)                            ! 28 Jun 2010
     & ,errRell(4),errRela(4),errRelb(4)                           ! 19 Nov 2010
     & ,n3,n4,nbig,nsmall                                          ! 15 Jul 2010
     & ,Hrel,errRel                                                ! 27 Jul 2010
      complex*16 relerr,relerra,relerrb                            ! 28 Jun 2010
      complex*16 errRell,errRela,errRelb                           ! 19 Nov 2010
      double precision kaysq,kaysqb,kaysqa                         ! 14 Jun 2010
     & ,kaysqi,kaysqbi,kaysqai                                     !  5 Aug 2010
      dimension dd(0:4)
     & ,ddz(0:4)                                                   ! 11 Jul 2017
!     & ,sol(4),soli(4),solb(4),solbi(4),sola(4),solai(4)          ! 15 Jun 2010

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)
     & ,(RAY,W(2))                                                 ! 28 Jun 2010
      REAL(KIND=8) :: maxstp,maxerr,intyp,llat,llong
c      double precision maxstp,maxerr,intyp,llat,llong              !  6 Jul 2017
c      REAL(KIND=8) :: ks,kv,kperp2
      double precision ks,kv,kperp2                                !  6 Jul 2017
      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,owi,ow2                              !  6 Jul 2017
      EQUIVALENCE(kvect,kay2),(vset,v),(ow,w(6))
     & ,(W(606),factor606)                                         ! 29 Jul 2016
      data errorsave /0.d0/

      ENTRY disper                                                 ! 25 Mar 2010

      space = .false.
      height = r - earthr
      ks = kr**2 + kth**2 + kph**2

C     WIND VELOCITY
      CALL WINDR                                                   ! 14 Apr 2010
      KV=KR*VR+KTH*VTH+KPH*VPH                                     ! 14 Apr 2010
      VLS=KV*KV/KS                                                 ! 14 Apr 2010
! OWI = intrinsic frequency
      OWI=OW-KV                                                    ! 14 Apr 2010
!      ows = ow**2
      ows = owi**2                                                 ! 14 Apr 2010
      ow2 = ow**2

C     SOUND SPEED 
c      print *, ' before speed '
      CALL speed
c      print *, ' after speed '

!  Acoustic cutoff frequency
c      print *, ' before cutoff'
      CALL cutoff(kA2)                                             ! 24 Feb 2006

!  Brunt-Vaisala frequency
c      print *, ' before Brunt '
      CALL Brunt(Ns)                                               ! 24 Feb 2006

      call coriolis1(term1)  !  6 Jul 2017
      call coriolis2(term2)  !  6 Jul 2017
      call coriolis3(term3)  !  6 Jul 2017
      call coriolis4(term4)  !  6 Jul 2017
      call coriolis5(term5)  !  6 Jul 2017
      call coriolis6(term6)  !  6 Jul 2017

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

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

      if(kv.eq.0.d0) then                                          ! 27 Apr 2010

         nbig = 0                                                  ! 26 Jul 2010
         n3 = 1                                                    ! 27 Jul 2010
         n4 = 0                                                    ! 26 Jul 2010
         nsmall = 0                                                ! 26 Jul 2010

         do i=1,4
         kaysq(i)=0.d0                                             ! 26 Jul 2010
         kaysqb(i)=0.d0                                            ! 26 Jul 2010
         kaysqa(i)=0.d0                                            ! 26 Jul 2010
         kaysqi(i)=0.d0                                            !  5 Aug 2010
         kaysqbi(i)=0.d0                                           !  5 Aug 2010
         kaysqai(i)=0.d0                                           !  5 Aug 2010
         sol(i)=0.d0                                               ! 26 Jul 2010
         soli(i)=0.d0                                              ! 26 Jul 2010
         solb(i)=0.d0                                              ! 26 Jul 2010
         solbi(i)=0.d0                                             ! 26 Jul 2010
         sola(i)=0.d0                                              ! 26 Jul 2010
         solai(i)=0.d0                                             ! 26 Jul 2010
         relerr(i) =(0.d0,0.d0)                                    ! 26 Jul 2010
         relerra(i)=(0.d0,0.d0)                                    ! 26 Jul 2010
         relerrb(i)=(0.d0,0.d0)                                    ! 26 Jul 2010
         enddo                                                     ! 26 Jul 2010

         kaysqb(n3)=                                               ! 27 Jul 2010
     &      ((ows/cs)-kA2%value)/((ows-Ns%value)*
     &      (kth**2+kph**2)/(ows*ks) +
     &      (kr**2/ks))
      c=dsqrt(cs)

      if(kaysqb(n3).ge.0.d0) then
         solb(n3)=c*dsqrt(kaysqb(n3))
         solbi(n3)=0.d0
      else
         solb(n3) =0.d0
         solbi(n3)=c*dsqrt(-kaysqb(n3))
      endif

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

         kaysq(n3)=kaysqb(n3)                                      ! 27 Jul 2010
         sol(n3)= solb(n3)                                         ! 27 Jul 2010
         soli(n3)=solbi(n3)                                        ! 27 Jul 2010

         if(ray.eq.0.d0) kay2 = kaysq(n3)            ! 27 Jul 2010 !  9 Mar 2012

         if(ray.ne.0.d0) kay2=-1.d5                  ! 26 Jul 2010 !  9 Mar 2012

            kay2i=0.d0

!      xxx = cmplx(solb(n3),solbi(n3))                             ! 27 Jul 2010
!      testb(n3) = p4(dd,xxx,relerrb(n3))                          ! 27 Jul 2010

!      xxx = cmplx(sol(n3),soli(n3))                               ! 27 Jul 2010
!      test(n3) = p4(dd,xxx,relerr(n3))                            ! 27 Jul 2010

!            print *, 'kay2 nowind disper= ', kay2
!            print *, 'kay2i disper= ', kay2i

      else                                                         ! 27 Apr 2010
!         if(tpulse.eq.0.d0) drdt(20)=23.d0                        ! 27 May 2010
         error = errRel                                            ! 27 Jul 2010
!         errorden = drdt(18)                                      ! 28 May 2010
!         errornum = drdt(17)                                      ! 28 May 2010
!         if(dabs(error).gt.w(42)*5.d3) drdt(20)=23.d0             ! 28 May 2010
!         if(dabs(error).gt.0.05d0) drdt(20)=23.d0                 ! 28 May 2010
!         if(dabs(error).gt.w(42)*4.6507705d3) drdt(20)=24.d0      ! 28 May 2010
!         if(dabs(error).gt.w(42)*1.d4) drdt(20)=24.d0             ! 21 Jul 2010
!         if(dabs(error).gt.w(42)*1.d4) drdt(20)=23.d0             ! 21 Jul 2010
!         if(dabs(error).gt.w(42)*1.d3) drdt(20)=23.d0             ! 27 May 2010
!         if(dabs(error).gt.w(42)*1.d2) drdt(20)=23.d0             !  1 Jun 2010
!         if(dabs(error).gt.w(42)*1.d1) drdt(20)=23.d0             !  1 Jun 2010
!         if(dabs(error).gt.w(42)*5.d1) drdt(20)=23.d0             !  1 Jun 2010
!         if(dabs(error).gt.w(42)*2.d1) drdt(20)=23.d0             !  1 Jun 2010
!         if(dabs(error).gt.w(42)*1.57d1) drdt(20)=23.d0           !  1 Jun 2010
!         drdt(20)=23.d0           
         if(drdt(20).eq.24.d0
     &   .and.error.ne.errorsave
     &   ) 
     &   print *, 'error = ', error 
!     &   , 'kay2 = ', errorden
!     &   , 'k2 = ', errornum

         if(drdt(20).eq.24.d0
     &   .and.error.ne.errorsave
     &   ) 
     &   drdt(20)=23.d0                 ! 28 May 2010

      errorsave=error

! Coefficients of the quartic for Ck
! dd(0)+dd(1)*(Ck)+dd(2)*(Ck)**2+dd(3)*(Ck)**3+dd(4)*(Ck)**4=0     !  5 Jul 2017
         c=dsqrt(cs)
         alpha = vls/cs
         alphalf=kv/(c*dsqrt(ks))
         if(drdt(20).eq.23.d0) print *, 'alpha = ', alpha
     &                                , 'alphalf  = ', alphalf
         kperp2 = kth*kth+kph*kph
         dd(0) = -ow2*(ow2-cs*kA2%value)
     & - (term1%f0+term2%f0+term3%f0+term4%f0+term5%f0+term6%f0)*cs! 21 Jul 2017
         dd(1) = 2.d0*ow*(2.d0*ow2-cs*kA2%value)*alphalf
     & - (term1%f1+term2%f1+term3%f1+term4%f1+term5%f1+term6%f1)*c ! 21 Jul 2017
         dd(2) = ((1.d0-6.d0*alpha)*ow2
     &   +kA2%value*cs*alpha-kperp2*Ns%value/ks)
     & - (term1%f2+term2%f2+term3%f2+term4%f2+term5%f2+term6%f2)   ! 21 Jul 2017
         dd(3) = -2.d0*ow*(1.d0-2.d0*alpha)*alphalf
         dd(4) = (1.d0-alpha)*alpha

         drdt20 = drdt(20)

!      print *, ' ks before quart = ', ks
         call quart(dd,sol,soli,Nsol
     & ,sola,solai,solb,solbi
     & ,test,testa,testb                                           ! 24 Jun 2010
     & ,relerr,relerra,relerrb                                     ! 28 Jun 2010
     & ,n3,n4,nbig,nsmall                                          ! 15 Jul 2010
     & ,height,ks,cs,c,error,errorsave
     & ,kaysq,kaysqa,kaysqb
     & ,kay2,drdt20                                                ! 15 Jun 2010
     & ,kv,propagate)                                              ! 24 May 2013
!      print *, ' ks after quart = ', ks
! Nsol is the number of real solutions of the quartic. kv is k dot v
! propagate(i) is zero if the ith root is evanescent, non-zero if propagating

      drdt20save=drdt20
!      drdt20=23.d0

! Ordinary ray
      IF(ray.eq.0.d0) THEN
         kay2=kaysq(n3)                              ! 28 Jun 2010 !  9 Mar 2012
         IF(propagate(n3).eq.0) kay2 = -abs(kay2)                  ! 29 May 2013
      ENDIF

! Small extraordinary ray
      if(ray.eq.-1.d0) then            ! 28 Jun 2010 ! 22 Nov 2010 !  9 Mar 2012
         kay2=kaysq(nsmall)                                        ! 28 Jun 2010
         if(sol(nsmall).le.0.d0
     &    .or.propagate(nsmall).eq.0                               ! 29 May 2013
     &    ) kay2=-abs(kay2)                                        ! 28 Jun 2010
      endif                                                        ! 28 Jun 2010

! Large extraordinary ray
      if(ray.eq.-2.d0) then            ! 28 Jun 2010 ! 22 Nov 2010 !  9 Mar 2012
         kay2=kaysq(nbig)                                          ! 28 Jun 2010
         if(sol(nbig).le.0.d0
     &    .or.propagate(nbig).eq.0                                 ! 29 May 2013
     &    ) kay2=-abs(kay2)                                        ! 28 Jun 2010
      endif                                                        ! 28 Jun 2010

! Choosing root number
      if(ray.gt.0.d0) then                           !  2 Mar 2012 !  9 Mar 2012
         kay2=kaysq(int(ray))                        !  2 Mar 2012 !  9 Mar 2012
         if(sol(int(ray)).le.0.d0
     &    .or.propagate(int(ray)).eq.0                             ! 29 May 2013
     &    ) kay2=-abs(kay2)   !                         2 Mar 2012 !  9 Mar 2012
      endif                                                        !  2 Mar 2012

!      if(kay2.le.0.d0) print *, 'kay2 = ', kay2

C ----------------------------------------------------------------------!-------
      error=errorcalc(ks,kay2)                                     ! 11 Mar 2011
      if(drdt20.eq.23.d0) print *, 'error = ', error
      if(drdt20.eq.23.d0) print *, 'kay2 = ', kay2

! 28 Jun 2010            if(sol(nsmall).le.0.d0) kay2=kaysq(n3)    ! 16 Jun 2010

      error=errorcalc(ks,kay2)                                     ! 11 Mar 2011
      if(drdt20.eq.23.d0) print *, 'error = ', error
      if(drdt20.eq.23.d0) print *, 'kay2 = ', kay2
      if(drdt20.eq.23.d0) print *, 'leaving subroutine quart'
!      drdt20=drdt20save

         drdt(20)=0.d0

!         stop

      endif                                                        ! 27 Apr 2010

!      print *, 'kay2 disper= ', kay2

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

! Additions for ground reflection calculations
      azz = 1.d0                                                   !  5 Aug 2016
      bzz = 0.d0                                                   !  5 Aug 2016
      czz = -1.d0                                                  !  5 Aug 2016

! Additions for ground reflection calculations with Coriolis terms ! 11 Jul 2017
! For ground reflection, we begin by taking the simple case where the ground is
! horizontal.  Also, we assume that the wind is parallel to the surface, so 
! we take vr=0 at the surface.  That means we neglect vr in the following 
! formulas.                                                        ! 28 Jul 2017
! Coefficients of the quartic for kr
! ddz(0)+ddz(1)*kr+ddz(2)*dr**2+ddz(3)*kr**3+ddz(4)*kr**4=0        !  5 Jul 2017
      kperpVperp = kth*vth+kph*vph                        
      owiperp=ow-kperpvperp                                    
      owiperp2=owiperp*owiperp      
      ddz(0)=-owiperp2*(kperp2+kA2%value)+OWIperp2*OWIperp2/cs
     &       +kperp2*NS%value        
     &+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                                                        ! 31 Jul 2017
     &+(term1%fph+term2%fph+term3%fph+term4%fph+term5%fph+term6%fph)
     & *kph                                                        ! 31 Jul 2017
     &+(term1%fthth+term2%fthth+term3%fthth+term4%fthth+term5%fthth
     &+term6%fthth)
     & *kth*kth                                                    ! 31 Jul 2017
     &+(term1%fphph+term2%fphph+term3%fphph+term4%fphph+term5%fphph
     &+term6%fphph)
     & *kph*kph                                                    ! 31 Jul 2017
     &+(term1%fthph+term2%fthph+term3%fthph+term4%fthph+term5%fthph
     &+term6%fthph)
     & *kth*kph                                                    ! 31 Jul 2017
      ddz(1)=0.d0
!     &      +2.d0*vr*OWIperp*(kperp2+kA2%value-2.d0*OWIperp2/cs)  ! 28 Jul 2017
     &+term1%fr+term2%fr+term3%fr+term4%fr+term5%fr+term6%fr
     &+(term1%frth+term2%frth+term3%frth+term4%frth+term5%frth
     &+term6%frth)
     & *kth                                                        ! 31 Jul 2017
     &+(term1%frph+term2%frph+term3%frph+term4%frph+term5%frph
     &+term6%frph)
     & *kph                                                        ! 31 Jul 2017
      ddz(2)=-OWIperp2
     & +term1%frr+term2%frr+term3%frr+term4%frr+term5%frr+term6%frr
!     &       -vr*vr*(kperp2+kA2%value)+6.d0*OWIperp2*vr*vr/cs     ! 28 Jul 2017
      ddz(3)=0.d0 
!     &      +2.d0*vr*OWIperp*(1.d0-2.d0*vr*vr/cs)                 ! 28 Jul 2017
      ddz(4)=0.d0
!     &      -vr*vr*(1.d0-vr*vr/cs)                                ! 28 Jul 2017

      azz=-ddz(2)/ow2                                              ! 28 Jul 2017
      bzz=ddz(1)/(2.d0*ow2)                                        ! 28 Jul 2017
      czz=ddz(0)/ow2                                               ! 28 Jul 2017

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


!      kz2a=((ows/cs)-kA2%value)-(kth2+kph2)*(ows-Ns%value)/ows     ! 13 Apr 2011
!      kz2a=((owiperp2/cs)-kA2%value)
!     &     -(kth2+kph2)*(owiperp2-Ns%value)/owiperp2               ! 31 Jul 2017

c      H=(kth**2+kph**2)*(Ns%value-ows) - ows*(kr**2+kA2%value-(ows/cs))
      H1=(kth**2+kph**2)*(Ns%value)                                ! 27 Jul 2010
      H2=(kth**2+kph**2)*(-ows)                                    ! 27 Jul 2010
      H3=-ows*(kr**2)                                              ! 17 Jul 2017
      H4=-ows*(kA2%value)                                          ! 17 Jul 2017
      H5=-ows*(-(ows/cs))                                          ! 17 Jul 2017
      H6 = term1%value                                             ! 11 Jul 2017
c      print*,'Coriolis term 1 = ',term1%value
      H7 = term2%value                                             ! 11 Jul 2017
      H8 = term3%value                                             ! 11 Jul 2017
      H9 = term4%value                                             ! 11 Jul 2017
      H10 = term5%value                                            ! 11 Jul 2017
      H11 = term6%value                                            ! 11 Jul 2017
      H=H1+H2+H3+H4+H5+H6+H7+H8+H9+H10+H11                         ! 11 Jul 2017

      Hrel=H/max(abs(H1),abs(H2),abs(H3),abs(H4),abs(H5)           ! 27 Jul 2010
     & ,abs(H6),abs(H7),abs(H8),abs(H9),abs(H10),abs(H11))         ! 11 Jul 2017

!      kz2a=((owiperp2/cs-kA2%value)*owiperp2
!     &     -(kth2+kph2)*(owiperp2-Ns%value)
!     &      +kth*kth*term2%fthth+kph*kph*term2%fphph 
!     &      +H6+H9+H10+H11)/(owiperp2-term2%frr-term3%frr)         ! 28 Jun 2018
      kz2a = 0.d0                                                  !  6 Jul 2018
      kz2ai = 0.d0                                                 ! 16 Mar 2006
      kzb = 0.d0                                                   ! 16 Mar 2006
      kzbi = 0.d0                                                  ! 16 Mar 2006
      if(azz.eq.0.d0.and.bzz.eq.0.d0) return                       !  6 Jul 2018
      if(azz.eq.0.d0) then                                         !  6 Jul 2018
      kz2a = czz**2/(4.d0*bzz**2)                                  
      return                                                       
      endif                                                        !  6 Jul 2018
      rad = bzz**2 + azz*czz                                       !  6 Jul 2018
      if(rad.ge.0.d0) then                                         !  6 Jul 2018
      kz2a = (bzz**2+rad+2.d0*bzz*sqrt(rad))/azz**2
      else
      kz2a = (bzz**2+rad)/azz**2
      kz2ai = 2.d0*bzz*sqrt(-rad)/azz**2
      Endif                                                        !  6 Jul 2018

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)*owi                             ! 29 Apr 2010
     &       -2.d0*(kr**2+kA2%value)*owi+4*owi*ows/cs              ! 29 Apr 2010
     &       +term1%pom+term2%pom+term3%pom+term4%pom+term5%pom    ! 11 Jul 2017
     &       +term6%pom                                            ! 11 Jul 2017
      pHkr = -2.d0*ows*kr
     & -pHow*VR                                                    ! 14 Apr 2010
     &       +term1%pkr+term2%pkr+term3%pkr+term4%pkr+term5%pkr    ! 11 Jul 2017
     &       +term6%pkr                                            ! 11 Jul 2017
      pHkth = 2.d0*(Ns%value-ows)*kth 
     & -pHow*VTH                                                   ! 14 Apr 2010
     &     +term1%pkth+term2%pkth+term3%pkth+term4%pkth+term5%pkth ! 11 Jul 2017
     &     +term6%pkth                                             ! 11 Jul 2017
      pHkph = 2.d0*(Ns%value-ows)*kph 
     & -pHow*VPH                                                   ! 14 Apr 2010
     &     +term1%pkph+term2%pkph+term3%pkph+term4%pkph+term5%pkph ! 11 Jul 2017
     &     +term6%pkph                                             ! 11 Jul 2017
      kpHk = -2.d0*ows*ks + 2.d0*Ns%value*(kth**2+kph**2)
     & -pHow*KV                                                    ! 14 Apr 2010
     &     +kr*(term1%pkr+term2%pkr+term3%pkr+term4%pkr+term5%pkr  ! 11 Jul 2017
     &     +term6%pkr)                                             ! 11 Jul 2017
     &+kth*(term1%pkth+term2%pkth+term3%pkth+term4%pkth+term5%pkth ! 11 Jul 2017
     &     +term6%pkth)                                            ! 11 Jul 2017
     &+kph*(term1%pkph+term2%pkph+term3%pkph+term4%pkph+term5%pkph ! 11 Jul 2017
     &     +term6%pkph)                                            ! 11 Jul 2017

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
     & -pHow*(KR*PVRT+KTH*PVTHT+KPH*PVPHT)                         ! 14 Apr 2010
     & +term1%pt+term2%pt+term3%pt+term4%pt+term5%pt+term6%pt      ! 11 Jul 2017
      pHr = (kth**2+kph**2)*Ns%pr - ((ows/cs)**2)*pcsr - ows*kA2%pr
     & -pHow*(KR*PVRR+KTH*PVTHR+KPH*PVPHR)                         ! 14 Apr 2010
     & +term1%pr+term2%pr+term3%pr+term4%pr+term5%pr+term6%pr      ! 11 Jul 2017
      pHth = (kth**2+kph**2)*Ns%pth - ((ows/cs)**2)*pcsth - ows*kA2%pth
     & -pHow*(KR*PVRTH+KTH*PVTHTH+KPH*PVPHTH)                      ! 14 Apr 2010
     & +term1%pth+term2%pth+term3%pth+term4%pth+term5%pth+term6%pth! 11 Jul 2017
      pHph = (kth**2+kph**2)*Ns%pph - ((ows/cs)**2)*pcsph - ows*kA2%pph
     & -pHow*(KR*PVRPH+KTH*PVTHPH+KPH*PVPHPH)                      ! 14 Apr 2010
     & +term1%pph+term2%pph+term3%pph+term4%pph+term5%pph+term6%pph! 11 Jul 2017

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

!      print *, 'gwwnl n3 = ', n3, 'n4 = ', n4
!     & , 'nbig = ', nbig, 'nsmall = ', nsmall

c      errRel=1.d10                                                ! 27 Jul 2010
c      IF (KAY2.NE.0.d0) errRel =ks/KAY2-1.d0        ! 27 Jul 2010 ! 19 Nov 2010
      errRel = errorcalc(ks,kay2)                                  ! 11 Mar 2011
C ----------------------------------------------------------------------!-------

      DO j=1,4                                                     ! 19 Nov 2010
      errRell(j) = errorcalc(ks,kaysq(j))                          ! 11 Mar 2011
      errRelb(j) = errorcalc(ks,kaysqb(j))                         ! 11 Mar 2011
      errRela(j) = errorcalc(ks,kaysqa(j))                         ! 11 Mar 2011
C ----------------------------------------------------------------------!-------
      ENDDO                                                        ! 19 Nov 2010

! 14 Apr 2010      ENTRY WINDR                                     ! 27 Feb 2006
      RETURN
C ----------------------------------------------------------------------!-------
      END SUBROUTINE GWWNL 
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

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

      EQUIVALENCE (w(1),earthr)
     & ,(RAY,W(2))                                                 ! 28 Jun 2010

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

C     COMMON DECK "FLAG" INSERTED HERE
      LOGICAL NEWWR,NEWWP,NEWTRC,PENET
      COMMON /FLG/ NTYP,NEWWR,NEWWP,NEWTRC,PENET,LINES,IHOP,HPUNCH
      COMMON/FLGP/NSET

! 6 Aug 2010      ntyp = ray                                       !  5 Aug 2010

      CALL iwindr                                                  ! 14 Apr 2010
      CALL ispeed
      CALL irecvr   
      CALL itopog   
      CALL isurfac
      CALL ipresg(modpg)                                           ! 24 Jan 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

! 14 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
      character*8 type11 ,type21 ,type12 ,type22 ,type13 ,type23   !  5 Aug 2010
      REAL(KIND=8),DIMENSION(3) :: tbpe
      LOGICAL :: space

C     COMMON DECK "FLAG" INSERTED HERE
      LOGICAL NEWWR,NEWWP,NEWTRC,PENET
      COMMON /FLG/ NTYP,NEWWR,NEWWP,NEWTRC,PENET,LINES,IHOP,HPUNCH
      COMMON/FLGP/NSET

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

      DATA dlspm /'GWWNL'/
      DATA tbpe /3*'1'/                                            ! 14 Apr 2010

      type13 = 'ordinary'                                          !  9 Mar 2012
      type23 = ' wave   '                                          !  9 Mar 2012
      type12 = '       X'                                          !  9 Mar 2012
      type22 = ' wave   '                                          !  9 Mar 2012
      type11 = 'large X '                                          !  9 Mar 2012
      type21 = 'wave    '                                          !  9 Mar 2012
      CALL SET2(RAYNAME(1,1),type11, 1)                            !  5 Aug 2010
      CALL SET2(RAYNAME(2,1),type21, 1)                            !  5 Aug 2010
      CALL SET2(RAYNAME(1,2),type12, 1)                            !  5 Aug 2010
      CALL SET2(RAYNAME(2,2),type22, 1)                            !  5 Aug 2010
      CALL SET2(RAYNAME(1,3),type13, 1)                            !  5 Aug 2010
      CALL SET2(RAYNAME(2,3),type23, 1)                            !  5 Aug 2010
      dispm = dlspm
      CALL rmove(type,tbpe,3)
      modrin = 'ACOUSTIC GRAVITY WAVE --- WITH 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 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

      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 QUART(dd,sol,soli,Nsol
     & ,sola,solai,solb,solbi
     & ,test,testa,testb                                           ! 24 Jun 2010
     & ,relerr,relerra,relerrb                                     ! 28 Jun 2010
     & ,n3,n4,nbig,nsmall                                          ! 15 Jul 2010
     & ,height,ks,cs,c,error,errorsave
     & ,kaysq,kaysqa,kaysqb
     & ,kay2,drdt20                                                ! 15 Jun 2010
     & ,kv,propagate)                                              ! 24 May 2013
C ----------------------------------------------------------------------!-------
! Nsol is the number of real solutions of the quartic. kv is k dot v
! propagate(i) is zero if the ith root is evanescent, non-zero if propagating
      IMPLICIT double precision (a-h,o-z)
      dimension dd(0:4),sol(4),soli(4)
     & ,solb(4),solbi(4),sola(4),solai(4)                          ! 15 Jun 2010
      complex*16 test(4),testa(4),testb(4)                         ! 24 Jun 2010
      REAL(KIND=8) :: KAY2,KAY2I,kay 
     & , kaysq(4) , kaysqb(4) , kaysqa(4)                          ! 14 Jun 2010
      REAL(KIND=8) :: ks,kv,kperp2
      dimension nn3(4,4), nn4(4,4)
!      double precision D1mach
      complex*16 p4,relerr(4),relerra(4),relerrb(4),xxx
      integer propagate(4)                                         ! 24 May 2013

      data nn3 /0,3,2,2,3,0,1,1,2,1,0,1,2,1,1,0/                   ! 16 Jun 2010
      data nn4 /0,4,4,3,4,0,4,4,4,4,0,3,3,3,2,0/                   ! 16 Jun 2010

      if(drdt20.eq.23.d0) print *, ' dd = ', dd
      call quartic(dd,sol,soli,Nsol)

      do i=1,4
         sola(i) = sol(i)
         solai(i) = soli(i)
      enddo

      call identifyroot(
     & Nsol,kv,                                                    ! 24 May 2013
     & dd,sol,soli,n3a,n4a,nbiga,nsmalla                           ! 28 Jul 2010
     & ,nsmall2,propagate)                                         ! 24 May 2013

      call identifyroot2(Nsol,kv,dd,sol,soli,n3,n4,nbig,nsmall 
     & ,nsmall2,propagate)                                         ! 28 May 2013
c      IF(n3a.ne.n3) 
c     & print *, 'identification differs, n3 = ', n3, ' n3a = ', n3a
c     &        , ' sol = ', sol
c     & , 'soli = ', soli, ' Nsol = ', Nsol, ' kv = ', kv

!      IF(n4a.ne.n4) 
!     & print *, 'identification differs, n4 = ', n4, ' n4a = ', n4a
!     &        , ' sol = ', sol
c      IF(nsmalla.ne.nsmall) 
c     & print *, 'identification differs, nsmall = '
c     & , nsmall, ' nsmalla = ', nsmalla, ' sol = ', sol
c     & , 'soli = ', soli, ' Nsol = ', Nsol, ' kv = ', kv
c      IF(nbiga.ne.nbig) 
c     & print *, 'identification differs, nbig = '
c     & , nbig, ' nbiga = ', nbiga, ' sol = ', sol
c     & , 'soli = ', soli, ' Nsol = ', Nsol, ' kv = ', kv

      if(drdt20.eq.23.d0) print *, ' Nsol = ', Nsol,'height = ',height
!      if(drdt20.eq.23.d0) print *, ' sol = ', sol
!      if(drdt20.eq.23.d0) print *, ' soli = ', soli
      if(drdt20.eq.23.d0) print *, 'errorsave = ', errorsave
      if(drdt20.eq.23.d0) print *, 'error = ', error
      if(drdt20.eq.23.d0) print *, 'n3 = ', n3, 'n4 = ', n4
     & , 'nbig = ', nbig, 'nsmall = ', nsmall

C ----------------------------------------------------------------------!-------
C ----------------------------------------------------------------------!-------
      IF(nbig.ne.0) call bigroot(dd,sol,soli,Nsol
     & ,sola,solai,solb,solbi
     & ,test,testa,testb                                          
     & ,relerr,relerra,relerrb                                     ! 28 Jun 2010
     & ,nbig,nsmall,n3,n4                                         
     & ,height,ks,cs,c,error,errorsave
     & ,kaysq,kaysqa,kaysqb
     & ,kay2,drdt20)                                              
C ----------------------------------------------------------------------!-------
      IF(n4.ne.0) call root4(dd,sol,soli,Nsol
     & ,sola,solai,solb,solbi
     & ,test,testa,testb                                          
     & ,relerr,relerra,relerrb                                     ! 28 Jun 2010
     & ,nbig,nsmall,n3,n4                                         
     & ,height,ks,cs,c,error,errorsave
     & ,kaysq,kaysqa,kaysqb
     & ,kay2,drdt20)                                              
C ----------------------------------------------------------------------!-------
      IF(nsmall.ne.0) call smallroot(dd,sol,soli,Nsol
     & ,sola,solai,solb,solbi
     & ,test,testa,testb                                          
     & ,relerr,relerra,relerrb                                     ! 28 Jun 2010
     & ,nbig,nsmall,n3,n4                                         
     & ,height,ks,cs,c,error,errorsave
     & ,kaysq,kaysqa,kaysqb
     & ,kay2,drdt20)                                              
C ----------------------------------------------------------------------!-------
      IF(nsmall2.ne.0) call smallroot(dd,sol,soli,Nsol             ! 29 May 2013
     & ,sola,solai,solb,solbi
     & ,test,testa,testb                                          
     & ,relerr,relerra,relerrb                                     ! 28 Jun 2010
     & ,nbig,nsmall2,n3,n4                                         
     & ,height,ks,cs,c,error,errorsave
     & ,kaysq,kaysqa,kaysqb
     & ,kay2,drdt20)                                              
C ----------------------------------------------------------------------!-------
      IF(n3.ne.0) call root3(dd,sol,soli,Nsol
     & ,sola,solai,solb,solbi
     & ,test,testa,testb                                          
     & ,relerr,relerra,relerrb                                     ! 28 Jun 2010
     & ,nbig,nsmall,n3,n4                                         
     & ,height,ks,cs,c,error,errorsave
     & ,kaysq,kaysqa,kaysqb
     & ,kay2,drdt20)                                              
C ----------------------------------------------------------------------!-------
C ----------------------------------------------------------------------!-------

!      if(n3.eq.0) then
!      print*, 'power =', power, 'powdif=',powdif
!     & , 'powdifmin=',powdifmin
!      stop
!      endif

!      print *, 'quart n3 = ', n3, 'n4 = ', n4
!     & , 'nbig = ', nbig, 'nsmall = ', nsmall
      return
C ----------------------------------------------------------------------!-------
      end subroutine quart
C ----------------------------------------------------------------------!-------
C ----------------------------------------------------------------------!-------
      subroutine identifyroot(
     & Nsol,kv,                                                    ! 24 May 2013
     & dd,sol,soli,n3,n4,nbig,nsmall                               ! 28 Jul 2010
     & ,nsmall2,propagate)                                         ! 24 May 2013
C ----------------------------------------------------------------------!-------
! Nsol is the number of real solutions of the quartic. kv is k dot v
! propagate(i) is zero if the ith root is evanescent, non-zero if propagating
      IMPLICIT double precision (a-h,o-z)
      dimension dd(0:4),sol(4),soli(4)
     & ,solb(4),solbi(4),sola(4),solai(4)                          ! 15 Jun 2010
      double precision power(4), powdif(4,4)                       ! 15 Jul 2010
     & ,prod(4,4)                                                  ! 30 Jul 2010
      REAL(KIND=8) :: kv                                           ! 24 May 2013
      integer propagate(4)                                         ! 24 May 2013
      powdifmax=0.d0                                               ! 15 Jul 2010
      do i=1,4
         power(i) = sol(i)**2+soli(i)**2                           ! 15 Jul 2010
         powdifmax=powdifmax+power(i)                              ! 15 Jul 2010
         powdifmin=powdifmax                                       ! 15 Jul 2010
      enddo

         nbig = 1
         n3 = 2
         n4 = 3
         nsmall = 4

! Find the largest of the roots                                    ! 30 Jul 2010
      powmax=0.d0                                                  ! 15 Jul 2010
      do i=1,4                                                    
         if(power(i).gt.powmax) then                               ! 30 Jul 2010
            nbig=i                                                 ! 15 Jul 2010
            powmax=power(i)                                        ! 15 Jul 2010
         endif
      enddo                                                       

      if(    nbig.eq.1.and.abs(dd(4)).lt.2.d-7
     &   .or.nbig.eq.2.and.abs(dd(4)).lt.1.d-6                     !  3 Aug 2010
     &   .or.nbig.eq.3.and.abs(dd(4)).lt.1.d-5                     !  3 Aug 2010
     &   .or.nbig.eq.4.and.abs(dd(4)).lt.7.d-6) then               !  9 Aug 2010

! Find the pair of small real or complex roots                     ! 15 Jul 2010
      n3 = 1                                                       ! 15 Jul 2010
      n4 = 1                                                       ! 15 Jul 2010
      do i=1,3                                                     ! 26 Jul 2010
         do j=i+1,4
            powdif(i,j) = abs(power(i)-power(j))                   ! 15 Jul 2010
            prod(i,j) = sol(i)*sol(j)+soli(i)*soli(j)              ! 30 Jul 2010
            if(powdif(i,j).lt.powdifmin
     &      .and.prod(i,j).lt.0.d0                  ! 30 Jul 2010  !  9 Aug 2010
     &      ) then                                   
               n3 = i                                              ! 15 Jul 2010
               n4 = j                                              ! 15 Jul 2010
               powdifmin=powdif(i,j)                               ! 30 Jul 2010
            endif
!            if(powdif(i,j).lt.powdifmin) powdifmin=powdif(i,j)    ! 30 Jul 2010
         enddo
      enddo

! Find the other root                                              ! 30 Jul 2010
      do i=1,4                                                    
         if(i.ne.n3.and.i.ne.n4.and.i.ne.nbig) nsmall=i            ! 15 Jul 2010
      enddo                                                       

      elseif(nbig.eq.1) then                                       ! 30 Jul 2010

         n3 = 2
         n4 = 3
         nsmall = 4

      elseif(nbig.eq.2) then                                       !  3 Aug 2010

         n3 = 3
         n4 = 4
         nsmall = 1

      elseif(nbig.eq.3) then                                       !  3 Aug 2010

         n3 = 1
         n4 = 2
         nsmall = 4

      else                                                         ! 23 Jul 2010

         n3 = 2                                                    ! 10 Aug 2010
         n4 = 3                                                    !  9 Aug 2010
         nsmall = 1                                                ! 10 Aug 2010

!         n3 = 1                                                   ! 29 Nov 2010
!         nsmall = 2                                               ! 29 Nov 2010

      endif                                                        ! 23 Jul 2010

      return
C ----------------------------------------------------------------------!-------
      end subroutine identifyroot
C ----------------------------------------------------------------------!-------
C ----------------------------------------------------------------------!-------

      subroutine identifyroot2(Nsol,kv,dd,sol,soli,n3,n4,nbig,nsmall
     & ,nsmall2,propagate)                                         ! 28 May 2013
C ----------------------------------------------------------------------!-------
! The identification of roots of the quartic is still not perfect  !  1 Jun 2013
! Nsol is the number of real solutions of the quartic. kv is k dot v
! propagate(i) is zero if the ith root is evanescent, non-zero if propagating
      IMPLICIT double precision (a-h,o-z)
      dimension dd(0:4),sol(4),soli(4)
     & ,solb(4),solbi(4),sola(4),solai(4)
      double precision power(4), powdif(4,4)
     & ,prod(4,4)
      REAL(KIND=8) :: kv
      integer propagate(4),m(4)
      n3=0
      n4=0
      nsmall=0
      nsmall2=0
      nbig=0
      DO i=1,4
         propagate(i)=0
      ENDDO
C ----------------------------------------------------------------------!-------
!      Ordering roots
      m(1)=1
      m(2)=2
      m(3)=3
      m(4)=4
      IF(sol(1).le.sol(2).and.sol(2).le.sol(3).and.sol(3).le.sol(4))
     & GOTO 1
      m(1)=1
      DO i=1,4
         IF(sol(i).lt.sol(m(1))) m(1)=i
      ENDDO
      m(2)=1
      IF(m(2).eq.m(1)) m(2)=2
      DO i=1,4
         IF(i.ne.m(1).and.sol(i).lt.sol(m(2))) m(2)=i
      ENDDO
      m(3)=1
      IF(m(3).eq.m(1)) m(3)=2
      IF(m(3).eq.m(2)) m(3)=3
      DO i=1,4
         IF(i.ne.m(1).and.i.ne.m(2).and.sol(i).lt.sol(m(3))) m(3)=i
      ENDDO
      DO i=1,4
         IF(i.ne.m(1).and.i.ne.m(2).and.i.ne.m(3)) m(4)=i
      ENDDO
  1   CONTINUE
C ----------------------------------------------------------------------!-------
      IF(sol(m(1)).gt.sol(m(2)).or.sol(m(2)).gt.sol(m(3))
     & .or.sol(m(3)).gt.sol(m(4))) THEN
         print *, 'Roots out of order, sol = ', 
     &   sol(m(1)), sol(m(2)), sol(m(3)), sol(m(4))
     &   , soli(m(1)), soli(m(2)), soli(m(3)), soli(m(4))
     &   , ' m(1) = ', m(1), ' m(2) = ', m(2)
     &   , ' m(3) = ', m(3), ' m(4) = ', m(4) 
         return
      ENDIF
C ----------------------------------------------------------------------!-------
c      IF(m(4).eq.2) THEN
c         print *, 'm(4) = 2, sol = ', 
c     &   sol(m(1)), sol(m(2)), sol(m(3)), sol(m(4))
c     &   , soli(m(1)), soli(m(2)), soli(m(3)), soli(m(4))
c     &   , ' m(1) = ', m(1), ' m(2) = ', m(2)
c     &   , ' m(3) = ', m(3), ' m(4) = ', m(4)
c     &   , ' kv = ', kv
c      ENDIF
C ----------------------------------------------------------------------!-------
      IF(Nsol.ne.2.and.Nsol.ne.4) THEN
c         print *, ' Nsol is not 2 or 4, Nsol = ', Nsol
c         print *, ' dd = ', dd
c         print *, ' sol = ', sol
c         print *, ' soli = ', soli
         IF(Nsol.ne.0) return
      ENDIF
C ----------------------------------------------------------------------!-------
      IF(kv.lt.0.) THEN
!        Upwind case
         IF(Nsol.eq.4) THEN
!           case of 4 propagating roots
            nbig=m(1)
            n4=m(2)
            n3=m(3)
            nsmall=m(4)
            nsmall2=0
            DO i=1,4
               propagate(i)=1
            ENDDO
         ELSEIF(Nsol.eq.2) THEN
!           case of 2 propagating roots and 2 evanescent roots
C ----------------------------------------------------------------------!-------
!            test1 = dd(3)*dd(1)**3/(dd(0)*dd(2)**3)*1.d10          ! 30 May 2013
!            test2 = dd(0)*dd(3)**2/dd(2)**3*1.d20                  ! 30 May 2013
!            test3 = dd(4)*dd(1)**4/(dd(0)*dd(2)**4)*1.d10          ! 30 May 2013
!            test4 = dd(4)*dd(0)/dd(2)**2*1.d10                     ! 30 May 2013
!            IF(abs(test1).lt.1.d0.and.abs(test2).lt.1.d0
!     &         .and.abs(test3).lt.1.d0.and.abs(test4).lt.1.d0) THEN
!               print *, ' Nsol is 2 special case'
!               print *, ' dd = ', dd
!               print *, ' sol = ', sol
!               print *, ' soli = ', soli
!               print *, ' kv = ', kv
!            ENDIF
C ----------------------------------------------------------------------!-------
! Find the 2 real roots in increasing order                        ! 31 May 2013
            DO i=1,4
               IF(soli(m(i)).eq.0.d0) THEN
                  nr1=m(i)
                  inext=i+1
                  GO TO 2
               ENDIF
            ENDDO
   2        CONTINUE
            DO i=inext,4
               IF(soli(m(i)).eq.0.d0) THEN
                  nr2=m(i)
                  GO TO 3
               ENDIF
            ENDDO
   3        CONTINUE
C ----------------------------------------------------------------------!-------
            IF(sol(nr1)*sol(nr2).gt.0.d0) THEN                     ! 31 May 2013
! case of both real roots having the same sign
               nbig=m(1)
               nsmall2=m(2)
               n3=m(3)
               nsmall=m(4)
               n4=0
               propagate(m(1))=1
               propagate(m(2))=1
               propagate(m(3))=0
               propagate(m(4))=0
            ELSE
! case of both real roots having the opposite sign                 ! 31 May 2013
               nbig=nr1
               nsmall=nr2
               n3=m(3)
               n4=m(2)
               nsmall2=0
               propagate(m(1))=1
               propagate(m(2))=1
               propagate(m(3))=1
               propagate(m(4))=1
            ENDIF
         ELSE
!           case of 4 evanescent roots (Nsol=0)                    ! 31 May 2013
               nbig=0
               nsmall=0
               n3=0
               n4=0
               nsmall2=0
               propagate(m(1))=0
               propagate(m(2))=0
               propagate(m(3))=0
               propagate(m(4))=0
         ENDIF
      ELSE
C ----------------------------------------------------------------------!-------
!        Downwind case
         IF(Nsol.eq.4) THEN
!           case of 4 propagating roots
            nsmall=m(1)
            n4=m(2)
            n3=m(3)
            nbig=m(4)
            nsmall2=0
            DO i=1,4
               propagate(i)=1
            ENDDO
         ELSEIF(Nsol.eq.2) THEN
!           case of 2 propagating roots and 2 evanescent roots
C ----------------------------------------------------------------------!-------
c            test1 = dd(3)*dd(1)**3/(dd(0)*dd(2)**3)*1.d10          ! 30 May 2013
c            test2 = dd(0)*dd(3)**2/dd(2)**3*1.d20                  ! 30 May 2013
c            test3 = dd(4)*dd(1)**4/(dd(0)*dd(2)**4)*1.d10          ! 30 May 2013
c            test4 = dd(4)*dd(0)/dd(2)**2*1.d10                     ! 30 May 2013
c            IF(abs(test1).lt.1.d0.and.abs(test2).lt.1.d0
c     &         .and.abs(test3).lt.1.d0.and.abs(test4).lt.1.d0) THEN
c               print *, ' Nsol is 2 special case'
c               print *, ' dd = ', dd
c               print *, ' sol = ', sol
c               print *, ' soli = ', soli
c               print *, ' kv = ', kv
c            ENDIF
C ----------------------------------------------------------------------!-------
! Find the 2 real roots in increasing order                        ! 31 May 2013
            DO i=1,4
               IF(soli(m(i)).eq.0.d0) THEN
                  nr1=m(i)
                  inext=i+1
                  GO TO 4
               ENDIF
            ENDDO
   4        CONTINUE
            DO i=inext,4
               IF(soli(m(i)).eq.0.d0) THEN
                  nr2=m(i)
                  GO TO 5
               ENDIF
            ENDDO
   5        CONTINUE
C ----------------------------------------------------------------------!-------
            IF(sol(nr1)*sol(nr2).gt.0.d0) THEN                     ! 31 May 2013
! case of both real roots having the same sign
               nsmall2=m(1)
               n4=m(2)
               nsmall=m(3)
               nbig=m(4)
               n3=0
               n4=0
               n3=m(2)
               propagate(m(1))=0
               propagate(m(2))=0
               propagate(m(3))=1
               propagate(m(4))=1
            ELSE
! case of both real roots having the opposite sign                 ! 31 May 2013
               nbig=nr2
               nsmall=nr1
               n3=m(3)
               n4=m(2)
               nsmall2=0
               propagate(m(1))=1
               propagate(m(2))=1
               propagate(m(3))=1
               propagate(m(4))=1
            ENDIF
         ELSE
!           case of 4 evanescent roots (Nsol=0)                    ! 31 May 2013
               nbig=0
               nsmall=0
               n3=0
               n4=0
               nsmall2=0
               propagate(m(1))=0
               propagate(m(2))=0
               propagate(m(3))=0
               propagate(m(4))=0
         ENDIF
      ENDIF
      return
C ----------------------------------------------------------------------!-------
      end subroutine identifyroot2
C ----------------------------------------------------------------------!-------
C ----------------------------------------------------------------------!-------

      SUBROUTINE bigroot(dd,sol,soli,Nsol
     & ,sola,solai,solb,solbi
     & ,test,testa,testb                                          
     & ,relerr,relerra,relerrb                                     ! 28 Jun 2010
     & ,nbig,nsmall,n3,n4                                         
     & ,height,ks,cs,c,error,errorsave
     & ,kaysq,kaysqa,kaysqb
     & ,kay2,drdt20)                                              
C ----------------------------------------------------------------------!-------
      IMPLICIT double precision (a-h,o-z)
      dimension dd(0:4),sol(4),soli(4)
     & ,solb(4),solbi(4),sola(4),solai(4)                          ! 15 Jun 2010
      complex*16 test(4),testa(4),testb(4)                         ! 24 Jun 2010
      REAL(KIND=8) :: KAY2,KAY2I,kay 
     & , kaysq(4) , kaysqb(4) , kaysqa(4)                          ! 14 Jun 2010
      REAL(KIND=8) :: ks,kv,kperp2
      dimension nn3(4,4), nn4(4,4)
!      double precision D1mach
      double precision errorcalc                                   ! 11 Mar 2011
      external errorcalc                                           ! 11 Mar 2011
      complex*16 p4,relerr(4),relerra(4),relerrb(4),xxx

C ----------------------------------------------------------------------!-------
      IF(nbig.eq.0) return                                         ! 29 May 2013
      sola(nbig) = sol(nbig)                                       ! 28 Jun 2010
      solai(nbig) = soli(nbig)                                     ! 28 Jun 2010
!      testa(nbig) = dd(0) + cmplx(sola(nbig),solai(nbig))*(
!     &        dd(1) + cmplx(sola(nbig),solai(nbig))*(
!     &        dd(2) + cmplx(sola(nbig),solai(nbig))*(
!     &        dd(3) + cmplx(sola(nbig),solai(nbig))*dd(4))))
!      print *, 'test 1 = ', testa(nbig)

      xxx = cmplx(sola(nbig),solai(nbig))
      testa(nbig) = p4(dd,xxx,relerra(nbig))
!      print *, 'test 2 = ', testa(nbig)
!      print *, 'rel err = ', relerr(nbig)

      if(drdt20.eq.23.d0) print *, 'root 1 testa(nbig) = ', testa(nbig)
      kay2 = (sola(nbig)**2-solai(nbig)**2)/cs
      kay2i=2.d0*sola(nbig)*solai(nbig)/cs
      error=errorcalc(ks,kay2)                                     ! 11 Mar 2011
      if(drdt20.eq.23.d0) print *, 'root 1 error = ', error
      if(drdt20.eq.23.d0) print *, 'root 1 kay2 = ', kay2
     &                       , 'sola(nbig) = ', sola(nbig), solai(nbig)
      if(drdt20.eq.23.d0) print *, 'root 1 kay2i = ', kay2i
      kaysq(nbig)=kay2                                             ! 11 Jun 2010
      kaysqa(nbig)=kay2                                            ! 14 Jun 2010
      test(nbig)=testa(nbig)                                       ! 28 Jun 2010
        
C ----------------------------------------------------------------------!-------
      start = -dd(3)/(2.d0*dd(4))                                  ! 18 Jun 2010
      rad = start*start + (-dd(2)
     &      )/dd(4)                                                ! 18 Jun 2010
      if(rad.lt.0.d0) then
         solb(nbig) = start                    
         solbi(nbig) = sign(1.d0,start)*dsqrt(-rad)

      testb(nbig) = dd(0) + cmplx(solb(nbig),solbi(nbig))*(
     &        dd(1) + cmplx(solb(nbig),solbi(nbig))*(
     &        dd(2) + cmplx(solb(nbig),solbi(nbig))*(
     &        dd(3) + cmplx(solb(nbig),solbi(nbig))*dd(4))))

      xxx = cmplx(solb(nbig),solbi(nbig))

      if(isnan(abs(xxx))) print *, 'xxx = ', xxx, 'dd = ', dd

      testb(nbig) = p4(dd,xxx,relerrb(nbig))

      if(drdt20.eq.23.d0) print *, 'root 1b testb(nbig) = ', testb(nbig)
      kay2 = (solb(nbig)**2-solbi(nbig)**2)/cs
      kay2i=2.d0*solb(nbig)*solbi(nbig)/cs
      error=errorcalc(ks,kay2)                                     ! 11 Mar 2011
      if(drdt20.eq.23.d0) print *, 'root 1b error = ', error
      if(drdt20.eq.23.d0) print *, 'root 1b kay2 = ', kay2
     &              , 'solb(nbig) = ', solb(nbig), solbi(nbig)
      if(drdt20.eq.23.d0) print *, 'root 1b kay2i = ', kay2i
        
      else
   
      solb(nbig) = sign(1.d0,start)*(dabs(start)+dsqrt(rad))       ! 18 Jun 2010
      solbi(nbig) = 0.d0

      testb(nbig) = dd(0) + cmplx(solb(nbig),solbi(nbig))*(
     &        dd(1) + cmplx(solb(nbig),solbi(nbig))*(
     &        dd(2) + cmplx(solb(nbig),solbi(nbig))*(
     &        dd(3) + cmplx(solb(nbig),solbi(nbig))*dd(4))))

      xxx = cmplx(solb(nbig),solbi(nbig))

      if(isnan(abs(xxx))) print *, 'xxx = ', xxx, 'dd = ', dd

      testb(nbig) = p4(dd,xxx,relerrb(nbig))

      if(drdt20.eq.23.d0) print *, 'root 1b testb(nbig) = ', testb(nbig)
      kay2 = (solb(nbig)**2-solbi(nbig)**2)/cs
      kay2i=2.d0*solb(nbig)*solbi(nbig)/cs
      error=errorcalc(ks,kay2)                                     ! 11 Mar 2011
      if(drdt20.eq.23.d0) print *, 'root 1b error = ', error
      if(drdt20.eq.23.d0) print *, 'root 1b kay2 = ', kay2
     &              , 'solb(nbig) = ', solb(nbig), solbi(nbig)
      if(drdt20.eq.23.d0) print *, 'root 1b kay2i = ', kay2i
        
      rad = start*start + (-dd(2)
     &      +(dd(1)+dd(0)/solb(nbig))/solb(nbig)
     &      )/dd(4)                                                ! 18 Jun 2010

      solb(nbig) = sign(1.d0,start)*(dabs(start)+dsqrt(rad))       ! 18 Jun 2010
      solbi(nbig) = 0.d0

      testb(nbig) = dd(0) + cmplx(solb(nbig),solbi(nbig))*(
     &        dd(1) + cmplx(solb(nbig),solbi(nbig))*(
     &        dd(2) + cmplx(solb(nbig),solbi(nbig))*(
     &        dd(3) + cmplx(solb(nbig),solbi(nbig))*dd(4))))

      xxx = cmplx(solb(nbig),solbi(nbig))

      if(isnan(abs(xxx))) print *, 'xxx = ', xxx, 'dd = ', dd

      testb(nbig) = p4(dd,xxx,relerrb(nbig))

      if(drdt20.eq.23.d0) print *, 'root 1b testb(nbig) = ', testb(nbig)
      kay2 = (solb(nbig)**2-solbi(nbig)**2)/cs
      kay2i=2.d0*solb(nbig)*solbi(nbig)/cs
      error=errorcalc(ks,kay2)                                     ! 11 Mar 2011
      if(drdt20.eq.23.d0) print *, 'root 1b error = ', error
      if(drdt20.eq.23.d0) print *, 'root 1b kay2 = ', kay2
     &               , 'solb(nbig) = ', solb(nbig), solbi(nbig)
      if(drdt20.eq.23.d0) print *, 'root 1b kay2i = ', kay2i
        
      rad = start*start + (-dd(2)
     &      +(dd(1)+dd(0)/solb(nbig))/solb(nbig)
     &      )/dd(4)                                                ! 18 Jun 2010

      solb(nbig) = sign(1.d0,start)*(dabs(start)+dsqrt(rad))       ! 18 Jun 2010
      solbi(nbig) = 0.d0

      testb(nbig) = dd(0) + cmplx(solb(nbig),solbi(nbig))*(
     &        dd(1) + cmplx(solb(nbig),solbi(nbig))*(
     &        dd(2) + cmplx(solb(nbig),solbi(nbig))*(
     &        dd(3) + cmplx(solb(nbig),solbi(nbig))*dd(4))))

      xxx = cmplx(solb(nbig),solbi(nbig))

      if(isnan(abs(xxx))) print *, 'xxx = ', xxx, 'dd = ', dd

      testb(nbig) = p4(dd,xxx,relerrb(nbig))

      if(drdt20.eq.23.d0) print *, 'root 1b testb(nbig) = ', testb(nbig)
      kay2 = (solb(nbig)**2-solbi(nbig)**2)/cs
      kay2i=2.d0*solb(nbig)*solbi(nbig)/cs
      error=errorcalc(ks,kay2)                                     ! 11 Mar 2011
      if(drdt20.eq.23.d0) print *, 'root 1b error = ', error
      if(drdt20.eq.23.d0) print *, 'root 1b kay2 = ', kay2
     &              , 'solb(nbig) = ', solb(nbig), solbi(nbig)
      if(drdt20.eq.23.d0) print *, 'root 1b kay2i = ', kay2i
        
      endif                                                        ! 18 Jun 2010

      kaysqb(nbig)=kay2                                            ! 11 Jun 2010

C ----------------------------------------------------------------------!-------
      relerr(nbig) = relerra(nbig)                                 ! 28 Jun 2010
      if(.not.(isnan(abs(testb(nbig)))).and..not.(isnan(kaysqb(nbig)))
     & .and.abs(testb(nbig)).lt.abs(testa(nbig)) 
     &   ) then
         kaysq(nbig) = kaysqb(nbig)                                ! 24 Jun 2010
         test(nbig)=testb(nbig)                                    ! 28 Jun 2010
         relerr(nbig) = relerrb(nbig)                              ! 28 Jun 2010
      endif

      kay2=kaysq(nbig)
      kay2i=0.d0

C ----------------------------------------------------------------------!-------
      end subroutine bigroot
C ----------------------------------------------------------------------!-------
C ----------------------------------------------------------------------!-------
      SUBROUTINE root3(dd,sol,soli,Nsol
     & ,sola,solai,solb,solbi
     & ,test,testa,testb                                          
     & ,relerr,relerra,relerrb                                     ! 28 Jun 2010
     & ,nbig,nsmall,n3,n4                                         
     & ,height,ks,cs,c,error,errorsave
     & ,kaysq,kaysqa,kaysqb
     & ,kay2,drdt20)                                              
C ----------------------------------------------------------------------!-------
      IMPLICIT double precision (a-h,o-z)
      dimension dd(0:4),sol(4),soli(4)
     & ,solb(4),solbi(4),sola(4),solai(4)                          ! 15 Jun 2010
      complex*16 test(4),testa(4),testb(4)                         ! 24 Jun 2010
      REAL(KIND=8) :: KAY2,KAY2I,kay 
     & , kaysq(4) , kaysqb(4) , kaysqa(4)                          ! 14 Jun 2010
      REAL(KIND=8) :: ks,kv,kperp2
      dimension nn3(4,4), nn4(4,4)
!      double precision D1mach
      double precision errorcalc                                   ! 11 Mar 2011
      external errorcalc                                           ! 11 Mar 2011
      complex*16 p4,relerr(4),relerra(4),relerrb(4),xxx

C ----------------------------------------------------------------------!-------
      IF(n3.eq.0) return                                           ! 29 May 2013
      sola(n3) = sol(n3)                                           ! 28 Jun 2010
      solai(n3) = soli(n3)                                         ! 28 Jun 2010
      drdt20save=drdt20
!      drdt20=23.d0
      testa(n3) = dd(0) + cmplx(sola(n3),solai(n3))*(
     &        dd(1) + cmplx(sola(n3),solai(n3))*(
     &        dd(2) + cmplx(sola(n3),solai(n3))*(
     &        dd(3) + cmplx(sola(n3),solai(n3))*dd(4))))
      if(drdt20.eq.23.d0) print *, 'root 3 testa(n3) = ', testa(n3)

      xxx = cmplx(sola(n3),solai(n3))

      if(isnan(abs(xxx))) print *, 'xxx = ', xxx, 'dd = ', dd

      testa(n3) = p4(dd,xxx,relerra(n3))

      if(drdt20.eq.23.d0) print *, 'root 3 testa(n3) = ', testa(n3)
      kay2 = (sola(n3)**2-solai(n3)**2)/cs
      kay2i=2.d0*sola(n3)*solai(n3)/cs
      error=errorcalc(ks,kay2)                                     ! 11 Mar 2011
      if(drdt20.eq.23.d0) print *, 'root 3 error = ', error
      if(drdt20.eq.23.d0) print *, 'root 3 kay2 = ', kay2
     &                       , 'sola(n3) = ', sola(n3), solai(n3)
      if(drdt20.eq.23.d0) print *, 'root 3 kay2i = ', kay2i
      kaysq(n3)=kay2                                               ! 11 Jun 2010
      kaysqa(n3)=kay2                                              ! 14 Jun 2010
      test(n3)=testa(n3)                                           ! 28 Jun 2010
!      drdt20=drdt20save
        
C ----------------------------------------------------------------------!-------
      drdt20save=drdt20
!      drdt20=23.d0
      if(drdt20.eq.23.d0) print *, 'dd(0) = ', dd(0) 
      if(drdt20.eq.23.d0) print *, 'dd(2) = ', dd(2) 

      if(dd(2).eq.0.d0) then

         solb(n3)=sola(n3)
         solbi(n3)=solai(n3)
         testb(n3)=testa(n3)
         kaysqb(n3) = kaysqa(n3) 
         relerrb(n3) = relerra(n3)                                    
         relerr(n3) = relerra(n3)                                    

      else

      kay2 = -dd(0)/(cs*dd(2))
      kay2i=0.d0

      if(kay2.lt.0.d0) then
!      drdt20=23.d0
         solb(n3)=0.d0
         solbi(n3)=c*dsqrt(-kay2)
      testb(n3) = dd(0) + cmplx(solb(n3),solbi(n3))*(
     &        dd(1) + cmplx(solb(n3),solbi(n3))*(
     &        dd(2) + cmplx(solb(n3),solbi(n3))*(
     &        dd(3) + cmplx(solb(n3),solbi(n3))*dd(4))))
      if(isnan(abs(testb(n3)))) print *, 'xxx = ', xxx

      if(drdt20.eq.23.d0) print *, 'testb(n3) = ', testb(n3)

      xxx = cmplx(solb(n3),solbi(n3))

      if(isnan(abs(xxx))) print *, 'xxx1 = ', xxx, 'dd = ', dd

      testb(n3) = p4(dd,xxx,relerrb(n3))

      if(isnan(abs(testb(n3)))) print *, 'xxx2 = ', xxx

      if(drdt20.eq.23.d0) print *, 'testb(n3) = ', testb(n3)
      error=errorcalc(ks,kay2)                                     ! 11 Mar 2011
      if(drdt20.eq.23.d0) print *, 'error = ', error
      if(drdt20.eq.23.d0) print *, 'kay2 = ', kay2
     &                       , 'solb(n3) = ', solb(n3), solbi(n3)
      drdt20save=drdt20

      else

      kay=dsqrt(kay2)
      solb(n3) = c*kay
      solbi(n3) = 0.d0
      testb(n3) = dd(0) + cmplx(solb(n3),solbi(n3))*(
     &        dd(1) + cmplx(solb(n3),solbi(n3))*(
     &        dd(2) + cmplx(solb(n3),solbi(n3))*(
     &        dd(3) + cmplx(solb(n3),solbi(n3))*dd(4))))
      if(isnan(abs(testb(n3)))) print *, 'xxx3 = ', xxx

      if(drdt20.eq.23.d0) print *, 'testb(n3) = ', testb(n3)

      xxx = cmplx(solb(n3),solbi(n3))

      if(isnan(abs(xxx))) print *, 'xxx4 = ', xxx, 'dd = ', dd

      testb(n3) = p4(dd,xxx,relerrb(n3))

      if(isnan(abs(testb(n3)))) print *, 'xxx5 = ', xxx

      if(drdt20.eq.23.d0) print *, 'testb(n3) = ', testb(n3)
      error=errorcalc(ks,kay2)                                     ! 11 Mar 2011
      if(drdt20.eq.23.d0) print *, 'error = ', error
      if(drdt20.eq.23.d0) print *, 'kay2 = ', kay2
     &                       , 'solb(n3) = ', solb(n3), solbi(n3)

      sgn = 1.d0
      do loop=1,6
      call iterate (dd,n3,sgn,solb,solbi,testb,relerrb,error
     & ,drdt20,loop,c,cs,ks,kay2,kay2i)                            ! 20 Jul 2010
      enddo

!      drdt20=drdt20save
      endif

      kaysqb(n3) = kay2                                            ! 11 Jun 2010

C ----------------------------------------------------------------------!-------
      relerr(n3) = relerra(n3)                                     ! 28 Jun 2010
      if(.not.(isnan(abs(testb(n3)))).and..not.(isnan(kaysqb(n3)))
     & .and.abs(testb(n3)).lt.abs(testa(n3)) 
     &   ) then
         kaysq(n3) = kaysqb(n3)                                    ! 24 Jun 2010
         test(n3)=testb(n3)                                        ! 28 Jun 2010
         relerr(n3) = relerrb(n3)                                  ! 28 Jun 2010
      endif

      endif
      kay2=kaysq(n3)

      return
C ----------------------------------------------------------------------!-------
C ----------------------------------------------------------------------!-------
      end subroutine root3
C ----------------------------------------------------------------------!-------
C ----------------------------------------------------------------------!-------
      subroutine iterate (dd,n3,sgn,solb,solbi,testb,relerrb,error
     & ,drdt20,loop,c,cs,ks,kay2,kay2i)                            ! 20 Jul 2010
C ----------------------------------------------------------------------!-------
      IMPLICIT double precision (a-h,o-z)
      dimension dd(0:4),solb(4),solbi(4)
      double precision ks,KAY2,KAY2I,kaysqb(4)              
      complex*16 testb(4),p4,relerrb(4),xxx,xxx2,addon

C ----------------------------------------------------------------------!-------
      xxx = cmplx(solb(n3),solbi(n3))
      xxx2=xxx*xxx

      addon=dd(1)*xxx+dd(3)*xxx2*xxx+dd(4)*xxx2*xxx2 
      if(dd(2).eq.0.d0.or.abs(addon).ge.abs(dd(0))) return

      xxx2=-(dd(0)+addon)/dd(2)
      xxx = sgn*sqrt(xxx2)

      if(isnan(abs(xxx))) print *, loop,'xxx = ', xxx, 'dd = ', dd

      testb(n3) = p4(dd,xxx,relerrb(n3))

      if(isnan(abs(testb(n3)))) print *, loop,'xxx = ', xxx

      solb(n3)=real(xxx)
      solbi(n3)=aimag(xxx)

      kay2=real(xxx2)/cs
      kay2i=aimag(xxx2)/cs

      if(drdt20.eq.23.d0) print *, 'testb(n3) = ', testb(n3)
      error=errorcalc(ks,kay2)                                     ! 11 Mar 2011
      if(drdt20.eq.23.d0) print *, 'error = ', error
      if(drdt20.eq.23.d0) print *, 'kay2 = ', kay2
     &                       , 'solb(n3) = ', solb(n3), solbi(n3)

      return
C ----------------------------------------------------------------------!-------
      end subroutine iterate                                       ! 20 Jul 2010
C ----------------------------------------------------------------------!-------
C ----------------------------------------------------------------------!-------
      SUBROUTINE root4(dd,sol,soli,Nsol
     & ,sola,solai,solb,solbi
     & ,test,testa,testb                                          
     & ,relerr,relerra,relerrb                                     ! 28 Jun 2010
     & ,nbig,nsmall,n3,n4                                         
     & ,height,ks,cs,c,error,errorsave
     & ,kaysq,kaysqa,kaysqb
     & ,kay2,drdt20)                                              
C ----------------------------------------------------------------------!-------
      IMPLICIT double precision (a-h,o-z)
      dimension dd(0:4),sol(4),soli(4)
     & ,solb(4),solbi(4),sola(4),solai(4)                          ! 15 Jun 2010
      complex*16 test(4),testa(4),testb(4)                         ! 24 Jun 2010
      REAL(KIND=8) :: KAY2,KAY2I,kay 
     & , kaysq(4) , kaysqb(4) , kaysqa(4)                          ! 14 Jun 2010
      REAL(KIND=8) :: ks,kv,kperp2
      dimension nn3(4,4), nn4(4,4)
!      double precision D1mach
      double precision errorcalc                                   ! 11 Mar 2011
      external errorcalc                                           ! 11 Mar 2011
      complex*16 p4,relerr(4),relerra(4),relerrb(4),xxx

C ----------------------------------------------------------------------!-------
      IF(n4.eq.0) return                                           ! 29 May 2013
      sola(n4) = sol(n4)                                           ! 28 Jun 2010
      solai(n4) = soli(n4)                                         ! 28 Jun 2010
      testa(n4) = dd(0) + cmplx(sola(n4),solai(n4))*(
     &        dd(1) + cmplx(sola(n4),solai(n4))*(
     &        dd(2) + cmplx(sola(n4),solai(n4))*(
     &        dd(3) + cmplx(sola(n4),solai(n4))*dd(4))))

      xxx = cmplx(sola(n4),solai(n4))
      testa(n4) = p4(dd,xxx,relerra(n4))

      if(drdt20.eq.23.d0) print *, 'root 4 testa(n4) = ', testa(n4)
      kay2 = (sola(n4)**2-solai(n4)**2)/cs
      kay2i=2.d0*sola(n4)*solai(n4)/cs
      error=errorcalc(ks,kay2)                                     ! 11 Mar 2011
      if(drdt20.eq.23.d0) print *, 'root 4 error = ', error
      if(drdt20.eq.23.d0) print *, 'root 4 kay2 = ', kay2
     &                       , 'sola(n4) = ', sola(n4), solai(n4)
      if(drdt20.eq.23.d0) print *, 'root 4 kay2i = ', kay2i
      kaysq(n4)=kay2                                               ! 11 Jun 2010
      kaysqa(n4)=kay2                                              ! 14 Jun 2010
      test(n4)=testa(n4)                                           ! 28 Jun 2010

C ----------------------------------------------------------------------!-------
      kay2 = -dd(0)/(cs*dd(2))
      kay2i=0.d0

      if(kay2.lt.0.d0) then
         solb(n4)=0.d0
         solbi(n4)=-c*dsqrt(-kay2)
      testb(n4) = dd(0) + cmplx(solb(n4),solbi(n4))*(
     &        dd(1) + cmplx(solb(n4),solbi(n4))*(
     &        dd(2) + cmplx(solb(n4),solbi(n4))*(
     &        dd(4) + cmplx(solb(n4),solbi(n4))*dd(4))))

      xxx = cmplx(solb(n4),solbi(n4))
      testb(n4) = p4(dd,xxx,relerrb(n4))

      if(drdt20.eq.23.d0) print *, 'testb(n4) = ', testb(n4)
      error=errorcalc(ks,kay2)                                     ! 11 Mar 2011
      if(drdt20.eq.23.d0) print *, 'error = ', error
      if(drdt20.eq.23.d0) print *, 'kay2 = ', kay2
     &                       , 'solb(n4) = ', solb(n4), solbi(n4)

      else

      kay=-dsqrt(kay2)
      solb(n4) = c*kay
      solbi(n4) = 0.d0
      testb(n4) = dd(0) + cmplx(solb(n4),solbi(n4))*(
     &        dd(1) + cmplx(solb(n4),solbi(n4))*(
     &        dd(2) + cmplx(solb(n4),solbi(n4))*(
     &        dd(3) + cmplx(solb(n4),solbi(n4))*dd(4))))

      xxx = cmplx(solb(n4),solbi(n4))
      testb(n4) = p4(dd,xxx,relerrb(n4))

      if(drdt20.eq.23.d0) print *, 'testb(n4) = ', testb(n4)
      error=errorcalc(ks,kay2)                                     ! 11 Mar 2011
      if(drdt20.eq.23.d0) print *, 'error = ', error
      if(drdt20.eq.23.d0) print *, 'kay2 = ', kay2
     &                       , 'solb(n4) = ', solb(n4), solbi(n4)

      sgn = -1.d0
      do loop=1,6
      call iterate (dd,n4,sgn,solb,solbi,testb,relerrb,error
     & ,drdt20,loop,c,cs,ks,kay2,kay2i)                            ! 20 Jul 2010
      enddo

      endif

      kaysqb(n4) = kay2                                            ! 11 Jun 2010

C ----------------------------------------------------------------------!-------
      relerr(n4) = relerra(n4)                                     ! 28 Jun 2010
      if(.not.(isnan(abs(testb(n4)))).and..not.(isnan(kaysqb(n4)))
     & .and.abs(testb(n4)).lt.abs(testa(n4)) 
     &   ) then
         kaysq(n4) = kaysqb(n4)                                    ! 24 Jun 2010
         test(n4)=testb(n4)                                        ! 28 Jun 2010
         relerr(n4) = relerrb(n4)                                  ! 28 Jun 2010
      endif

      kay2=kaysq(n4)
      return

C ----------------------------------------------------------------------!-------
C ----------------------------------------------------------------------!-------
      end subroutine root4
C ----------------------------------------------------------------------!-------
C ----------------------------------------------------------------------!-------
      SUBROUTINE smallroot(dd,sol,soli,Nsol
     & ,sola,solai,solb,solbi
     & ,test,testa,testb                                          
     & ,relerr,relerra,relerrb                                     ! 28 Jun 2010
     & ,nbig,nsmall,n3,n4                                         
     & ,height,ks,cs,c,error,errorsave
     & ,kaysq,kaysqa,kaysqb
     & ,kay2,drdt20)                                              
C ----------------------------------------------------------------------!-------
      IMPLICIT double precision (a-h,o-z)
      dimension dd(0:4),sol(4),soli(4)
     & ,solb(4),solbi(4),sola(4),solai(4)                          ! 15 Jun 2010
      complex*16 test(4),testa(4),testb(4)                         ! 24 Jun 2010
      REAL(KIND=8) :: KAY2,KAY2I,kay 
     & , kaysq(4) , kaysqb(4) , kaysqa(4)                          ! 14 Jun 2010
      REAL(KIND=8) :: ks,kv,kperp2
      dimension nn3(4,4), nn4(4,4)
!      double precision D1mach
      double precision errorcalc                                   ! 11 Mar 2011
      external errorcalc                                           ! 11 Mar 2011
      complex*16 p4,relerr(4),relerra(4),relerrb(4),xxx

C ----------------------------------------------------------------------!-------
      IF(nsmall.eq.0) return                                       ! 29 May 2013
      sola(nsmall) = sol(nsmall)                                   ! 28 Jun 2010
      solai(nsmall) = soli(nsmall)                                 ! 28 Jun 2010
      start = -dd(3)/(2.d0*dd(4))                                  ! 18 Jun 2010
      eps = (-dd(2)
     &      )/dd(4)                                                ! 18 Jun 2010
      rad = start*start + eps                                      ! 22 Jun 2010
      if(rad.lt.0.d0) then
         solb(nsmall) = start                    
         solbi(nsmall) = -sign(1.d0,start)*dsqrt(-rad)

      testb(nsmall) = dd(0) + cmplx(solb(nsmall),solbi(nsmall))*(
     &        dd(1) + cmplx(solb(nsmall),solbi(nsmall))*(
     &        dd(2) + cmplx(solb(nsmall),solbi(nsmall))*(
     &        dd(3) + cmplx(solb(nsmall),solbi(nsmall))*dd(4))))

      xxx = cmplx(solb(nsmall),solbi(nsmall))
      testb(nsmall) = p4(dd,xxx,relerrb(nsmall))

      if(drdt20.eq.23.d0) print *, 'root 2b test = ', testb(nsmall)
      kay2 = (solb(nsmall)**2-solbi(nsmall)**2)/cs
      kay2i=2.d0*solb(nsmall)*solbi(nsmall)/cs
      error=errorcalc(ks,kay2)                                     ! 11 Mar 2011
      if(drdt20.eq.23.d0) print *, 'root 2b error = ', error
      if(drdt20.eq.23.d0) print *, 'root 2b kay2 = ', kay2
     &        , 'solb(nsmall) = ', solb(nsmall), solbi(nsmall)
      if(drdt20.eq.23.d0) print *, 'root 2b kay2i = ', kay2i
        
      else
   
      solb(nsmall)=-sign(1.d0,start)*eps/(dabs(start)+dsqrt(rad))  ! 23 Jun 2010
      solbi(nsmall) = 0.d0
!      if(dabs(solb(nsmall))*eps**2/start**4.lt.8*dabs(start)*D1mach(3))
!     & solb(nsmall) = eps/(2.d0*dabs(start))                      ! 22 Jun 2010

      testb(nsmall) = dd(0) + cmplx(solb(nsmall),solbi(nsmall))*(
     &        dd(1) + cmplx(solb(nsmall),solbi(nsmall))*(
     &        dd(2) + cmplx(solb(nsmall),solbi(nsmall))*(
     &        dd(3) + cmplx(solb(nsmall),solbi(nsmall))*dd(n4))))

      xxx = cmplx(solb(nsmall),solbi(nsmall))
      testb(nsmall) = p4(dd,xxx,relerrb(nsmall))

      if(drdt20.eq.23.d0) print *, 'root 2b test = ', testb(nsmall)
      kay2 = (solb(nsmall)**2-solbi(nsmall)**2)/cs
      kay2i=2.d0*solb(nsmall)*solbi(nsmall)/cs
      error=errorcalc(ks,kay2)                                     ! 11 Mar 2011
      if(drdt20.eq.23.d0) print *, 'root 2b error = ', error
      if(drdt20.eq.23.d0) print *, 'root 2b kay2 = ', kay2
     &        , 'solb(nsmall) = ', solb(nsmall), solbi(nsmall)
      if(drdt20.eq.23.d0) print *, 'root 2b kay2i = ', kay2i
        
      eps = (-dd(2)
     &      +(dd(1)+dd(0)/solb(nsmall))/solb(nsmall)
     &      )/dd(4)                                                ! 18 Jun 2010
      rad = start*start + eps                                      ! 22 Jun 2010

      solb(nsmall)=-sign(1.d0,start)*eps/(dabs(start)+dsqrt(rad))  ! 23 Jun 2010
      solbi(nsmall) = 0.d0
!      if(dabs(solb(nsmall))*eps**2/start**4.lt.8*dabs(start)*D1mach(3))
!     & solb(nsmall) = eps/(2.d0*dabs(start))                      ! 22 Jun 2010

      testb(nsmall) = dd(0) + cmplx(solb(nsmall),solbi(nsmall))*(
     &        dd(1) + cmplx(solb(nsmall),solbi(nsmall))*(
     &        dd(2) + cmplx(solb(nsmall),solbi(nsmall))*(
     &        dd(3) + cmplx(solb(nsmall),solbi(nsmall))*dd(4))))

      xxx = cmplx(solb(nsmall),solbi(nsmall))
      testb(nsmall) = p4(dd,xxx,relerrb(nsmall))

      if(drdt20.eq.23.d0) print *, 'root 2b test = ', testb(nsmall)
      kay2 = (solb(nsmall)**2-solbi(nsmall)**2)/cs
      kay2i=2.d0*solb(nsmall)*solbi(nsmall)/cs
      error=errorcalc(ks,kay2)                                     ! 11 Mar 2011
      if(drdt20.eq.23.d0) print *, 'root 2b error = ', error
      if(drdt20.eq.23.d0) print *, 'root 2b kay2 = ', kay2
     &         , 'solb(nsmall) = ', solb(nsmall), solbi(nsmall)
      if(drdt20.eq.23.d0) print *, 'root 2b kay2i = ', kay2i
        
      eps = (-dd(2)
     &      +(dd(1)+dd(0)/solb(nsmall))/solb(nsmall)
     &      )/dd(4)                                                ! 18 Jun 2010
      rad = start*start + eps                                      ! 22 Jun 2010

      solb(nsmall)=-sign(1.d0,start)*eps/(dabs(start)+dsqrt(rad))  ! 23 Jun 2010
      solbi(nsmall) = 0.d0
!      if(dabs(solb(nsmall))*eps**2/start**4.lt.8*dabs(start)*D1mach(3))
!     & solb(nsmall) = eps/(2.d0*dabs(start))                      ! 22 Jun 2010

      testb(nsmall) = dd(0) + cmplx(solb(nsmall),solbi(nsmall))*(
     &        dd(1) + cmplx(solb(nsmall),solbi(nsmall))*(
     &        dd(2) + cmplx(solb(nsmall),solbi(nsmall))*(
     &        dd(3) + cmplx(solb(nsmall),solbi(nsmall))*dd(4))))

      xxx = cmplx(solb(nsmall),solbi(nsmall))
      testb(nsmall) = p4(dd,xxx,relerrb(nsmall))

      if(drdt20.eq.23.d0) print *, 'root 2b test = ', testb(nsmall)
      kay2 = (solb(nsmall)**2-solbi(nsmall)**2)/cs
      kay2i=2.d0*solb(nsmall)*solbi(nsmall)/cs
      error=errorcalc(ks,kay2)                                     ! 11 Mar 2011
      if(drdt20.eq.23.d0) print *, 'root 2b error = ', error
      if(drdt20.eq.23.d0) print *, 'root 2b kay2 = ', kay2
     &        , 'solb(nsmall) = ', solb(nsmall), solbi(nsmall)
      if(drdt20.eq.23.d0) print *, 'root 2b kay2i = ', kay2i
        
      endif                                                        ! 18 Jun 2010

      kaysqb(nsmall) = kay2                                        ! 11 Jun 2010

C ----------------------------------------------------------------------!-------
      testa(nsmall) = dd(0) + cmplx(sola(nsmall),solai(nsmall))*(
     &        dd(1) + cmplx(sola(nsmall),solai(nsmall))*(
     &        dd(2) + cmplx(sola(nsmall),solai(nsmall))*(
     &        dd(3) + cmplx(sola(nsmall),solai(nsmall))*dd(4))))

      xxx = cmplx(sola(nsmall),solai(nsmall))
      testa(nsmall) = p4(dd,xxx,relerra(nsmall))

      if(drdt20.eq.23.d0) print *, 'root 2 test = ', testa(nsmall)
      kay2 = (sola(nsmall)**2-solai(nsmall)**2)/cs
      kay2i=2.d0*sola(nsmall)*solai(nsmall)/cs
      error=errorcalc(ks,kay2)                                     ! 11 Mar 2011
      if(drdt20.eq.23.d0) print *, 'root 2 error = ', error
      if(drdt20.eq.23.d0) print *, 'root 2 kay2 = ', kay2
     &     , 'sola(nsmall) = ', sola(nsmall), solai(nsmall)
      if(drdt20.eq.23.d0) print *, 'root 2 kay2i = ', kay2i
      kaysq(nsmall)=kay2                                           ! 11 Jun 2010
      kaysqa(nsmall)=kay2                                          ! 14 Jun 2010
      test(nsmall)=testa(nsmall)                                   ! 28 Jun 2010
        
C ----------------------------------------------------------------------!-------
      relerr(nsmall) = relerra(nsmall)                             ! 28 Jun 2010
      if(.not.(isnan(abs(testb(nsmall))))
     & .and..not.(isnan(kaysqb(nsmall)))
     & .and.abs(testb(nsmall)).lt.abs(testa(nsmall)) 
     &   ) then
         kaysq(nsmall) = kaysqb(nsmall)                            ! 24 Jun 2010
         test(nsmall)=testb(nsmall)                                ! 28 Jun 2010
         relerr(nsmall) = relerrb(nsmall)                          ! 28 Jun 2010
      endif

      kay2=kaysq(nsmall)
      return

C ----------------------------------------------------------------------!-------
C ----------------------------------------------------------------------!-------
      end subroutine smallroot
C ----------------------------------------------------------------------!-------
C ----------------------------------------------------------------------!-------
      complex*16 function p4(dd,x,rel)
C ----------------------------------------------------------------------!-------
!      IMPLICIT double precision (a-h,o-z)
!      dimension dd(0:4)
      double precision dd(0:4)
      complex*16 x,rel

      p4=dd(0) + x*(dd(1) + x*(dd(2) + x*(dd(3) + x*(dd(4)))))

      if(isnan(abs(p4))) print *, 'x = ', x, 'dd = ', dd

      rel=p4/max(abs(dd(0)),abs(dd(1)*x),abs(dd(2)*x**2)
     & ,abs(dd(3)*x**3),abs(dd(4)*x**4))

      if(isnan(abs(rel))) print *, 'x = ', x, 'dd = ', dd, 'p4 = ', p4

C ----------------------------------------------------------------------!-------
      end function p4
C ----------------------------------------------------------------------!-------
C ----------------------------------------------------------------------!-------
      double precision function errorcalc(ks,kay2)                 ! 11 Mar 2011
C ----------------------------------------------------------------------!-------
      double precision ks,kay2
      if(ks.eq.kay2) then
         errorcalc = 0.d0
      else
         errorcalc = (ks-kay2)/max(dabs(ks),dabs(kay2))
      endif
      return
C ----------------------------------------------------------------------!-------
      end function errorcalc
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       ! 28 Jul 2017
      term2%frph = factor606*8.d0*coriolis%gsr*coriolis%gsth       ! 28 Jul 2017
      term2%fthth = factor606*4.d0*coriolis%gsth*coriolis%gsth
      term2%fthph = factor606*8.d0*coriolis%gsth*coriolis%gsph     ! 28 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,ks,kv                            !  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

C     COMMON DECK "UU" INSERTED HERE                               !  6 Jul 2017
      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

      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)              

      double precision ow,ows,owi,ow2                              !  6 Jul 2017

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

c      ows = ow**2

C     WIND VELOCITY
      CALL WINDR                                                   !  6 Jul 2017
      KV=KR*VR+KTH*VTH+KPH*VPH                                     !  6 Jul 2017
      VLS=KV*KV/KS                                                 !  6 Jul 2017
! OWI = intrinsic frequency
      OWI=OW-KV                                                    !  6 Jul 2017
!      ows = ow**2
      ows = owi**2                                                 !  6 Jul 2017
      ow2 = ow**2
      ks = kr**2 + kth**2 + kph**2                                 !  7 Jul 2017

      CALL GRAVg(gravity,coriolis)  

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

c      term4%f0 = -factor608*4.d0*ows*Omegasquare/cs                !  4 Aug 2016
      term4%f0 = -factor608*4.d0*ow2*Omegasquare/cs                !  7 Jul 2017
      term4%f1=0.d0                                                !  7 Jul 2017
      if(ks.ne.0.d0) 
     &term4%f1=factor608*8.d0*ow*kv*Omegasquare/(cs*sqrt(ks))      !  7 Jul 2017                                                !  7 Jul 2017
      term4%f2=0.d0                                                !  7 Jul 2017
      if(ks.ne.0.d0) 
     &term4%f2=-factor608*4.d0*kv*kv*Omegasquare/(cs*ks)           !  7 Jul 2017
      term4%fr=factor608*8.d0*vr*ow*Omegasquare/cs                 ! 10 Jul 2017
      term4%fth=factor608*8.d0*vth*ow*Omegasquare/cs               ! 10 Jul 2017
      term4%fph=factor608*8.d0*vph*ow*Omegasquare/cs               ! 10 Jul 2017
      term4%frr=-factor608*4.d0*vr*vr*Omegasquare/cs               ! 10 Jul 2017
      term4%frth=-factor608*8.d0*vr*vth*Omegasquare/cs             ! 10 Jul 2017
      term4%frph=-factor608*8.d0*vr*vph*Omegasquare/cs             ! 10 Jul 2017
      term4%fthth=-factor608*4.d0*vth*vth*Omegasquare/cs           ! 10 Jul 2017
      term4%fthph=-factor608*8.d0*vth*vph*Omegasquare/cs           ! 10 Jul 2017
      term4%fphph=-factor608*4.d0*vph*vph*Omegasquare/cs           ! 10 Jul 2017

      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
     & + factor608*8.d0*OWI*Omegasquare/cs
     & *(kr*pvrt+kth*pvtht+kph*pvpht)                              ! 10 Jul 2017

      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
     & + factor608*8.d0*OWI*Omegasquare/cs
     & *(kr*pvrr+kth*pvthr+kph*pvphr)                              ! 10 Jul 2017

      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
     & + factor608*8.d0*OWI*Omegasquare/cs
     & *(kr*pvrth+kth*pvthth+kph*pvphth)                           ! 10 Jul 2017

      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
     & + factor608*8.d0*OWI*Omegasquare/cs
     & *(kr*pvrph+kth*pvthph+kph*pvphph)                           ! 10 Jul 2017

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

      term4%pom=-factor608*8.d0*OWI*Omegasquare/cs                 ! 10 Jul 2017
      term4%pkr=factor608*8.d0*OWI*vr*Omegasquare/cs               ! 10 Jul 2017
      term4%pkth=factor608*8.d0*OWI*vth*Omegasquare/cs             ! 10 Jul 2017
      term4%pkph=factor608*8.d0*OWI*vph*Omegasquare/cs             ! 10 Jul 2017

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

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


C ----------------------------------------------------------------------!-------
      SUBROUTINE coriolis5(term5)                                  ! 22 Jul 2016
C ----------------------------------------------------------------------!-------
! 27 Jul 2017 - Corrected error in Coriolis5
! 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,ks,kv                            !  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

C     COMMON DECK "UU" INSERTED HERE                               !  6 Jul 2017
      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

      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)   

      double precision ow,ows,owi,ow2                              !  6 Jul 2017

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

c      ows = ow**2

C     WIND VELOCITY
      CALL WINDR                                                   !  6 Jul 2017
      KV=KR*VR+KTH*VTH+KPH*VPH                                     !  6 Jul 2017
      VLS=KV*KV/KS                                                 !  6 Jul 2017
! OWI = intrinsic frequency
      OWI=OW-KV                                                    !  6 Jul 2017
!      ows = ow**2
      ows = owi**2                                                 !  6 Jul 2017
      ow2 = ow**2
      ks = kr**2 + kth**2 + kph**2                                 !  7 Jul 2017

      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%f0 = -factor609*4.d0*ow2*Omegasquare/cs                !  7 Jul 2017
      term5%f1=0.d0                                                !  7 Jul 2017
      if(ks.ne.0.d0) 
!     & term5%f1=factor609*8.d0*kv*Omegasquare/(cs*sqrt(ks))        !  7 Jul 2017
     & term5%f1=factor609*8.d0*ow*kv*Omegasquare/(cs*sqrt(ks))     ! 27 Jul 2017
      term5%f2=0.d0
      if(ks.ne.0.d0) 
     & term5%f2=-factor609*4.d0*kv*kv*Omegasquare/(cs*ks)          !  7 Jul 2017
      term5%fr=factor609*8.d0*vr*ow*Omegasquare/cs                 ! 10 Jul 2017
      term5%fth=factor609*8.d0*vth*ow*Omegasquare/cs               ! 10 Jul 2017
      term5%fph=factor609*8.d0*vph*ow*Omegasquare/cs               ! 10 Jul 2017
      term5%frr=-factor609*4.d0*vr*vr*Omegasquare/cs               ! 10 Jul 2017
      term5%frth=-factor609*8.d0*vr*vth*Omegasquare/cs             ! 10 Jul 2017
      term5%frph=-factor609*8.d0*vr*vph*Omegasquare/cs             ! 10 Jul 2017
      term5%fthth=-factor609*4.d0*vth*vth*Omegasquare/cs           ! 10 Jul 2017
      term5%fthph=-factor609*8.d0*vth*vph*Omegasquare/cs           ! 10 Jul 2017
      term5%fphph=-factor609*4.d0*vph*vph*Omegasquare/cs           ! 10 Jul 2017

      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
     & + factor609*8.d0*OWI*Omegasquare/cs
     & *(kr*pvrt+kth*pvtht+kph*pvpht)                              ! 10 Jul 2017

      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
     & + factor609*8.d0*OWI*Omegasquare/cs
     & *(kr*pvrr+kth*pvthr+kph*pvphr)                              ! 10 Jul 2017

      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
     & + factor609*8.d0*OWI*Omegasquare/cs
     & *(kr*pvrth+kth*pvthth+kph*pvphth)                           ! 10 Jul 2017

      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
     & + factor609*8.d0*OWI*Omegasquare/cs
     & *(kr*pvrph+kth*pvthph+kph*pvphph)                           ! 10 Jul 2017

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

      term5%pom=-factor609*8.d0*OWI*Omegasquare/cs                 ! 10 Jul 2017
      term5%pkr=factor609*8.d0*OWI*vr*Omegasquare/cs               ! 10 Jul 2017
      term5%pkth=factor609*8.d0*OWI*vth*Omegasquare/cs             ! 10 Jul 2017
      term5%pkph=factor609*8.d0*OWI*vph*Omegasquare/cs             ! 10 Jul 2017

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,kv                            !  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

C     COMMON DECK "UU" INSERTED HERE                               !  6 Jul 2017
      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

      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

      double precision ow,ows,owi,ow2                              !  6 Jul 2017

C     WIND VELOCITY
      CALL WINDR                                                   !  6 Jul 2017
      KV=KR*VR+KTH*VTH+KPH*VPH                                     !  6 Jul 2017
      VLS=KV*KV/KS                                                 !  6 Jul 2017
! OWI = intrinsic frequency
      OWI=OW-KV                                                    !  6 Jul 2017
!      ows = ow**2
      ows = owi**2                                                 !  6 Jul 2017
      ow2 = ow**2
      ks = kr**2 + kth**2 + kph**2                                 !  7 Jul 2017

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 

      OmegaCrossGammaSubR= (
     &  coriolis%gsth*Gamma%gsph
     & -coriolis%gsph*Gamma%gsth
     & )                                                           !  7 Jul 2017
      OmegaCrossGammaSubTh= (
     &  coriolis%gsph*Gamma%gsr
     & -coriolis%gsr*Gamma%gsph
     & )                                                           !  7 Jul 2017
      OmegaCrossGammaSubPh= (
     &  coriolis%gsr*Gamma%gsth
     & -coriolis%gsth*Gamma%gsr
     & )                                                           !  7 Jul 2017
      OmegaCrossGammaDotK=OmegaCrossGammaSubR*kr
     & + OmegaCrossGammaSubTh*kth + OmegaCrossGammaSubPh*kph       !  7 Jul 2017

      term6%f0=0.d0
      term6%f1=0.d0                                                !  7 Jul 2017
      if(ks.ne.0.d0) 
     & term6%f1=factor610*4.d0*ow*OmegaCrossGammaDotK/sqrt(ks)      !  7 Jul 2017
      term6%f2=0.d0
      if(ks.ne.0.d0) 
     & term6%f2=-factor610*4.d0*kv*OmegaCrossGammaDotK/ks           !  7 Jul 2017
c      term6%fr=factor610*4.d0*ow*(
c     &  coriolis%gsth*Gamma%gsph
c     & -coriolis%gsph*Gamma%gsth
c     & )
      term6%fr=factor610*4.d0*ow*OmegaCrossGammaSubR               !  7 Jul 2017
c      term6%fth=factor610*4.d0*ow*(
c     &  coriolis%gsph*Gamma%gsr
c     & -coriolis%gsr*Gamma%gsph
c     & )
      term6%fth=factor610*4.d0*ow*OmegaCrossGammaSubTh             !  7 Jul 2017
c      term6%fph=factor610*4.d0*ow*(
c     &  coriolis%gsr*Gamma%gsth
c     & -coriolis%gsth*Gamma%gsr
c     & )
      term6%fph=factor610*4.d0*ow*OmegaCrossGammaSubPh             !  7 Jul 2017
      term6%frr=-factor610*4.d0*vr*OmegaCrossGammaSubR             ! 10 Jul 2017
      term6%frth=-factor610*4.d0*vr*OmegaCrossGammaSubTh
     & -factor610*4.d0*vth*OmegaCrossGammaSubR                     ! 10 Jul 2017
      term6%frph=-factor610*4.d0*vr*OmegaCrossGammaSubPh
     & -factor610*4.d0*vph*OmegaCrossGammaSubR                     ! 10 Jul 2017
      term6%fthth=-factor610*4.d0*vth*OmegaCrossGammaSubTh         ! 10 Jul 2017
      term6%fthph=-factor610*4.d0*vth*OmegaCrossGammaSubPh
     & -factor610*4.d0*vph*OmegaCrossGammaSubTh                    ! 10 Jul 2017
      term6%fphph=-factor610*4.d0*vph*OmegaCrossGammaSubPh         ! 10 Jul 2017

c      term6%value = term6%fr*kr+term6%fth*kth+term6%fph*kph
      term6%value = factor610*4.d0*OWI*OmegaCrossGammaDotK         !  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))
     & - factor610*4.d0*OmegaCrossGammaDotK
     & *(kr*pvrt+kth*pvtht+kph*pvpht)                              ! 10 Jul 2017

      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))
     & - factor610*4.d0*OmegaCrossGammaDotK
     & *(kr*pvrr+kth*pvthr+kph*pvphr)                              ! 10 Jul 2017

      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))
     & - factor610*4.d0*OmegaCrossGammaDotK
     & *(kr*pvrth+kth*pvthth+kph*pvphth)                           ! 10 Jul 2017

      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))
     & - factor610*4.d0*OmegaCrossGammaDotK
     & *(kr*pvrph+kth*pvthph+kph*pvphph)                           ! 10 Jul 2017

c      term6%pom=term6%value/ow
      term6%pom=factor610*4.d0*OmegaCrossGammaDotK                 ! 11 Jul 2017
c      term6%pkr=term6%fr
      term6%pkr=factor610*4.d0*OWI*OmegaCrossGammaSubR
     & -factor610*4.d0*vr*OmegaCrossGammaDotK                      ! 11 Jul 2017
c      term6%pkth=term6%fth
      term6%pkth=factor610*4.d0*OWI*OmegaCrossGammaSubTh
     & -factor610*4.d0*vth*OmegaCrossGammaDotK                     ! 11 Jul 2017
c      term6%pkph=term6%fph
      term6%pkph=factor610*4.d0*OWI*OmegaCrossGammaSubPh
     & -factor610*4.d0*vph*OmegaCrossGammaDotK                     ! 11 Jul 2017

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

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

