      SUBROUTINE GRVSTL(PEDNDT,PEDMDT,PIDFDT,PIDNDT,PIDMDT,PENUM,PEMAS, &
                        PEN0,PEPHI0,PEPSI,PEPHIS0,PEDPHI0,PEDDN, &
                        PEWETRC,PIFRC,PINUM,PIMAS,PIN0,PIPHI0, &
                        PIPSI,PIPHIS0,PIDPHI0,PIDDN,PIWETRC, &
                        T,P,DP,RHOA,GRAV,DT,ILGA,LEVA)
!-----------------------------------------------------------------------
!     PURPOSE:
!     --------
!     GRAVITATIONAL SETTLING OF AEROSOL PARTICLES. IT IS NECESSARY
!     TO PROVIDE CONSISTENT REPRESENTATIONS OF THE AEROSOL SIZE
!     DISTRIBUTIONS IN TERMS OF NUMBER AND MASS CONCENTRATIONS,
!     PLA PARAMETERS, AND DRY PARTICLE DENSITY AS INPUT TO THIS SUBROUTINE.
!
!     HISTORY:
!     --------
!     * FEB 19/2010 - K.VONSALZEN   CORRECT PIFRCS
!     * MAR 26/2007 - K.VONSALZEN   MODIFIED FOR EXT. AND INT. MIXED
!                                   AEROSOL SPECIES
!     * DEC 11/2005 - K.VONSALZEN   NEW.
!
!     ************************ INDEX OF VARIABLES **********************
!
!     PRIMARY INPUT/OUTPUT VARIABLES (ARGUMENT LIST):
!
! I   DP       DIFFERENCE IN AIR PRESSURE BETWEEN BOTTOM AND TOP OF GRID
!              CELL (PA)
! I   DT       TIME STEP (SEC)
! I   GRAV     GRAVTIATIONAL ACCELARATION (M/S2)
! I   ILGA     NUMBER OF GRID POINTS IN HORIZONTAL DIRECTION
! I   LEVA     NUMBER OF GRID POINTA IN VERTICAL DIRECTION
! I   P        AIR PRESSURE (PA)
! I   PEDDN    DENSITY OF DRY AEROSOL PARTICLE (KG/M3), EXT. MIXTURE
! O   PEDMDT   MASS TENDENCY (KG/KG/SEC), EXT. MIXTURE
! O   PEDNDT   NUMBER TENDENCY (1/KG/SEC), EXT. MIXTURE
! I   PEDPHI0  SECTION WIDTH, EXT. MIXTURE
! I   PEMAS    AEROSOL MASS CONCENTRATION (KG/KG), EXT. MIXTURE
! I   PEN0     1ST PLA SIZE DISTRIBUTION PARAMETER (AMPLITUDE), EXT.
!              MIXTURE
! I   PENUM    AEROSOL NUMBER CONCENTRATION (1/KG), EXT. MIXTURE
! I   PEPHI0   3RD PLA SIZE DISTRIBUTION PARAMETER (MODE SIZE), EXT.
!              MIXTURE
! I   PEPHIS0  PARTICLE SIZE AT THE LOWER SECTION BOUNDARY, EXT. MIXTURE
! I   PEPSI    2ND PLA SIZE DISTRIBUTION PARAMETER (WIDTH), EXT. MIXTURE
! I   PEWETRC  TOTAL/WET PARICLE RADIUS (M), EXT. MIXTURE
! I   PIDDN    DENSITY OF DRY AEROSOL PARTICLE (KG/M3), INT. MIXTURE
! O   PIDFDT   AEROSOL SPECIES MASS FRACTION TENDENCY (1/SEC), INT.
!              MIXTURE
! O   PIDMDT   MASS TENDENCY (KG/KG/SEC), INT. MIXTURE
! O   PIDNDT   NUMBER TENDENCY (1/KG/SEC), INT. MIXTURE
! I   PIDPHI0  SECTION WIDTH, INT. MIXTURE
! I   PIFRC    AEROSOL SPECIES MASS FRACTION
! I   PIMAS    AEROSOL MASS CONCENTRATION (KG/KG), INT. MIXTURE
! I   PIN0     1ST PLA SIZE DISTRIBUTION PARAMETER (AMPLITUDE), INT.
!              MIXTURE
! I   PINUM    AEROSOL NUMBER CONCENTRATION (1/KG), INT. MIXTURE
! I   PIPHI0   3RD PLA SIZE DISTRIBUTION PARAMETER (MODE SIZE), INT.
!              MIXTURE
! I   PIPHIS0  PARTICLE SIZE AT THE LOWER SECTION BOUNDARY, INT. MIXTURE
! I   PIPSI    2ND PLA SIZE DISTRIBUTION PARAMETER (WIDTH), INT. MIXTURE
! I   PIWETRC  TOTAL/WET PARICLE RADIUS (M), INT. MIXTURE
! I   RH       RELATIVE HUMIDITY (DIMENSIONLESS)
! I   T        AIR TEMPERATURE (K)
!
!     SECONDARY INPUT/OUTPUT VARIABLES (VIA MODULES AND COMMON BLOCKS):
!
! I   ISAEXT   NUMBER OF AEROSOL TRACERS, EXT. MIXTURE
! I   ISAINT   NUMBER OF AEROSOL TRACERS, INT. MIXTURE
! I   KINT     NUMBER OF INTERNALLY MIXED AEROSOL TYPES, INT. MIXTURE
! I   YTINY    SMALLEST VALID NUMBER
!
!-----------------------------------------------------------------------
!
      USE SDPARM
