C
C  NEW 2D MODEL Routine for computing  RAINOUT, sets trop values to Oort, 83  --  EF  4/18/94
C   NO MORE RELATIVE HUMIDITY FIX at TROPICAL TROPOPAUSE NEEDED WITH PPM
C   Also, Oort values at mid-high latitudes seem too wet much of the year compared to saturation
C   vapor pressure over ice based on the NMC temps, so rainout also at levels 3-4, 
C
C    NOTE: Oort values for NOW USED FOR ONLY LEVELS 1-2, although this is very similar to setting 
C          levels 1-4, since its determined by the RH Rainout limit......
C
C
        SUBROUTINE RAINOUT


#include "com2d.h"

C
C   ITRAIN IS THE INDEX OF THE TOP MODEL LEVEL TO COMPUTE RAINOUT FOR EACH LATITUDE, -- NO (OLD)
C
C   For BASE8SB, now have A RAINOUT PROFILE FOR the TROPICS (15S-15N - needs to be somewhat wet),
C   and outside the tropics (dry), and interpolate at 25-35N,S,
C   so relative humidity limit (RHRAIN) is a function of latitude and height
C   ****  NOW specified based on the daily climalogical trop hgts for up to 9 levels in the 
C         tropics, lower levels at mid-high latitudes, uses ITROP360(L$,360) in COMMON *******
C
C   Note: need to make stratosphere/mesosphere wetter with new PPM transport, so raise RH level, 
C   keep constant w/ height at 45% (60% too wet), note levels 1-2 are not used
C

         DIMENSION ITRAIN(L$), RHRAIN(L$,9), RHTR(9), RHML(9)
C
CC         DATA ITRAIN/7,7,7,7, 8,8, 9,9,9, 9,9,9, 8,8, 7,7,7,7/
CC  OLD       DATA ITRAIN/7,7,8,9,10,11,11,11,11,11,11,11,11,10,9,8,7,7/

         DATA RHTR/.48, .48, .48, .48, .48, .48, .48, .48, .48/
         DATA RHML/.20, .20, .20, .20, .20, .20, .20, .20, .20/


         DO 1001 ik=1,9
            DO 1005 ij=1,5
 1005          RHRAIN(ij,ik) = RHML(ik)

            DO 1010 ij=6,7
 1010          RHRAIN(ij,ik) = (RHTR(ik)-RHML(ik))/3.*(ij-5) + RHML(ik)

            DO 1015 ij=8,11
 1015          RHRAIN(ij,ik) = RHTR(ik)

            DO 1020 ij=12,13
 1020          RHRAIN(ij,ik) = (RHML(ik)-RHTR(ik))/3.*(ij-14) + RHML(ik)
   
            DO 1025 ij=14,18
 1025          RHRAIN(ij,ik) = RHML(ik)
 1001  CONTINUE

C
C   Area of the globe is in  CM-2
C
        GLOBE = 5.1E18 
        IDAY = IDAY360
	DRPLFT = DT       
C
C
C    SET TROPOSPHERIC VALUES (Level 1-2) TO THOSE OF OORT, 1983, from DAILY input values
C
        DO 8889 IK=1,2
        DO 8889 IJ=1,18
8889    CN(15,IJ,IK) = TROPWV(IJ,IK,IDAY)*M(IJ,IK)/1.E06


C  Now compute rainout, SET INDEX OF LOWEST MODEL LEVEL TO COMPUTE RAINOUT  EG, 3 = 498 mb
C   and rainout up to the NMC tropopause hgt (ITROP360) + 1 since with only a 2km grid spacing
C   we'll miss some of the troposphere, especially in the tropics
C
        IBW = 3


        DO 666 IJ = 1,18
        DO 666 IK = IBW, ITROP360(IJ,IDAY360)+1

