C-----------------------------------------------------------------------------
C file sbloss.for
C-----------------------------------------------------------------------------
C 19 Mar 2008 - Added the calculation of Tr
C 13 Mar 2008 - Corrected length in common block /RIN/
C     15 Mar 2007 - Corrected some errors
C     13 Mar 2007 - Finished programming loss model
C      9 Mar 2007 - Added molecular weights to /pconst/
C      5 Mar 2007 - Added molecular fractions to common block /MM/
C-----------------------------------------------------------------------------
C     BACKGROUND ABSORPTION FORMULA
C References:
C Sutherland, Louis C. and Henry E. Bass, JASA 115, 1012-1032, 2004.
C Sutherland, Louis C. and Henry E. Bass, JASA 120, 2985, 2006.
C Bass, Henry E. and Claus H. Hetzer, Inframatics, The newsletter of 
C   subaudible sound, Number 15 September 2006, pp. 1-5.
C-----------------------------------------------------------------------------
C     24 Jul 2006
C     18 Oct 2003
C     17 December 2002, version to use with pcharpo
C     6 May 2003 - error corrected in calculation of aph
C     13 June 2003 - corrected sqrt to dsqrt in calculation of aph
C-----------------------------------------------------------------------------
      SUBROUTINE SBLOSS
C-----------------------------------------------------------------------------
      implicit double precision (a-h,o-z)
C
      double precision nu,n,N2rat, N2ratSQ
      complex k1, k1SQ, k2, i

C     COMMON DECK "RKAM" INSERTED HERE                                  
      DOUBLE PRECISION KR,KTH,KPH,KV
      COMMON//R,TH,PH,KR,KTH,KPH,RKVARS(14),TPULSE,CSTEP,DRDT(20)       
C
C     COMMON DECK "WW" INSERTED HERE
      PARAMETER (NWARSZ=1000)
      COMMON/WW/ dum(10),MAXW,W(NWARSZ)
      EQUIVALENCE (EARTHR,W(1)),(OW,W(6))
C
      EQUIVALENCE (W(503),mu0),(W(504),S)
     1 ,(W(505),ZN2inf),(W(506),T3N2),(W(507),ZO2inf),(W(508),T3O2)
     2 ,(W(509),thO2),(W(510),thN2),(W(511),a1st),(W(512),a2st)
     3 ,(W(513),bst),(W(514),cst),(W(515),dst),(W(516),est)
     4 ,(W(517),gst),(W(518),c1),(W(519),c2),(W(520),c3)
     5 ,(W(521),c4),(W(522),c5),(W(523),c6),(W(524),c7)
      double precision mu0
C
      EQUIVALENCE (W(553),p0)

C     COMMON DECK "RINREAL" INSERTED HERE
      LOGICAL SPACE
      DOUBLE PRECISION LPOLAR,LPOLRI,KPHK,KPHKI,KAY2,KAY2I
     & , kz2a,kz2ai,kzb,kzbi                         ! 16 Mar 2006 ! 13 Mar 2008
      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 ! 13 Mar 2008
     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 
C     COMMON DECK "LL" INSERTED HERE                           
      DOUBLE PRECISION MODL
      COMMON/LL/MODL(4),APH,APHPT,APHPR,APHPTH,APHPPH 
C
C     COMMON DECK "AA" INSERTED HERE
      DOUBLE PRECISION MODA
      DOUBLE PRECISION MU,MUPT,MUPR,MUPTH,MUPPH
      DOUBLE PRECISION KAP,KAPPT,KAPPR,KAPPTH,KAPPPH
      COMMON/AA/MODA(4),MU,MUPT,MUPR,MUPTH,MUPPH
      COMMON/AA/KAP,KAPPT,KAPPR,KAPPTH,KAPPPH
C
C     COMMON DECK "CONST" INSERTED HERE
      COMMON/PCONST/CREF,RGAS,GAMMA
     1 ,MO2,MN2,MO,MN,MH2O,MO3                        ! 9 Mar 2007 ! 13 Mar 2008
      COMMON/MCONST/PI,PIT2,PID2,DEGS,RAD,ALN10
C
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                                  
C
C     COMMON DECK "CC" INSERTED HERE
      DOUBLE PRECISION MODC
      COMMON/CC/MODC(4),CsSQ,PCST,PCSR,PCSTH,PCSPH

C     COMMON DECK "PP" INSERTED HERE
      DOUBLE PRECISION MODP
      COMMON/PP/MODP(4),P,PPT,PPR,PPTH,PPPH
C
C     COMMON DECK "TT" INSERTED HERE
      DOUBLE PRECISION MODT
      COMMON/TT/MODT(4), T,PTT,PTR,PTTH,PTPH
C
C     COMMON DECK "MM" INSERTED HERE
      DOUBLE PRECISION M,MODM
     1 ,MO2,MN2,MO,MN,MH2O,MO3 ! 5 Mar 2007
      COMMON/MM/MODM(4),M,PMT,PMR,PMTH,PMPH
     1 ,XO2,XN2,XO,XN,XH2O,XO3 ! 5 Mar 2007