!
      IMPLICIT REAL (A-H,O-Z), INTEGER (I-N)
!
      REAL, INTENT(OUT), DIMENSION(ILGA,LEVA,ISAEXT) :: PEDNDT,PEDMDT
      REAL, INTENT(OUT), DIMENSION(ILGA,LEVA,ISAINT) :: PIDNDT,PIDMDT
      REAL, INTENT(OUT), DIMENSION(ILGA,LEVA,ISAINT,KINT) :: PIDFDT
      REAL, INTENT(IN), DIMENSION(ILGA,LEVA,ISAEXT) :: &
                                         PEN0,PEPHI0,PEPSI,PEWETRC, &
                                         PEDDN,PEPHIS0,PEDPHI0,PENUM, &
                                         PEMAS
      REAL, INTENT(IN), DIMENSION(ILGA,LEVA,ISAINT,KINT) :: PIFRC
      REAL, INTENT(IN), DIMENSION(ILGA,LEVA,ISAINT) :: &
                                         PIN0,PIPHI0,PIPSI,PIWETRC, &
                                         PIDDN,PIPHIS0,PIDPHI0,PINUM, &
                                         PIMAS
      REAL, INTENT(IN), DIMENSION(ILGA,LEVA) :: T,RHOA,DP,P
      INTEGER, INTENT(IN) :: ILGA,LEVA
!
      REAL, ALLOCATABLE, DIMENSION(:,:,:) :: PEFN,PEFM,PEGSVDN, &
                                         PEGSVDM
      REAL, ALLOCATABLE, DIMENSION(:,:,:) :: PIFN,PIFM,PIGSVDN, &
                                         PIGSVDM
      REAL, ALLOCATABLE, DIMENSION(:,:,:,:) :: PIFMS,PIDMDTS,PIMASS, &
                                         PIFRCS
      REAL, ALLOCATABLE, DIMENSION(:,:) :: PESFN,PESFM
      REAL, ALLOCATABLE, DIMENSION(:,:) :: PISFN,PISFM
      REAL, ALLOCATABLE, DIMENSION(:,:,:) :: PISFMS,PIMASST
      REAL, DIMENSION(ILGA,LEVA) :: AMU,AMFP,TERM1,VCFL
      LOGICAL, PARAMETER :: KSRMV=.TRUE.             ! SURFACE REMOVAL
!
!-----------------------------------------------------------------------
!     * ALLOCATE WORK SPACE.
!
      IF ( ISAEXT > 0 ) THEN
        ALLOCATE(PEFN   (ILGA,LEVA,ISAEXT))
        ALLOCATE(PEFM   (ILGA,LEVA,ISAEXT))
        ALLOCATE(PEGSVDN(ILGA,LEVA,ISAEXT))
        ALLOCATE(PEGSVDM(ILGA,LEVA,ISAEXT))
        ALLOCATE(PESFN  (ILGA,ISAEXT))
        ALLOCATE(PESFM  (ILGA,ISAEXT))
      ENDIF
      IF ( ISAINT > 0 ) THEN
        ALLOCATE(PIFRCS (ILGA,LEVA,ISAINT,KINT))
        ALLOCATE(PIMASS (ILGA,LEVA,ISAINT,KINT))
        ALLOCATE(PIDMDTS(ILGA,LEVA,ISAINT,KINT))
        ALLOCATE(PIFMS  (ILGA,LEVA,ISAINT,KINT))
        ALLOCATE(PIMASST(ILGA,LEVA,ISAINT))
        ALLOCATE(PISFMS (ILGA,ISAINT,KINT))
        ALLOCATE(PIFN   (ILGA,LEVA,ISAINT))
        ALLOCATE(PIFM   (ILGA,LEVA,ISAINT))
        ALLOCATE(PIGSVDN(ILGA,LEVA,ISAINT))
        ALLOCATE(PIGSVDM(ILGA,LEVA,ISAINT))
        ALLOCATE(PISFN  (ILGA,ISAINT))
        ALLOCATE(PISFM  (ILGA,ISAINT))
      ENDIF
