C-----------------------------------------------------------------------------
      SUBROUTINE MUARDC   
C-----------------------------------------------------------------------------
C     ARDC BACKGROUND ABSORPTION FORMULA
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     Compile with options /c (compile only) and /4Yb (error tracing)
      implicit double precision (a-h,o-z)
C
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),BETAV),(W(504),SUTH),(W(505),PRNDTL) 
C
C     COMMON DECK "RINREAL" INSERTED HERE
      LOGICAL SPACE
      DOUBLE PRECISION LPOLAR,LPOLRI,KPHK,KPHKI,KAY2,KAY2I
      CHARACTER DISPM*6,modrin*64
      COMMON/RINPL/DISPM 
      COMMON /RIN/ MODRIN,RAYNAME(2,3),TYPE(3),SPACE
      COMMON/RIN/OMEGMIN,OMEGMAX,KAY2,KAY2I,
     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
      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 "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
      COMMON/MM/MODM(4),M,PMT,PMR,PMTH,PMPH
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/ 
C                         
      ENTRY SETABS        
      VMX=VPX                       
      CALL IMOVE(VNTBL,VQTBL,10)    
      CALL IMOVE(VITBL,VLTBL,10)    
      CALL IMOVE(VFRMTBL,VIRMTBL,10)
	call SETPRES
      CALL SETPAB                   
      RETURN                        
C                                   
      ENTRY IABSRP                  
      IF(W(500) .NE. 2.) THEN
           PRINT*,' Model check number mismatch: '
           STOP ' W(500) should be 2 if you are using model MUARDC.'
      endif
C
      MODL(1)=8HMUARDC                             
      MODL(2)= w(502)
	call IPRES
      CALL IPABSRP                                 
      RETURN                                       
C                                                  
      ENTRY ABSRP                                  
      MU=BETAV*T**1.5d0/(SUTH+T)
      KAP=GAMMA*RGAS*MU/((GAMMA-1.d0)*M*PRNDTL)
      GMS=GAMMA-1.d0                                    
C     WIND VELOCITY                                                     
      CALL WINDR                                                        
      KV=KR*VR+KTH*VTH+KPH*VPH                                          
      OWI=OW-KV                                                         
C     PRESSURE
      CALL PRES
C The first version below will duplicate the published error for the sample
C case in the HARPA report.
C      APH=(OW/(GAMMA*P))*(4./3.*MU+GMS*GMS*M/(GAMMA*RGAS)*KAP)
	APH=(OWI/(GAMMA*P))*(4./3.*MU+GMS*GMS*M/(GAMMA*RGAS)*KAP)
     1 *dsqrt(kay2)
      CALL PABSRP       
      END               
