	PROGRAM SOUNDING
C*****************************************************************************
C
C  	This program reads in a sounding of temperature and specific humidity,
C  allowing the user to specify the tropospheric pressure levels for the
C  model. A standard atmosphere sounding may be used; this is employed in
C  any case to specify the inital sounding above 100 mb and to specify the
C  distribution of ozone. The radiative cooling profile can also be
C  specified, if interactive radiation is not used, and the vertical 
C  velocity can likewise be specified if desired.  
C
C*****************************************************************************
C
C  Dimension various arrays
C
	PARAMETER (NS=150)
	REAL P(NS), T(NS), R(NS), PRAD(NS), TRAD(NS), RRAD(NS)
	REAL CO2(NS), O3(NS), OMEGA(NS), RCOOL(NS)
C
	CHARACTER*1 STDA, SKIPA
C
	OPEN(UNIT=8,FILE='sounding.in', STATUS='NEW')
C
C  Initialize
C
	DO 2 I=1,NS
	 OMEGA(I)=0.0
	 RCOOL(I)=0.0
    2	CONTINUE
C
C  Ask the user whether to skip reading in the thermodynamic sounding
C
	PRINT*, 'Do you want to read in the thermodynamic sounding'
	PRINT*, 'from the file sounding.out and skip to the specification'
	PRINT*, 'of vertical velocity and/or radiative cooling? (y or n)'
	READ(*,4) SKIPA
    4	FORMAT(A1)
	IF(SKIPA.EQ.'y'.OR.SKIPA .EQ.'Y')THEN
	 OPEN(UNIT=9,FILE='sounding.out',STATUS='OLD')
	 READ(9,5)NP,TS
    5 	 FORMAT(T15,I3,T55,F5.2,////)
	 DO 8 I=1,NP
	  READ(9,7)P(I),T(I),R(I),O3(NP+1-I),OMEGA(I),RCOOL(I)
    7 	  FORMAT(3X,F6.1,T12,F8.3,T25,F10.6,T39,F6.3,T47,F6.2,
     1     T58,F7.3)
	  O3(NP+1-I)=1.0E-6*O3(NP+1-I)
    8	 CONTINUE 
	 CLOSE(UNIT=9,DISPOSE='SAVE')
	 GOTO 100
	END IF
C
C  Input sea surface temperature
C
	PRINT*, ' '
	PRINT*, 'Enter sea surface temperature (C)'
	READ*, TS
C
C  Ask the user to decide whether to use standard atmosphere
C
	PRINT*, ' '
	PRINT*, 'Do you want to use standard atmosphere for whole'
	PRINT*, '     sounding? (y or n)'
	READ(*,4) STDA
C
C  If yes, read in standard atmosphere
C
	IF(STDA.EQ.'y' .OR.STDA.EQ.'Y')THEN
	 N=0
   10	 N=N+1
	 PRINT*, ' '
	 PRINT*, 'Enter pressure level (mb) starting at surface,'
	 PRINT*, '   Enter any number less than 100.0 when finished'
	 READ*, PN
	 IF(PN.GE.100.0)THEN
	  P(N)=PN
	  GOTO 10
	 END IF 
	  N=N-1
	  NP=N+9 
	 PRINT*, ' '
	 PRINT*, 'Total number of levels=',NP
	 GOTO 40 
	END IF
C
C  If no, ask the user for sounding
C
	PRINT*, ' '
	PRINT*, 'Enter 0 if you want to enter specific humidity; 1 if you'
	PRINT*, 'prefer to enter dew point temperatures'
	READ*, MM
	N=0
	IF(MM.EQ.0)THEN
   20	 N=N+1
	 PRINT*, ' '
	 PRINT*, 'Enter pressure (mb), temperature (C) and specific'
	 PRINT*, '  humidity (g/kg) starting at lowest level.'
	 PRINT*, '  (Enter any number less than 100.0 for pressure to'
	 PRINT*, '  finish)'
	 READ*, PN, TN, RN
	 IF(PN.GE.100.0)THEN
	  P(N)=PN
	  R(N)=RN
	  T(N)=TN
	  GOTO 20
	 END IF
	 N=N-1
	 NP=N+9 
	 PRINT*, ' '
	 PRINT*, 'Total number of levels=',NP
	ELSE
   30	 N=N+1
	 PRINT*, ' '
	 PRINT*, 'Enter pressure (mb), temperature (C), and dew point'
	 PRINT*, '  temperature (C) starting at lowest level.'
	 PRINT*, '  (Enter any number less than 100.0 for pressure to'
	 PRINT*, '  finish)'
	 READ*, PN, TN, TDN
	 IF(PN.GE.100.0)THEN
	  P(N)=PN
	  T(N)=TN
	  EV=6.112*EXP(17.67*TDN/(243.5+TDN))
	  EV=MIN(EV,(PN-0.1))
	  R(N)=622.0*EV/(PN-EV)
	  GOTO 30
	 END IF
	 N=N-1
	 NP=N+9 
	 PRINT*, ' '
	 PRINT*, 'Total number of levels=',NP
	END IF
   40	CONTINUE
C
C  Predefine stratospheric pressure levels
C
	P(N+1)=80.0
	P(N+2)=60.0
	P(N+3)=40.0
	P(N+4)=30.0
	P(N+5)=25.0
	P(N+6)=20.0
	P(N+7)=15.0
	P(N+8)=10.0
	P(N+9)=5.0
C
C  Invert pressure array for standard atmosphere, which uses reverse order
C
	DO 50 I=1,NP
	 PRAD(I)=P(NP+1-I)
   50	CONTINUE
C
C  Call standard atmosphere subroutine to get standard T, R, CO2, O3
C
	PRINT*,' '
	PRINT*,' '
	PRINT*, 'Standard atmosphere: Enter 1 for tropical atmos.,'
	PRINT*, ' 2 for midlatitude summer, 3 for midlatitude'
	PRINT*, ' winter, 4 for subarctic summer, or 5 for'
	PRINT*, ' subarctic winter'
	READ*, IATM
	CALL STDATM(PRAD(1),PRAD(NP),NP,IATM,PRAD,TRAD,RRAD,CO2,O3)
	CO21=CO2(NP)
C
C  If standard atmosphere used in troposphere, reset T and R
C
	IF(STDA.EQ.'y'.OR. STDA.EQ.'Y')THEN
	 DO 60 I=1,NP
	  T(I)=TRAD(NP+1-I)-273.15
	  R(I)=1000.0*RRAD(NP+1-I)
   60	 CONTINUE
C
C  Otherwise, reset T and R only in stratosphere
C
	ELSE
	 DO 70 I=N+1,NP
	  T(I)=TRAD(NP+1-I)-273.15	
	  R(I)=1000.0*RRAD(NP+1-I)
   70	CONTINUE
	END IF
  100	CONTINUE
C
C  Ask user whether to input omega and/or radiative cooling profiles
C
	PRINT*, ' '
	PRINT*, 'Do you want to enter just an omega profile (enter 1),'
	PRINT*, ' just a radiative cooling profile (enter 2),'
	PRINT*,'  both (enter 3), or neither (enter 4)?'
	READ*, IPRO
C
C  Input omega and/or radiative cooling profiles
C
	IF(IPRO.EQ.1)THEN
	 DO 110 I=1,NP
	  PRINT*, ' '
	  PRINT*, 'Enter omega (mb/hr) at p=',P(I),'  mb'
	  READ*, OMEGA(I)
  110	 CONTINUE
	ELSE IF(IPRO.EQ.2)THEN
	 DO 120 I=1,NP
	  PRINT*, ' '
	  PRINT*, 'Enter radiative cooling (deg/day) at p=',P(I),'  mb' 
	  READ*, RCOOL(I)
  120	 CONTINUE
	ELSE IF(IPRO.EQ.3)THEN
	 DO 130 I=1,NP
	  PRINT*, ' '
	  PRINT*, 'Enter omega (mb/hr) and radiative cooling (deg/day)'
	  PRINT*, '  at p=',P(I),'  mb'
	  READ*, OMEGA(I), RCOOL(I)
  130	 CONTINUE
	END IF
C
C  Output to file sounding.in
C
	WRITE(8,140) NP,TS
  140	FORMAT(2X, 'N levels = ',T15,I3,
     1    T48,'SST=',T55,F5.2,' C',/)
C
	WRITE(8,150)
  150	FORMAT(2X,'Pressure',T14,'Temp.',T25,'Spec. humid.',T40,
     1    'O3 mr',T48, 'Omega',T56,'Rad. cooling')
	WRITE(8,160)
  160	FORMAT(4X,'(mb)',T15,'(C)',T28,'(g/kg)',T40,'10**-6',T47,
     1    '(mb/hr)',T57,'(deg/day)')
  	WRITE(8,170)
  170	FORMAT(2X,'-------',T14,'-----',T25,'------------',T40,
     1   '-----',T48, '-----',T56,'------------')
C
	DO 200 I=1,NP
	 OZ=1.0E6*O3(NP+1-I)
	 WRITE(8,180)P(I),T(I),R(I),OZ,OMEGA(I),RCOOL(I)
  180	 FORMAT(3X,F6.1,T12,F8.3,T25,F10.6,T39,F6.3,T47,F6.2,
     1   T58,F7.3)
  200	CONTINUE	
	CLOSE(UNIT=8,DISPOSE='SAVE')
	STOP
	END
************************************************************************
*  FILE: STDATM.FOR                          FIRST VERSION: JUNE/1984. *
*  ----------------                          MODIFIED     : SEPT/1990. *
*  Created by Ridgway (N. Renno)                                       *
*----------------------------------------------------------------------*
*     This  routine provides standard atmosphere profiles (AFGL).      *  
************************************************************************

      SUBROUTINE STDATM(PMIN,PMAX,NL,IATM,P,T,W,C,O)

*----------------------------------------------------------------------*
* REQUIRED-                                                            *
* --------                                                             *
*     PMIN : MINIMUM PRESSURE (MB) WHERE PROFILES ARE TO BE LOADED     *
*     PMAX : MAXIMUM PRESSURE (MB) WHERE PROFILES ARE TO BE LOADED     *
*     NL   : NUMBER OF LEVELS .LE. ML=100                              *
*     IATM : CHOICE OF ATMOSPHERIC PROFILE 1-5, THEY ARE:              *
*            IATM = 1  TROPICAL ATM                                    *
*                 = 2  MIDLATITUDE SUMMER                              *
*                 = 3  MIDLATITUDE WINTER                              *
*                 = 4  SUBARCTIC SUMMER                                *
*                 = 5  SUBARCTIC WINTER                                *
*                                                                      *
*     P(NL): PRESSURE (MB) LEVELS IN INCREASING ORDER                  *
*                                                                      *
* RETURNED-                                                            *
* --------                                                             *
*     T(NL): TEMPERATURES (K) AT PRESSURE LEVELS P()                   *
*     W(NL): WATER VAPOR (GM/GM) MASS MIXING RATIOS AT P()             *
*     C(NL): CO2   ( PPMV) MIXING RATIOS AT P()                        *
*     O(NL): OZONE (GM/GM) MASS MIXING RATIO AT P()                    *
*----------------------------------------------------------------------*

      PARAMETER (ML=100)
      REAL    P(*),T(*),W(*),O(*),C(*)
      REAL    PUP(ML),ZUP(ML),TUP(ML),WUP(ML),OUP(ML)
      DATA    CO2 / 330./
*
*     PRINT *,' .... STDATM.......PRESSURE LEVELS  ',PMIN,'  TO',PMAX
*
      IF(NL.GT.ML) STOP  'BECAUSE NL.GT.ML IN STDATM'
      IF(IATM.LE.0) STOP 'BECAUSE IATM.LE.0 IN STDATM'
      IF(IATM.GT.5) STOP 'BECAUSE IATM.GT.5 IN STDATM'
*
      DO 1 L=1,NL
1     PUP(L)=P(NL+1-L)
*
*  FDTEMP AND FDH2O COMPUTE ALTITUDE, TEMPERATURE, AND WATER
*  VAPOR DENSITIES AT SELECTED PRESSURE LEVELS.  FIRST SET ZUP(1).
*
      ZUP(1)=0.
      CALL FDTEMP(IATM,NL,PUP,ZUP,TUP)
      CALL FDH2O (IATM,NL,PUP,TUP,WUP)
      CALL FINDO3(IATM,NL,PUP,TUP,OUP)
*
*  PROFILE LEVELS ARE SURFACE-TO-SPACE.
*  ABSORPTION PRESSURE LEVELS ARE SPACE-TO-SURFACE.
*
      DO 100 L=1,NL
      IF(PMIN.LE.P(L).AND.P(L).LE.PMAX) THEN
      P(L)=PUP(NL+1-L)
      T(L)=TUP(NL+1-L)
      W(L)=WUP(NL+1-L)
      C(L)=CO2
      O(L)=OUP(NL+1-L)
      END IF
100   CONTINUE
      RETURN
      END

*----------------------------------------------------------------------*
      SUBROUTINE FDTEMP(L,NL,PRESS,ALT,TEMP)
*
*     FELS TEMPERATURE INTERPOLATION PROGRAM AS MODIFIED
*     FOR USE IN COMPUTING CLEAR ATMOSPHERE FLUXES
*
*     JUNE 15, 1984                    (RIDGWAY)
*
*
*  CALCULATES ALTITUDES AND TEMPERATURES AT SELECTED PRESSURES.
*  SIX ANALYTICAL FUNCTIONS APPROXIMATE AFGL STANDARD ATMS.
*
*
*  USAGE NOTES :
*
*     ARRAYS PRESS, ALT, AND TEMP MUST HAVE DIMENSIONS
*            OF AT LEAST (NL) IN CALLING PROGRAM.
*
*     ALT(1) SHOULD BE SET TO MATCH PRESS(1) IN CALLING PROGRAM.
*
*     PRESS ARRAY MUST BE IN ORDER OF DECREASING PRESSURE,
*            THAT IS, INCREASING ALTITUDE.
*
*
*
*  ARGUMENT DEFINITIONS :
*
*
*      L       :  ATMOSPHERIC MODEL INDEX.
*      NL      :  NUMBER OF PRESSURE LEVELS DEFINED.
*      PRESS(1):  WILL BE RESET EQUAL TO SURF PRESS IF GREATER
*      PRESS   :  CALLING PROGRAM-SUPPLIED SET OF PRESSURES.
*      ALT(1)  :  ALTITUDE CORRESPONDING TO FIRST PRESSURE LEVEL.
*      ALT     :  RETURNED ALTITUDE LEVELS IN KM.
*      TEMP    :  RETURNED MODEL TEMPERATURES AT CHOSEN LEVELS.
*
*
      DIMENSION PRESS(1),ALT(1),TEMP(1)
*
      DIMENSION PS(6)
      DATA PS /1013., 1013., 1018., 1010., 1013., 1013./
*
      DELZAP=0.5
      R=8.314
      G0=9.8
      ZMASS=28.9
      AA=6400.
*
*     ALT(1)=0.0
*     FIRST ALTITUDE NOW SET BY CALLING PROGRAM
*
      TEMP(1)=ANTEMP(L,ALT(1))
*     FIRST TEMPERATURE SET AT FIRST ALTITUDE.
*
*     SET FIRST PRESSURE LEVEL TO THAT AT SURFACE
      IF(PRESS(1).GT.PS(L)) PRESS(1)=PS(L)
*
      NLM1=NL-1
      DO 100 N=1,NLM1
*
*  PUT IN COMPUTATION LEVELS EVERY "DELZAP" KM.   ???
*
      DLOGP=7.0*LOG(PRESS(N)/PRESS(N+1))
      NINT=DLOGP/DELZAP
      NINT=NINT+1
      ZNINT=NINT
      G=G0*(1.0-2.0*ALT(N)/AA)
      DZ=R*DLOGP/(7.0*ZMASS*G*ZNINT)
      HT=ALT(N)
*
*  CALCULATE HT AT NEXT USER LEVEL BY R-K INTEGRATION
*
      DO 200 M=1,NINT
      RK1=ANTEMP(L,HT)*DZ
      RK2=ANTEMP(L,HT+0.5*RK1)*DZ
      RK3=ANTEMP(L,HT+0.5*RK2)*DZ
      RK4=ANTEMP(L,HT+RK3)*DZ
      HT=HT+(RK1+RK2+RK2+RK3+RK3+RK4)/6.0
200   CONTINUE
      ALT(N+1)=HT
      TEMP(N+1)=ANTEMP(L,HT)
100   CONTINUE
*
      RETURN
      END

*----------------------------------------------------------------------*
      FUNCTION ANTEMP(L,Z)
      DIMENSION ZB(10,6),C(11,6),DELTA(10,6),TSTAR(6)
      DIMENSION ZB1(10),ZB2(10),ZB3(10),ZB4(10),ZB5(10),ZB6(10)
      DIMENSION C1(11),C2(11),C3(11),C4(11),C5(11),C6(11)
      DIMENSION DELTA1(10),DELTA2(10),DELTA3(10),
     &          DELTA4(10),DELTA5(10),DELTA6(10)
      EQUIVALENCE (ZB(1,1),ZB1(1)),(C(1,1),C1(1)),
     &            (DELTA(1,1),DELTA1(1))
      EQUIVALENCE (ZB(1,2),ZB2(1)),(C(1,2),C2(1)),
     &            (DELTA(1,2),DELTA2(1))
      EQUIVALENCE (ZB(1,3),ZB3(1)),(C(1,3),C3(1)),
     &            (DELTA(1,3),DELTA3(1))
      EQUIVALENCE (ZB(1,4),ZB4(1)),(C(1,4),C4(1)),
     &            (DELTA(1,4),DELTA4(1))
      EQUIVALENCE (ZB(1,5),ZB5(1)),(C(1,5),C5(1)),
     &            (DELTA(1,5),DELTA5(1))
      EQUIVALENCE (ZB(1,6),ZB6(1)),(C(1,6),C6(1)),
     &            (DELTA(1,6),DELTA6(1))
*
* TROPICAL SOUNDING
*
      DATA ZB1
     &            /  2.0,  3.0, 16.5, 21.5, 45.0,
     &              51.0, 70.0,100.0,200.0,300.0   /
*
      DATA C1
     &            / -6.0, -4.0, -6.7,  4.0,  2.2,
     &               1.0, -2.8,-0.27,  0.0,  0.0,  0.0   /
*
      DATA DELTA1
     &            /  0.5,  0.5,  0.3,  0.5,  1.0,
     &               1.0,  1.0,  1.0,  1.0,  1.0   /
*
* SUB-TROPICAL SUMMER
*
      DATA ZB2
     &            /  1.5,  6.5, 13.0, 18.0, 26.0,
     &              36.0, 48.0, 50.0, 70.0,100.0   /
*
      DATA C2
     &            / -4.0, -6.0, -6.5,  0.0,  1.2,
     &               2.2,  2.5,  0.0, -3.0,-0.25, 0.0   /
*
      DATA DELTA2
     &            /  0.5,  1.0,  0.5,  0.5,  1.0,
     &               1.0,  2.5,  0.5,  1.0,  1.0   /
*
* SUB-TROPICAL WINTER
*
      DATA ZB3
     &            /  3.0, 10.0, 19.0, 25.0, 32.0,
     &              44.5, 50.0, 71.0, 98.0,200.0   /
*
      DATA C3
     &            / -3.5, -6.0, -0.5,  0.0,  0.4,
     &               3.2,  1.6, -1.8, -0.7,  0.0, 0.0   /
*
      DATA DELTA3
     &            /  0.5,  0.5,  1.0,  1.0,  1.0,
     &               1.0,  1.0,  1.0,  1.0,  1.0   /
*
* SUB-ARCTIC SUMMER
*
      DATA ZB4
     &            /  4.7, 10.0, 23.0, 31.8, 44.0,
     &              50.2, 69.2,100.0,102.0,103.0   /
*
      DATA C4
     &            / -5.3, -7.0,  0.0,  1.4,  3.0,
     &               0.7, -3.3, -0.2,  0.0,  0.0, 0.0   /
*
      DATA DELTA4
     &            /  0.5,  0.3,  1.0,  1.0,  2.0,
     &               1.0,  1.5,  1.0,  1.0,  1.0   /
*
* SUB-ARCTIC WINTER
*
      DATA ZB5
     &            /  1.0,  3.2,  8.5, 15.5, 25.0,
     &              30.0, 35.0, 50.0, 70.0,100.0   /
*
      DATA C5
     &            /  3.0, -3.2, -6.8,  0.0, -0.6,
     &               1.0,  1.2,  2.5, -0.7, -1.2, 0.0   /
*
      DATA DELTA5
     &            /  0.4,  1.5,  0.3,  0.5,  1.0,
     &               1.0,  1.0,  1.0,  1.0,  1.0   /
*
* US STANDARD 1976
*
      DATA ZB6
     &            / 11.0, 20.1, 32.1, 47.4, 51.4,
     &              71.7, 85.7, 90.0, 91.0, 92.0   /
*
      DATA C6
     &            / -6.5,  0.0,  1.0, 2.75,  0.0,
     &             -2.75,-1.97,  0.0,  0.0,  0.0, 0.0   /
*
      DATA DELTA6
     &            /  0.3,  1.0,  1.0,  1.0,  1.0,
     &               1.0,  1.0,  1.0,  1.0,  1.0   /
*
      DATA TSTAR  /  300.0, 294.0, 272.2, 287.0, 257.1, 288.15/
*
      IF(L.LE.6) GO TO 1
      ANTEMP=FLOAT(L)
*
*  L > 6  REPRESENTS AN INTEGER TEMPERATURE, NOT MODEL NUMBER.
*
      RETURN
*
1     NLAST=10
*
*  L  IS MODEL NUMBER  (1-6)
*
      TEMP=TSTAR(L)+C(1,L)*Z
      DO 20 N=1,NLAST
      EXPO=(Z-ZB(N,L))/DELTA(N,L)
      EXPO=MIN( 50.0,EXPO)
      EXPO=MAX(-50.0,EXPO)
      EXPP=ZB(N,L)/DELTA(N,L)
      EXPP=MIN( 50.0,EXPP)
      EXPP=MAX(-50.0,EXPP)
      FAC=EXP(EXPP)+EXP(-EXPP)
      TEMP=TEMP+(C(N+1,L)-C(N,L))*0.5*(Z+DELTA(N,L)*
     &      LOG((EXP(EXPO)+EXP(-EXPO))/FAC))
20    CONTINUE
      ANTEMP=TEMP
      RETURN
      END

*----------------------------------------------------------------------*
      SUBROUTINE FDH2O(L,NLEV,PRESS,TEMP,DEN)
*
*        L      =   MCCLATCHEY ATMOSPHERE INDEX
*        NLEV   =   NUMBER OF PRESSURE LEVELS USED
*        PRESS  =   ARRAY OF PRESSURES ANY ORDER
*        TEMP   =   TEMPERATURES AT PRESSURE LEVELS
*        DEN    =   WATER VAPOR MASS RATIO IN  G/KG
*
*  ARRAYS MUST BE SET UP FOR NLEV ELEMENTS IN CALLING PROGRAM
*
      DIMENSION PRESS(*),TEMP(*),DEN(*)
*
      DO 100 N=1,NLEV
      DEN(N)=0.
      IF(L.EQ.0) GO TO 100
      CALL QGET(PRESS(N),L,Q,QMR)
*
*  TO GET MIXING RATIO IN GM/KG USING IN CALL TO DENSTY
****  IF(Q.GT.0.) DEN(N)=Q/DENSTY(PRESS(N),TEMP(N),0.)
*  TO GET MIXING RATIO IN GM/GM USING IN CALL TO DENSTY
      IF(Q.GT.0.) DEN(N)=0.001*Q/DENSTY(PRESS(N),TEMP(N),0.)
*
*  ACTUAL MIXING RATIO (G/G) RETURNED IN STRATOSPHERE :
*
*  TO CONVERT TO GM / KG :
****  IF(QMR.GT.0.) DEN(N)=1.E+3*QMR
*  TO REMAIN IN  GM / GM :
      IF(QMR.GT.0)  DEN(N)= QMR
100   CONTINUE
*
      RETURN
      END

*----------------------------------------------------------------------*
      SUBROUTINE QGET(P,L,Q,QMR)
*
*  P   : PRESSURE IN MBAR
*  L   : MCLATCHEY ATMOSPHERE INDEX
*      = 1  TROPICAL ATM
*        2  MIDLATITUDE SUMMER
*        3  MIDLATITUDE WINTER
*        4  SUBARCTIC SUMMER
*        5  SUBARCTIC WINTER
*  Q   : WATER VAPOR DENSITY IN G/M**3
*  QMR : WATER VAPOR STRATOSPHERIC MASS MIXING RATIO MH2O/MTOTAL
*            QMR IS DIMENSIONLESS (GM/GM)
*
      DIMENSION A(6,5),B(6,5),QS(5),PRAT(6,5),PS(5)
      DIMENSION PMIN(5),QCONST(5)
*
      DATA QS   / 19., 14., 3.5, 9.1, 1.2  /
*
      DATA A    / 3.55, 1.37, 3.99, 2.79, 8.66, 4.66,
     &            3.31, 4.77, 4.19, 1.30, 8.21, 4.52,
     &            2.72, 2.11, 3.41, 5.44, 3.91, 0.00,
     &            3.76, 2.95, 0.83, 4.71, 7.43, 4.28,
     &            0.00,-0.59, 4.99, 4.18, 3.30, 0.00   /
*
      DATA B    / 1.91,-7.58,-0.41,-1.52, 1.81, 0.00,
     &           -1.83, 0.67, 0.00,-2.26, 1.67, 0.00,
     &            0.39,-2.02,-0.81, 0.71, 0.00, 0.00,
     &            2.37,-1.01,-3.17, 0.00, 1.41, 0.00,
     &            0.00,-4.46, 0.86, 0.38, 0.00, 0.00   /
*
      DATA PRAT /0.795, 0.694, 0.339, 0.172, 0.110, 0.000,
     &            0.558, 0.421, 0.278, 0.172, 0.110, 0.000,
     &            0.776, 0.342, 0.263, 0.100, 0.000, 0.000,
     &            0.787, 0.375, 0.294, 0.145, 0.107, 0.000,
     &            0.876, 0.350, 0.185, 0.100, 0.000, 0.000  /
*
      DATA PS   / 1013., 1013., 1018., 1010., 1013.  /
*
      DATA PMIN / 110.99, 152.99, 117.79, 124.99, 110.29  /
*
      DATA QCONST / 3.25, 4.0, 4.0, 4.0, 4.0  /
*
      Q=0.
      QMR=0.
      PR=P/PS(L)
      N=1
      IF(P.LE.PMIN(L)) GO TO 25
15    IF(PR.GE.PRAT(N,L)) GO TO 20
      N=N+1
      GO TO 15
20    ETA=A(N,L) + B(N,L)*LOG(PR)
      Q=QS(L)*PR**ETA
      RETURN
25    QMR=1.E-6*QCONST(L)
      RETURN
      END

*----------------------------------------------------------------------*
      FUNCTION DENSTY(P,T,R)
*
*     DENSTY   :   ATMOSPHERIC DENSITY (KG/M**3)
*     P        :   TOTAL PRESSURE (MBAR)
*     T        :   TEMPERATURE (K)
*     R        :   DIMENSIONLESS MASS MIXING RATIO (GM/GM)
*
      DATA RD /287.05 /
      DENSTY = 100.*P/(RD*VIRTEM(T,R))
      RETURN
      END

*----------------------------------------------------------------------*
      FUNCTION VIRTEM(T,R)
*
*     VIRTEM   :   VIRTUAL TEMPERATURE (K)
*     R        :   DIMENSIONLESS MASS MIXING RATIO (GM/GM)
*
      DATA EPS /0.622 /
      VIRTEM = T * (1.+R/EPS)/(1.+R)
      RETURN
      END

*----------------------------------------------------------------------*
      SUBROUTINE FINDO3(L,NLEVEL,PDUM,TDUM,RDUM)
*
*    VERSION OF JULY 17, 1984             (RIDGWAY)
*
*  RETURNS THE OZONE MIXING RATIO  R  IN  G/KG USING A
*  LINEAR INTERPOLATION IN MIXING RATIO VS. LOG(P).
*
*              L  :   ATMOSPHERIC MODEL NUMBER 1-6
*                     L=0  NO OZONE
*
*              P  :   PRESSURE IN MBAR
*              T  :   TEMPERATURE IN DEGREES KELVIN
*
*              R  :   OZONE MIXING RATIO   G/KG
*
*
* **************************************************************
*
*  ******* ONLY MODEL 2 *** MIDLATTITUDE SUMMER IS AVAILABLE
*
* **************************************************************
*
      IMPLICIT REAL (A-H,O-Z)
      REAL  PDUM(1), TDUM(1), RDUM(1)
*
      REAL PREF(33,6)
      EQUIVALENCE (PREF(1,1),PREF1(1))
      EQUIVALENCE (PREF(1,2),PREF2(1))
      EQUIVALENCE (PREF(1,3),PREF3(1))
      EQUIVALENCE (PREF(1,4),PREF4(1))
      EQUIVALENCE (PREF(1,5),PREF5(1))
      EQUIVALENCE (PREF(1,6),PREF6(1))
*
      REAL PREF1(33) /
     &             1013.,   904.,   805.,   715.,   633.,
     &              559.,   492.,   432.,   378.,   329.,
     &              286.,   247.,   213.,   182.,   156.,
     &              132.,   111.,   93.7,   78.9,   66.0,
     &              56.5,   48.0,   40.9,   35.0,   30.0,
     &              25.7,   12.2,   6.00,   3.05,   1.59,
     &              0.854,  5.79E-2,  3.00E-4               /
*
      REAL PREF2(33) /
     &             1013.,   902.,   802.,   710.,   628.,
     &              554.,   487.,   426.,   372.,   324.,
     &              281.,   243.,   209.,   179.,   153.,
     &              130.,   111.,   95.0,   81.2,   69.5,
     &              59.5,   51.0,   43.7,   37.6,   32.2,
     &              27.7,   13.2,   6.52,   3.33,   1.76,
     &              .951,6.71E-2, 3.0E-4                   /
*
      REAL PREF3(33) /
     &             1018.,  897.3,  789.7,  693.8,  608.1,
     &             531.3,  462.7,  401.6,  347.3,  299.2,
     &             256.8,  219.9,  188.2,  161.0,  137.8,
     &             117.8,  100.7,  86.10,  73.50,  62.80,
     &             53.70,  45.80,  39.10,  33.40,  28.60,
     &             24.30,  11.10,  5.180,  2.530,  1.290,
     &             .6820,4.67E-2, 3.0E-4                   /
*
      REAL PREF4(33) /
     &             1010.,  896.0,  792.9,  700.0,  616.0,
     &             541.0,  473.0,  413.0,  359.0,  310.7,
     &             267.7,  230.0,  197.7,  170.0,  146.0,
     &             125.0,  108.0,  92.80,  79.80,  68.60,
     &             58.90,  50.70,  43.60,  37.50,  32.27,
     &             27.80,  13.40,  6.610,  3.400,  1.810,
     &             .9870,7.07E-2, 3.0E-4                   /
*
      REAL PREF5(33) /
     &             1013.,  887.8,  777.5,  679.8,  593.2,
     &             515.8,  446.7,  385.3,  330.8,  282.9,
     &             241.8,  206.7,  176.6,  151.0,  129.1,
     &             110.3,  94.31,  80.58,  68.82,  58.75,
     &             50.14,  42.77,  36.47,  31.09,  26.49,
     &             22.56,  10.20,  4.701,  2.243,  1.113,
     &             .5719,  4.016E-2,  3.0E-4               /
*
*   MODEL  6   NOT CODED YET ********
*
      REAL PREF6(33) /
     &             1013.,   902.,   802.,   710.,   628.,
     &              554.,   487.,   426.,   372.,   324.,
     &              281.,   243.,   209.,   179.,   153.,
     &              130.,   111.,   95.0,   81.2,   69.5,
     &              59.5,   51.0,   43.7,   37.6,   32.2,
     &              27.7,   13.2,   6.52,   3.33,   1.76,
     &              .951,6.71E-2, 3.0E-4                   /
*
      REAL AMTO3(33,6)
      EQUIVALENCE (AMTO3(1,1),AMT1(1))
      EQUIVALENCE (AMTO3(1,2),AMT2(1))
      EQUIVALENCE (AMTO3(1,3),AMT3(1))
      EQUIVALENCE (AMTO3(1,4),AMT4(1))
      EQUIVALENCE (AMTO3(1,5),AMT5(1))
      EQUIVALENCE (AMTO3(1,6),AMT6(1))
*
      REAL AMT1(33) /
     &            5.6E-5, 5.6E-5, 5.4E-5, 5.1E-5, 4.7E-5,
     &            4.5E-5, 4.3E-5, 4.1E-5, 3.9E-5, 3.9E-5,
     &            3.9E-5, 4.1E-5, 4.3E-5, 4.5E-5, 4.5E-5,
     &            4.7E-5, 4.7E-5, 6.9E-5, 9.0E-5, 1.4E-4,
     &            1.9E-4, 2.4E-4, 2.8E-4, 3.2E-4, 3.4E-4,
     &            3.4E-4, 2.4E-4, 9.2E-5, 4.1E-5, 1.3E-5,
     &            4.3E-6, 8.6E-8, 4.3E-11                   /
*
      REAL AMT2(33) /
     &            6.0E-5, 6.0E-5, 6.0E-5, 6.2E-5, 6.4E-5,
     &            6.6E-5, 6.9E-5, 7.5E-5, 7.9E-5, 8.6E-5,
     &            9.0E-5, 1.1E-4, 1.2E-4, 1.5E-4, 1.8E-4,
     &            1.9E-4, 2.1E-4, 2.4E-4, 2.8E-4, 3.2E-4,
     &            3.4E-4, 3.6E-4, 3.6E-4, 3.4E-4, 3.2E-4,
     &            3.0E-4, 2.0E-4, 9.2E-5, 4.1E-5, 1.3E-5,
     &            4.3E-6, 8.6E-8, 4.3E-11                   /
*
      REAL AMT3(33) /
     &            6.0E-5, 5.4E-5, 4.9E-5, 4.9E-5, 4.9E-5,
     &            5.8E-5, 6.4E-5, 7.7E-5, 9.0E-5, 1.2E-4,
     &            1.6E-4, 2.1E-4, 2.6E-4, 3.0E-4, 3.2E-4,
     &            3.4E-4, 3.6E-4, 3.9E-4, 4.1E-4, 4.3E-4,
     &            4.5E-4, 4.3E-4, 4.3E-4, 3.9E-4, 3.6E-4,
     &            3.4E-4, 1.9E-4, 9.2E-5, 4.1E-5, 1.3E-5,
     &            4.3E-6, 8.6E-8, 4.3E-11                   /
*
      REAL AMT4(33) /
     &            4.9E-5, 5.4E-5, 5.6E-5, 5.8E-5, 6.0E-5,
     &            6.4E-5, 7.1E-5, 7.5E-5, 7.9E-5, 1.1E-4,
     &            1.3E-4, 1.8E-4, 2.1E-4, 2.6E-4, 2.8E-4,
     &            3.2E-4, 3.4E-4, 3.9E-4, 4.1E-4, 4.1E-4,
     &            3.9E-4, 3.6E-4, 3.2E-4, 3.0E-4, 2.8E-4,
     &            2.6E-4, 1.4E-4, 9.2E-5, 4.1E-5, 1.3E-5,
     &            4.3E-6, 8.6E-8, 4.3E-11                   /
*
      REAL AMT5(33) /
     &            4.1E-5, 4.1E-5, 4.1E-5, 4.3E-5, 4.5E-5,
     &            4.7E-5, 4.9E-5, 7.1E-5, 9.0E-5, 1.6E-4,
     &            2.4E-4, 3.2E-4, 4.3E-4, 4.7E-4, 4.9E-4,
     &            5.6E-4, 6.2E-4, 6.2E-4, 6.2E-4, 6.0E-4,
     &            5.6E-4, 5.1E-4, 4.7E-4, 4.3E-4, 3.6E-4,
     &            3.2E-4, 1.5E-4, 9.2E-5, 4.1E-5, 1.3E-5,
     &            4.3E-6, 8.6E-8, 4.3E-11                   /
*
*   MODEL   6    NOT CODED YET ********
*
      REAL AMT6(33) /
     &            6.0E-5, 6.0E-5, 6.0E-5, 6.2E-5, 6.4E-5,
     &            6.6E-5, 6.9E-5, 7.5E-5, 7.9E-5, 8.6E-5,
     &            9.0E-5, 1.1E-4, 1.2E-4, 1.5E-4, 1.8E-4,
     &            1.9E-4, 2.1E-4, 2.4E-4, 2.8E-4, 3.2E-4,
     &            3.4E-4, 3.6E-4, 3.6E-4, 3.4E-4, 3.2E-4,
     &            3.0E-4, 2.0E-4, 9.2E-5, 4.1E-5, 1.3E-5,
     &            4.3E-6, 8.6E-8, 4.3E-11                   /
*
      IF(L.GT.5) STOP ' STDATM CALLED FINDO3 : MODEL NUMBER 1 - 5 ONLY'
*
      DO 30 JLEVEL=1,NLEVEL
      RDUM(JLEVEL)=0.D0
      IF(L.LE.0) GO TO 30
      P=PDUM(JLEVEL)
      T=TDUM(JLEVEL)
      J2=1
      DO 10 J=1,32
      IF(P.GE.PREF(J,L)) GO TO 20
10    J2=J
20    J1=J2+1
*
*  LOGORYTHMIC PRESSURE SCALE
*
      X1=LOG(PREF(J1,L))
      X2=LOG(PREF(J2,L))
      X =LOG(P)
*
*  LINEAR INTERPOLATION OF MIXING RATIO AT SINGLE TEMPERATURE
*
      R1=AMTO3(J1,L)/DENSTY(PREF(J1,L),T,0.)
      R2=AMTO3(J2,L)/DENSTY(PREF(J2,L),T,0.)
*
      R=R1+(X-X1)*(R2-R1)/(X2-X1)
*
*   LOAD ARRAY AND CONVERT   GM -- KG     1.E-3
*                            PPMM         1.E+6
*                            PPMV         28.9644 / 47.9982
*
* TO GET AMOUNT IN PPMV :
*     RDUM(JLEVEL)=R*1.E+3*28.9644/47.9982
* TO GET AMOUNT IN GM/GM:
      RDUM(JLEVEL)=R*1.E-3
*
30    CONTINUE
      RETURN
      END
*----------------------------------------------------------------------*