!
!-----------------------------------------------------------------------
!     * DYNAMIC VISCOSITY OF AIR.
!
      TERM1=SQRT(T)
      TERM1=TERM1*T
      AMU=145.8*1.E-8*TERM1/(T+110.4)
!
!     * MEAN FREE PATH (K.V. BEARD, 1976) AND CUNNINGHAM SLIP
!     * CORRECTION FACTOR.
!
      TERM1=SQRT(T/293.15)
      AMFP=6.54E-8*(AMU/1.818E-5)*(1.013E5/P)*TERM1
!
!-----------------------------------------------------------------------
!     * MAXIMUM FALL VELOCITY ACCORDING TO CFL CRITERION.
!
      VCFL=.999*DP/(RHOA*GRAV*DT)
!
!-----------------------------------------------------------------------
!     * TERMINAL FALL VELOCITIES.
!
      IF ( ISAEXT > 0 ) &
        CALL GSPAR(PEGSVDN,PEGSVDM,RHOA,VCFL,AMU,AMFP,PENUM,PEMAS,PEN0, &
                   PEPHI0,PEPSI,PEPHIS0,PEDPHI0,PEDDN,PEWETRC, &
                   PEDRYRC,GRAV,ILGA,LEVA,ISAEXT)
      IF ( ISAINT > 0 ) &
        CALL GSPAR(PIGSVDN,PIGSVDM,RHOA,VCFL,AMU,AMFP,PINUM,PIMAS,PIN0, &
                   PIPHI0,PIPSI,PIPHIS0,PIDPHI0,PIDDN,PIWETRC, &
                   PIDRYRC,GRAV,ILGA,LEVA,ISAINT)
!
!-----------------------------------------------------------------------
!     * TRACER FLUXES AT GRID CELL INTERFACES.
!
      L=LEVA
      DO IS=1,ISAEXT
        PESFN(:,  IS)=RHOA(:,L)*PEGSVDN(:,L,IS)*PENUM(:,L,IS)
        PEFN (:,:,IS)=RHOA(:,:)*PEGSVDN(:,:,IS)*PENUM(:,:,IS)
        PESFM(:,  IS)=RHOA(:,L)*PEGSVDM(:,L,IS)*PEMAS(:,L,IS)
        PEFM (:,:,IS)=RHOA(:,:)*PEGSVDM(:,:,IS)*PEMAS(:,:,IS)
      ENDDO
      L=LEVA
      DO IS=1,ISAINT
        PISFN(:,  IS)=RHOA(:,L)*PIGSVDN(:,L,IS)*PINUM(:,L,IS)
        PIFN (:,:,IS)=RHOA(:,:)*PIGSVDN(:,:,IS)*PINUM(:,:,IS)
        DO IT=1,KINT
          PISFMS(:,  IS,IT)=RHOA(:,L)*PIGSVDM(:,L,IS) &
                                        *PIFRC(:,L,IS,IT)*PIMAS(:,L,IS)
          PIFMS (:,:,IS,IT)=RHOA(:,:)*PIGSVDM(:,:,IS) &
                                        *PIFRC(:,:,IS,IT)*PIMAS(:,:,IS)
        ENDDO
      ENDDO
