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