C file ANCNL.FOR
C ----------------------------------------------------------------------!-------
!  5 Aug 2016 - Added common block /rink3/ to help with reflection of wave
! 13 Apr 2011 - Corrected calculation of kz2a for case where kr2, kth2, or kph2 
!               might be negative
! 13 Apr 2011 - Added common block /rink2/
!  9 Dec 2010 - Added calculation of Hrel
!  9 Dec 2010 - Added calculation of errRel
!  9 Dec 2010 - Added common block /rink/
C 13 Mar 2008 - Corrected length in common block /PCONST/
C 19 Jul 2006 - Added calls to igrav and idensity
C 19 Jul 2006 - Added dummy subroutines for brunt and cutoff
C 19 Jul 2006 - Added calls to setgrav and setdensity
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 ----------------------------------------------------------------------!-------
      SUBROUTINE ANCNL
C ----------------------------------------------------------------------!-------
C     DISPERSION RELATION FOR ACOUSTIC WAVES NO CURRENT, NO LOSSES

C     COMMON DECK "RKAM" INSERTED HERE 
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DOUBLE PRECISION KR,KTH,KPH
      COMMON//R,TH,PH,KR,KTH,KPH,RKVARS(14),TPULSE,CSTEP,DRDT(20)

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

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
      LOGICAL SPACE
      DOUBLE PRECISION LPOLAR,LPOLRI,KPHK,KPHKI,KAY2,KAY2I 
     & , kz2a,kz2ai,kzb,kzbi  !16 Mar 2006
      CHARACTER DISPM*6,modrin*64
      COMMON/RINPL/DISPM
      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)                   !  9 Dec 2010
     & ,kaysqi(4),kaysqbi(4),kaysqai(4)                           
     & ,sol(4),soli(4),solb(4),solbi(4),sola(4),solai(4)          
     & ,alpha,dd                                                  
     & ,relerr(4),relerra(4),relerrb(4)                           
     & ,errRell(4),errRela(4),errRelb(4)                          
     & ,n3,n4,nbig,nsmall                                         
     & ,Hrel,errRel                                               
      complex*16 relerr,relerra,relerrb                            !  9 Dec 2010
      complex*16 errRell,errRela,errRelb                          
      double precision kaysq,kaysqb,kaysqa                        
     & ,kaysqi,kaysqbi,kaysqai                                    
      dimension dd(0:4)
!     & ,sol(4),soli(4),solb(4),solbi(4),sola(4),solai(4)          

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

      COMMON/WW/dum(10),MAXW,W(1000)
      DOUBLE PRECISION MAXSTP,MAXERR,INTYP,LLAT,LLON
      DOUBLE PRECISION KS,MKS
      DOUBLE PRECISION KVECT(24)
      DOUBLE PRECISION VSET(20)
      EQUIVALENCE(KVECT,KAY2),(VSET,V),(ow,w(6))
      DOUBLE PRECISION MRDRIN(8)    ,TBPE( 3)
      CHARACTER DLSPM*5
      DATA DLSPM/'ANCNL'/
      DATA TBPE/3*1H0 / 
C 
C 
      ENTRY SETDSP
      DISPM=DLSPM 
      CALL RMOVE(TYPE,TBPE, 3)
	modrin = 'ACOUSTIC WAVE --- NO CURRENT --- NO LOSSES'
C
      CALL SETSPD 
      CALL SETRCV 
      CALL SETTOP 
      CALL SETSUR 
      CALL setgrav ! 19 Jul 2006
      CALL setdensity ! 19 Jul 2006
      RETURN      
C
      ENTRY IDISPER 
      CALL ISPEED   
      CALL IRECVR   
      CALL ITOPOG   
      CALL ISURFAC  
      call igrav ! 19 Jul 2006
      call idensity ! 19 Jul 2006
      RETURN   
c
      ENTRY DISPER
      ENTRY RINDEX
      SPACE=.FALSE. 
      KS=KR*KR+KTH*KTH+KPH*KPH
C     SOUND SPEED   
      CALL SPEED  
      OWS=OW*OW 
      KAY2=OWS/CS

      kz2a=(ows/cs)-(kth2+kph2)                      ! 16 Mar 2006 ! 13 Apr 2011
      kz2ai = 0.d0 ! 16 Mar 2006
      kzb = 0.d0 ! 16 Mar 2006
      kzbi = 0.d0 ! 16 Mar 2006

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

      H=OW*OW - CS*KS

      H1= OW*OW                                                    !  9 Dec 2010
      H2= CS*KS                                                    !  9 Dec 2010
      Hrel=H/max(abs(H1),abs(H2))                                  !  9 Dec 2010

      MKS=-KS
      PHT=MKS*PCST
      PHR=MKS*PCSR
      PHTH=MKS*PCSTH
      PHPH=MKS*PCSPH
C
      PHOW=2.d0*OW  
      CS2=-2.d0*CS  
      PHKR=CS2*KR   
      PHKTH=CS2*KTH 
      PHKPH=CS2*KPH 
      KPHK=CS2*KS   

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

      RETURN        
C ----------------------------------------------------------------------!-------
      END SUBROUTINE ANCNL          
C ----------------------------------------------------------------------!-------
      SUBROUTINE cutoff(kA2)
C ----------------------------------------------------------------------!-------

C Just a dummy routine to satisfy loader

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 

      RETURN 
C ----------------------------------------------------------------------!-------
      END SUBROUTINE cutoff
C ----------------------------------------------------------------------!-------
      SUBROUTINE Brunt(Ns)
C ----------------------------------------------------------------------!-------

C Just a dummy routine to satisfy loader

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

      TYPE (value_and_derivs) :: Ns

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

                                                                                