!
!-----------------------------------------------------------------------
!     * TENDENCIES FROM UPSTREAM CALCULATIONS.
!
      DO IS=1,ISAEXT
        L=1
        PEDNDT(:,L,IS)=GRAV*(-PEFN(:,L,IS))/DP(:,L)
        PEDMDT(:,L,IS)=GRAV*(-PEFM(:,L,IS))/DP(:,L)
        DO L=2,LEVA-1
          PEDNDT(:,L,IS)=GRAV*(PEFN(:,L-1,IS)-PEFN(:,L,IS))/DP(:,L)
          PEDMDT(:,L,IS)=GRAV*(PEFM(:,L-1,IS)-PEFM(:,L,IS))/DP(:,L)
        ENDDO
      ENDDO
      DO IS=1,ISAINT
        L=1
        PIDNDT(:,L,IS)=GRAV*(-PIFN(:,L,IS))/DP(:,L)
        DO IT=1,KINT
          PIDMDTS(:,L,IS,IT)=GRAV*(-PIFMS(:,L,IS,IT))/DP(:,L)
        ENDDO
        DO L=2,LEVA-1
          PIDNDT(:,L,IS)=GRAV*(PIFN(:,L-1,IS)-PIFN(:,L,IS))/DP(:,L)
          DO IT=1,KINT
            PIDMDTS(:,L,IS,IT)= &
                     GRAV*(PIFMS(:,L-1,IS,IT)-PIFMS(:,L,IS,IT))/DP(:,L)
          ENDDO
        ENDDO
      ENDDO
!
!-----------------------------------------------------------------------
!     * SURFACE REMOVAL. NOTE THAT THIS PART SHOULD BE INCLUDED IN
!     * THE PARAMETERIZATION OF DRY DEPOSITION IN LATER MODEL VERSIONS!
!
      IF ( KSRMV ) THEN
        L=LEVA
        DO IS=1,ISAEXT
          PEDNDT(:,L,IS)=GRAV*(PEFN(:,L-1,IS)-PESFN(:,IS))/DP(:,L)
          PEDMDT(:,L,IS)=GRAV*(PEFM(:,L-1,IS)-PESFM(:,IS))/DP(:,L)
        ENDDO
        DO IS=1,ISAINT
          PIDNDT(:,L,IS)=GRAV*(PIFN(:,L-1,IS)-PISFN(:,IS))/DP(:,L)
          DO IT=1,KINT
            PIDMDTS(:,L,IS,IT)= &
                      GRAV*(PIFMS(:,L-1,IS,IT)-PISFMS(:,IS,IT))/DP(:,L)
          ENDDO
        ENDDO
      ENDIF
!
!-----------------------------------------------------------------------
!     * NEW MASS CONCENTRATIONS FOR INTERNALLY MIXED AEROSOL AFTER
!     * GRAVITATIONAL SETTLING.
!
      IF ( ISAINT > 0 ) THEN
        DO IT=1,KINT
          PIMASS(:,:,:,IT)=PIFRC(:,:,:,IT)*PIMAS(:,:,:) &
                                                  +DT*PIDMDTS(:,:,:,IT)
        ENDDO
        DO IT=1,KINT
          PIMASST=SUM(PIMASS,DIM=4)
          WHERE ( PIMASST > MAX(PIMASS(:,:,:,IT)/YLARGE,YTINY) )
            PIFRCS(:,:,:,IT)=PIMASS(:,:,:,IT)/PIMASST
          ELSEWHERE
            PIFRCS(:,:,:,IT)=PIFRC(:,:,:,IT)
          ENDWHERE
        ENDDO
        PIDMDT=(PIMASST-PIMAS)/DT
        PIDFDT=(PIFRCS-PIFRC)/DT
      ENDIF
!
!-----------------------------------------------------------------------
!     * DEALLOCATE WORK SPACE.
!
      IF ( ISAEXT > 0 ) THEN
        DEALLOCATE(PEFN   )
        DEALLOCATE(PEFM   )
        DEALLOCATE(PEGSVDN)
        DEALLOCATE(PEGSVDM)
        DEALLOCATE(PESFN  )
        DEALLOCATE(PESFM  )
      ENDIF
      IF ( ISAINT > 0 ) THEN
        DEALLOCATE(PIFRCS )
        DEALLOCATE(PIMASS )
        DEALLOCATE(PIDMDTS)
        DEALLOCATE(PIFMS  )
        DEALLOCATE(PIMASST)
        DEALLOCATE(PISFMS )
        DEALLOCATE(PIFN   )
        DEALLOCATE(PIFM   )
        DEALLOCATE(PIGSVDN)
        DEALLOCATE(PIGSVDM)
        DEALLOCATE(PISFN  )
        DEALLOCATE(PISFM  )
      ENDIF
!
      END SUBROUTINE GRVSTL