C
C     COMMON DECK "CB17" INSERTED HERE                
      INTEGER VMX,VNTBL,VITBL,VFRMTBL,IDSV(10)        
      COMMON/B17/VMX,VNTBL(10),VITBL(10),VFRMTBL(10),VGP(53)
      EQUIVALENCE (VGP,IDSV),(ANV,VGP(11))                  
      INTEGER VPX    ,VQTBL(10),VLTBL(10),VIRMTBL(10)       
      DATA VPX/1/         
      DATA VQTBL/1,11,8*0/
      DATA VLTBL/1,9*0/   
      DATA VIRMTBL/1,9*0/ 

      data i /(0.0,1.0)/
      data T0,F /293.15d0,60000.d0/
C                         
      ENTRY SETABS        
      VMX=VPX                       
      CALL IMOVE(VNTBL,VQTBL,10)    
      CALL IMOVE(VITBL,VLTBL,10)    
      CALL IMOVE(VFRMTBL,VIRMTBL,10)
        call setspd ! sound speed
	call SETPRES ! pressure
      CALL SETPAB ! absorption perturbation                  
      RETURN                        
C                                   
      ENTRY IABSRP                  
      IF(W(500) .NE. 3.) THEN
           PRINT*,' Model check number mismatch: '
           STOP ' W(500) should be 3 if you are using model LOSSBASS.'
      endif
C
      MODL(1)=6HSBLOSS                           
      MODL(2)= w(502)
        call ispeed ! sound speed
	call IPRES
      CALL IPABSRP                                 
      RETURN                                       
C                                                  
      ENTRY ABSRP                                  
C     WIND VELOCITY                                                     
      CALL WINDR                                                        
      KV=KR*VR+KTH*VTH+KPH*VPH                                          
      OWI=OW-KV                                                         
C     sound speed
      call speed ! sound speed
      cs = dsqrt(csSQ)
C     PRESSURE
      CALL PRES
C-----------------------------------------------------------------------------

      mu = mu0*(T/T0)**1.5*(T0+S)/(T+S)
c      print *, 'mu = ', mu
      nu = 4.d0 * OWI*mu/(3.d0*p)      
c      print *, 'nu = ', nu
      k1SQ = -1.d0/(1.d0+i*nu)
      k1 = csqrt(k1SQ)
      sigma = 5.d0/dsqrt(21.d0)
      sigmaSQ = sigma**2
      ZrotN2 = ZN2inf*Dexp(-T3N2/T**(1.d0/3.d0))
      ZrotO2 = ZO2inf*Dexp(-T3O2/T**(1.d0/3.d0))
      Zrot = 1.d0/(XN2/ZrotN2 + XO2/ZrotO2)
      n = 0.8d0*dsqrt(3.d0/7.d0)*Zrot
c      print *, 'n = ', n
      x = 3.d0*n*nu/4.d0
c      print *, 'x = ', x
      xprime = c1*x
      xpSQ = xprime**2
      k2 = ((sigmaSQ-1.d0)*x+1.d0*i*sigma*(1.d0+xpSQ))/
     1 (2.d0*sigma*dsqrt((1.d0+xpSQ)*(1.d0+sigmaSQ*xpSQ)))
      thN2dT = thN2/T
      exN2 = dexp(-thN2dT)
      CpN2dR = thN2dT**2*exN2/(1.d0-exN2)**2
      thO2dT = thO2/T
      exO2 = dexp(-thO2dT)
      CpO2dR = thO2dT**2*exO2/(1.d0-exO2)**2
c      print *, cpo2dr
      AmaxN2dp = (XN2/2.d0)*(CpN2dR)/(3.5d0*(2.5d0+CpN2dR))
      AmaxO2dp = (XO2/2.d0)*(CpO2dR)/(3.5d0*(2.5d0+CpO2dR))
      Tr = (T/T0)**(-1.d0/3.d0)-1.d0                               ! 19 Mar 2008
      A1 = (XO2+XN2)*a1st*dexp(-c2*Tr)
      A2 = (XO+XN)*a2st
      B = bst*dexp(c3*Tr)
c        B = 0.d0
      C = cst*dexp(-c4*Tr)
      D = dst*dexp(c5*Tr)
      E = est*dexp(-c6*Tr)
c        F = 0.d0
      G = gst*dexp(-c7*Tr)
c        G = 0.d0
      hprime = 100.d0*(XH2O+XO3)
      omvibN2 = PIT2 *p/p0*mu0/mu*(E+F*XO3+G*XH2O)
      omvibO2 = PIT2 *p/p0*mu0/mu*(A1+A2+B*hprime*(C+hprime)*(D+hprime))
c      print *, omvibo2
      N2rat = OWI/omvibN2
      N2ratSQ = N2rat**2
      vibN2 = AmaxN2dp*N2rat/(1.d0+N2ratSQ)
c      vibN2 = 0.d0
      O2rat = OWI/omvibO2
      O2ratSQ = O2rat**2
      vibO2 = AmaxO2dp*O2rat/(1.d0+O2ratSQ)
c      vibO2 = 0.d0
c         k2=(0.d0,0.d0)
c         k2=0.d0
      APH = (OWI/Cs)*(REAL(-i*k1*k2) + vibN2 + vibO2)
c         db = 20.d0*aph/aln10
c         thelog=0.d0
c         if(db.gt.0.d0) thelog = dlog10(db)
c         height = R-earthr
c         print *, height, thelog
c      if(nu.gt.0.1d0) print *, k1SQ, k1, aph
C-----------------------------------------------------------------------------
      CALL PABSRP       
      END               
