! PROGRAM MassDistribution
! PURPOSE - Determine the fraction of the total atmosphere mass that lies
!   below a given altitude
! AUTHOR - Ralph L. Carmichael, Public Domain Aeronautical Software
! REVISION HISTORY
!   DATE  VERS PERSON  STATEMENT OF CHANGES
! 11Oct02  1.0   RLC   Initial release to the Web site

!+
MODULE Extra
! ---------------------------------------------------------------------------
! PURPOSE - Procedures used by main program
IMPLICIT NONE
!  INTEGER,PARAMETER,PRIVATE:: WP = SELECTED_REAL_KIND(10,50)
!----------------------------------------------------------------------------

CONTAINS


SUBROUTINE Atmosphere(alt, sigma, delta, theta)
! ---------------------------------------------------------------------------
! PURPOSE - Compute the properties of the 1976 standard atmosphere to 86 km.

  IMPLICIT NONE
  REAL,INTENT(IN)::  alt    ! geometric altitude, km.                        
  REAL,INTENT(OUT):: sigma  ! density/sea-level standard density              
  REAL,INTENT(OUT):: delta  ! pressure/sea-level standard pressure           
  REAL,INTENT(OUT):: theta  ! temperature/sea-level standard temperature

  REAL,PARAMETER:: REARTH = 6369.0                 ! radius of the Earth (km)
  REAL,PARAMETER:: GMR = 34.163195                             ! gas constant
  INTEGER,PARAMETER:: NTAB=8       ! number of entries in the defining tables

  INTEGER:: i,j,k                                                  ! counters
  REAL:: h                                       ! geopotential altitude (km)
  REAL:: tgrad, tbase      ! temperature gradient and base temp of this layer
  REAL:: tlocal                                           ! local temperature
  REAL:: deltah                             ! height above base of this layer

  REAL,DIMENSION(NTAB),PARAMETER:: htab= &
                          (/0.0, 11.0, 20.0, 32.0, 47.0, 51.0, 71.0, 84.852/)
  REAL,DIMENSION(NTAB),PARAMETER:: ttab= &
          (/288.15, 216.65, 216.65, 228.65, 270.65, 270.65, 214.65, 186.946/)
  REAL,DIMENSION(NTAB),PARAMETER:: ptab= &
               (/1.0, 2.233611E-1, 5.403295E-2, 8.5666784E-3, 1.0945601E-3, &
                                     6.6063531E-4, 3.9046834E-5, 3.68501E-6/)
  REAL,DIMENSION(NTAB),PARAMETER:: gtab= &
                                (/-6.5, 0.0, 1.0, 2.8, 0.0, -2.8, -2.0, 0.0/)
!----------------------------------------------------------------------------
  h=alt*REARTH/(alt+REARTH)      ! convert geometric to geopotential altitude

  i=1 
  j=NTAB                                       ! setting up for binary search
  DO
    k=(i+j)/2                                              ! integer division
    IF (h < htab(k)) THEN
      j=k
    ELSE
      i=k
    END IF   
    IF (j <= i+1) EXIT
  END DO

  tgrad=gtab(i)                                     ! i will be in 1...NTAB-1
  tbase=ttab(i)
  deltah=h-htab(i)
  tlocal=tbase+tgrad*deltah
  theta=tlocal/ttab(1)                                    ! temperature ratio

  IF (tgrad == 0.0) THEN                                     ! pressure ratio
    delta=ptab(i)*EXP(-GMR*deltah/tbase)
  ELSE
    delta=ptab(i)*(tbase/tlocal)**(GMR/tgrad)
  END IF

  sigma=delta/theta                                           ! density ratio
  RETURN
END Subroutine Atmosphere   ! -----------------------------------------------
  
!+
FUNCTION Integrand(z) RESULT(f)
! ---------------------------------------------------------------------------
! PURPOSE - Return the value of the integrand for mass of atmosphere.
! NOTE - Since we are using kilometers as the unit of length, RHOZERO is the
!  sea-level density in g/km**3, which is a rather large number.

  REAL,INTENT(IN):: z   ! altitude in kilometers
  REAL:: f
  
  REAL,PARAMETER:: PI=3.14159265, FOURPI=4.0*PI
  REAL,PARAMETER:: REARTH = 6369.0                 ! radius of the Earth (km)
  REAL,PARAMETER:: RHOZERO = 1.2250E12        ! density at sealevel, g/cu.km
  REAL:: sigma,delta,theta,z4
!----------------------------------------------------------------------------
  z4=z
  CALL Atmosphere(z4, sigma,delta,theta) 
  f=FOURPI*sigma*RHOZERO*(REARTH+z)**2
  RETURN
END Function Integrand   ! --------------------------------------------------

END Module Extra   ! ========================================================
  

PROGRAM MassDistribution                             ! \atmos\mass\masssD.f90
! ---------------------------------------------------------------------------
USE Extra,ONLY: Integrand
USE ForsytheMalcolmMoler,ONLY: Quanc8
IMPLICIT NONE

!  INTEGER,PARAMETER:: WP = SELECTED_REAL_KIND(10,50)  
  REAL:: abserr = 1E-8
  REAL:: errest
  REAL:: flag
  CHARACTER(LEN=*),PARAMETER:: FMT = '(F6.1,2ES12.4, F6.3, I4,F4.0)'
  INTEGER:: k
  REAL:: mass
  REAL,PARAMETER:: MM = 5.29E21 ! limit value of mass (kg)
  INTEGER:: nofun
  REAL:: relerr = 1E-8
  REAL:: z
  REAL,PARAMETER:: ZERO=0.0
!----------------------------------------------------------------------------
  DO k=1,80
    z=REAL(k)
    CALL Quanc8(Integrand, ZERO,z, abserr,relerr, mass, errest, nofun, flag)
    WRITE(*,FMT) z, mass, errest, mass/MM, nofun, flag
  END DO  
  STOP
END Program MassDistribution   ! --------------------------------------------