C file awcwl.for
C ----------------------------------------------------------------------!-------
      SUBROUTINE AWCWL
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
! 19 Nov 2010 - Added errRell,errRelb,errRela to common block /rink/
! 13 Sep 2010 - Added calculation of errRel
! 13 Sep 2010 - Added common block /rink/
C 13 Mar 2008 - Corrected length in common block /PCONST/
C 27 Mar 2007 - Temporary call to densitya
C 23 Feb 2007 - Changed the definition of APH to be an attenuation coefficient
C               in amplitude instead of in power
C 26 Jul 2006 - Added CALL SETPRES and CALL IMOLWT and SETMW
C 25 Jul 2006 - Added CALL SETTMP and call itemp
C 19 Jul 2006 - Added calls to igrav and idensity
C 13 Jul 2006 - Added dummy subroutines for brunt and cutoff
C 13 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 ----------------------------------------------------------------------!-------
C     DISPERSION RELATION FOR ACOUSTIC WAVES WITH CURRENT, WITH 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)                  ! 13 Sep 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)                          ! 19 Nov 2010
     & ,n3,n4,nbig,nsmall                                         
     & ,Hrel,errRel                                               
      complex*16 relerr,relerra,relerrb                           
      complex*16 errRell,errRela,errRelb                          ! 19 Nov 2010
      double precision kaysq,kaysqb,kaysqa                        
     & ,kaysqi,kaysqbi,kaysqai                                    
      dimension dd(0: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 KS,KV
      DOUBLE PRECISION KVECT(24)
      EQUIVALENCE(KVECT,KAY2),(ow,w(6)) 
      DOUBLE PRECISION  TBPE( 3),modl
      CHARACTER DLSPM*5 
      COMMON/LL/MODL(4),APH,APHPT,APHPR,APHPTH,APHPPH 
      DATA DLSPM/'AWCWL'/
      DATA TBPE/3*1H3 / 
C
C
      ENTRY SETDSP  
	modrin = 'ACOUSTIC WAVE --- WITH CURRENT --- WITH LOSSES'
      DISPM=DLSPM 
      CALL RMOVE(TYPE,TBPE, 3) 
c      print*,'before setspd'
      CALL SETSPD
c      print*,'before setwnd'
      CALL SETWND
c      print*,'before setrcv'
      CALL SETRCV
c      print*,'before settop'
      CALL SETTOP
c      print*,'before setsur'
      CALL SETSUR
c      print*,'before setabs'
      CALL SETABS
c      print*,'before setmw'
      CALL SETMW ! 26 Jul 2006
c      print*,'before setpres'
      CALL SETPRES ! 26 Jul 2006
c      print*,'before settmp'
      CALL SETTMP ! 25 Jul 2006
c      print*,'before setgrav'
      CALL setgrav ! 13 Jul 2006
c      print*,'before setdensity'
      CALL setdensity ! 13 Jul 2006
      RETURN
C 
      ENTRY IDISPER
      CALL ISPEED
      CALL IWINDR
      CALL IRECVR
      CALL ITOPOG
      CALL ISURFAC
      CALL IABSRP
      CALL IPRES ! 26 Jul 2006
      CALL IMOLWT ! 26 Jul 2006
      CALL ITEMP ! 25 Jul 2006
      call igrav ! 19 Jul 2006
      call idensity ! 19 Jul 2006
c      print*,'After call to idensity'
      RETURN
C
      ENTRY DISPER
c      print*,'Beginning of disper'
      SPACE=.FALSE.
      KS=KR*KR+KTH*KTH+KPH*KPH
C     WIND VELOCITY 
      CALL WINDR
c      print *, ' vr = ', vr
c      print *, ' vth = ', vth
c      print *, ' vph = ', vph
c      print *, ' kr = ', kr
c      print *, ' kth = ', kth
c      print *, ' kph = ', kph
      KV=KR*VR+KTH*VTH+KPH*VPH
      VLS=KV*KV/KS
C     SOUND SPEED
      CALL SPEED
C     DENSITY                                           
      CALL DENSITYA                                                ! 27 Mar 2007
      OWS=OW*OW
c      print *, ' ows = ', ows
c      print *, ' cs = ', cs
c      print *, ' kv = ', kv
c      print *, ' ks = ', ks
      KAY2=OWS/(dSQRT(CS)+KV/dSQRT(KS))**2  
c      print *, ' kay2 = ', kay2

      kz2a=(ows*cs-cs*(kth2+kph2)*(cs-vr*vr))/(cs-vr*vr)!16 Mar 2006!13 Apr 2011
      kz2ai = 0.d0 ! 16 Mar 2006
      kzb = -vr*ow/(cs-vr*vr) ! 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

      OWI=OW-KV

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

      POWIT=-KR*PVRT - KTH*PVTHT - KPH*PVPHT 
      PHT=2.d0*OWI*POWIT - KS*PCST 
      POWIR=-KR*PVRR - KTH*PVTHR - KPH*PVPHR 
      PHR=2.d0*OWI*POWIR - KS*PCSR
      POWITH=-KR*PVRTH - KTH*PVTHTH - KPH*PVPHTH
      PHTH=2.d0*OWI*POWITH - KS*PCSTH
      POWIPH=-KR*PVRPH - KTH*PVTHPH - KPH*PVPHPH  
      PHPH=2.d0*OWI*POWIPH - KS*PCSPH
      PHOW=2.d0*OWI
      PHKR=-2.d0*(OWI*VR + CS*KR)   
      PHKTH=-2.d0*(OWI*VTH + CS*KTH)
      PHKPH=-2.d0*(OWI*VPH + CS*KPH)
      KPHK=-2.d0*(OWI*KV + CS*KS)   
      CALL ABSRP
      GMS=GAMMA-1.d0
C      KAY2I=-dSQRT(KAY2)*APH
      KAY2I=-2.d0*dSQRT(KAY2)*APH ! 23 Feb 2007 

      errRel=1.d10                                                 ! 13 Sep 2010
      IF (KAY2.NE.0.d0) errRel =(KR**2+KTH**2+KPH**2)/KAY2-1.d0    ! 13 Sep 2010
c      print *, ' errRel = ', errRel

      RETURN 
C ----------------------------------------------------------------------!-------
      END SUBROUTINE AWCWL
C ----------------------------------------------------------------------!-------
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 ----------------------------------------------------------------------!-------
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 ----------------------------------------------------------------------!-------