CC        DO 666 IK = IBW,ITRAIN(IJ)
c
c   calculate LOCAL saturation vapor pressure over ice from analytic form since temps in region
C      of rainout are 180-239K, and convert from Torr (mm of Hg) to mb => 1013./760.

	TVAP = TEMP(IJ,IK)
        VPS = 1013./760.*exp(-88.44375824 + 0.7641223392*tvap -
     c             0.002285011125*(tvap**2) + 2.539292763e-6*(tvap**3))
C
C   PRECIP IS THE NUMBER OF H2O MOLECULES WHICH ARE MORE THAN the Rel. hum. limit OF THE
C      H2O MOLECULES IN SATURATED AIR - then subtract from the total density,
C      .622 = ratio of molecular weight of water to dry air
C 
       	XSAT = .622*VPS/PRESS(IK)
      	WVLIM = RHRAIN(ij,ik)*XSAT*M(IJ,IK)

      	PRECIP = CN(15,IJ,IK) - WVLIM
        PRECIP = AMAX1(PRECIP,0.0)
        
        CN(15,IJ,IK) = CN(15,IJ,IK) - PRECIP

	RAIN=PRECIP/DRPLFT
        PRECIPYR(IJ,IK,IDAY) = PRECIP/1.E20

666     CONTINUE

C
C   Now compute rainfall for a year for each latitude, and total # of H2O 
C       molecules rained out in a year for the entire globe for LAST 2 YEARS OF RUN
C          
       IF (ITDAY .GE. INDAYS-360-359) THEN

           DO 201 ILAT=1,18
           SUMZ = 0.0
           DO 200 IZ = IBW,11
200        SUMZ = SUMZ + (PRECIPYR(ILAT,IZ,IDAY)*DELTAZ(ILAT,IZ)*1.E5)
201        RAINDAY(ILAT,IDAY) = SUMZ
C
C        SUML = 0.0
C        DO 203 ILAT=1,18
C203     SUML = SUML + (RAINDAY(ILAT,IDAY)*AREA(ILAT)*GLOBE*0.5)
C        TOTRAINDAY = SUML
C        WRITE(52,709)
C709    FORMAT(/,5X,'LATITUDE  RAINOUT  1.E+20 # mole / cm2')
C       WRITE(52,707) (RAINDAY(IL,IDAY),IL=1,18)
C        ENDIF
C

      IF (IDAY .EQ. 360) THEN
        DO 668 JLAT=1,18
C 
        SUM = 0.0
        DO 667 JDAY=1,360
667     SUM = SUM + RAINDAY(JLAT,JDAY)
C
668     RAINYR(JLAT) = SUM
C
        TOTH2O = 0.0
        DO 670 ILAT=1,18
670     TOTH2O = TOTH2O + (RAINYR(ILAT)*AREA(ILAT)*GLOBE*0.5)
C
       WRITE(52,810) IYR, DAY360, IDAY, ITDAY
810    FORMAT(//,10X,'IYR =',I4,7X,'DAY360 =',F6.1,5X,
     c      'IDAY=',I4,5X,'ITDAY =',I6)
C
        WRITE(52,811) 
811     FORMAT(/,3X,'TOTAL RAINOUT for the year at each Latitude is',
     *  '  --  E+20 # molecules H2O/cm-2 :')
        WRITE(52,707) (RAINYR(IJ),IJ=1,18)
707     FORMAT(5X,1P6E14.4)
C
        WRITE(52,813) TOTH2O 
813   FORMAT(/,1X,'Total  RAINOUT for the globe',
     *' for 1 year  -- E+20 # molecules H2O = ',1PE14.4,/)
C
        DO 901 JDAY=1,360
        DO 901 ILAT=1,18
 901    RAINDAY(ILAT,JDAY) = 0.0
        
        DO 902 ILAT=1,18
 902    RAINYR(ILAT) = 0.0
C

      ENDIF
      ENDIF


      IF (IDAY .EQ. 360) THEN
        DO 910 KD = 1,360 
        DO 910 JD = 1,46
        DO 910 ID = 1,18
 910    PRECIPYR(ID,JD,KD) = 0.0
      ENDIF
C

        SAVE
	RETURN
	END