      SUBROUTINE COAG(PIDNDT,PIDMDT,PIDFDT,PINUM,PIMAS,PIFRC,PIN0, &
                      PIPHI0,PIPSI,PIDDN,PIWETRC,TA,PA,DT,ILGA,LEVA)
!-----------------------------------------------------------------------
!     PURPOSE:
!     --------
!     BROWNIAN COAGULATION.
!
!     HISTORY:
!     --------
!     * DEC  6/2017 - K.VONSALZEN   IMPROVED SCALING OF LINEARIZED
!                                   SIZE DISTRIBUTION
!     * FEB 10/2010 - K.VONSALZEN   NEWLY INSTALLED.
!
!-----------------------------------------------------------------------
!
      USE SDPARM
      USE SDCODE
      USE COADAT
      USE FPDEF
!
      IMPLICIT REAL (A-H,O-Z), INTEGER (I-N)
!
      REAL, INTENT(OUT), DIMENSION(ILGA,LEVA,ISAINT) :: &
                                                   PIDNDT,PIDMDT
      REAL, INTENT(OUT), DIMENSION(ILGA,LEVA,ISAINT,KINT) :: &
                                                   PIDFDT
      REAL, INTENT(IN), DIMENSION(ILGA,LEVA,ISAINT) :: &
                                                   PIN0,PIPHI0,PIPSI, &
                                                   PIDDN,PIWETRC, &
                                                   PINUM,PIMAS
      REAL, INTENT(IN), DIMENSION(ILGA,LEVA,ISAINT,KINT) :: &
                                                   PIFRC
      REAL, INTENT(IN), DIMENSION(ILGA,LEVA) ::    TA,PA
      REAL, INTENT(IN) :: DT
      INTEGER, PARAMETER :: ISS=1
      INTEGER, PARAMETER ::    IMOM3=1.
      INTEGER, PARAMETER ::    IMOM6=2.
      REAL, DIMENSION(ILGA,LEVA,ISAINT) ::         PIL1,PIL2
      REAL, DIMENSION(ILGA,LEVA,ISAINT) ::         TIMAST,TINUM
      REAL, DIMENSION(ILGA,LEVA,ISAINT,KINT) ::    TIMAS
      REAL(R8), DIMENSION(ILGA,LEVA,ISFINT) ::     FMOM0,FMOM3,FMOM6
      REAL, DIMENSION(ILGA,LEVA,ISFINT) ::         FIDDN,TIMATN,TIMAT, &
                                                   FIFRCS
      REAL, DIMENSION(ILGA,LEVA,ISFINT,KINT) ::    TIMA,FIFRC
      REAL, DIMENSION(ILGA,LEVA,ISFTRI) ::         KERN,CTM
      REAL(R8), DIMENSION(ILGA,LEVA,ISS) ::        FIL1S,FIL2S,VOLLS, &
                                                   VOLRS,SDTMP
      REAL, DIMENSION(ILGA,LEVA) ::                VOLGF,DTC,DNDTT
      REAL, DIMENSION(ILGA,LEVA) ::                ADIFF,ADT,TEMP,PRES, &
                                                   RATM,TOTM1,TOTM2
      REAL(R8), ALLOCATABLE, DIMENSION(:,:,:) ::   FMOM0S,FMOM3S
      REAL(R8), ALLOCATABLE, DIMENSION(:,:,:) ::   FMOM0I,FMOM0D, &
                                                   FMOM3I,FMOM3D
      REAL(R8), ALLOCATABLE, DIMENSION(:,:,:) ::   FIL1,FIL2,FINSD, &
                                                   FIVOLL,FIVOLR,PAR1N, &
                                                   PAR2N,PAR1M,PAR2M, &
                                                   PISCL1,PISCL2
      REAL, ALLOCATABLE, DIMENSION(:,:,:) ::       FIPHIC,TERM,TERM1, &
                                                   FIVOLA,FIVOLW
      REAL, ALLOCATABLE, DIMENSION(:,:,:,:) ::     DNDTS,DNDT,DMDTS
      REAL, ALLOCATABLE, DIMENSION(:,:,:,:,:) ::   DMDT
      LOGICAL, ALLOCATABLE, DIMENSION(:,:,:) ::    PITST
      INTEGER, DIMENSION(ILGA,LEVA) :: IDT,IDP
      INTEGER, DIMENSION(ILGA,LEVA,ISFINT) :: IDR
      REAL(R8) :: AT1,AT2,AT3
!
      REAL, PARAMETER :: TUNE=1.
      INTEGER, PARAMETER :: IMETH=2
!
!-----------------------------------------------------------------------
!     * INITIALIZATION.
!
      IF ( ISAINT > 0 ) THEN
        PIDNDT=0.
        PIDMDT=0.
        PIDFDT=0.
      ENDIF
!
      IF ( KCOAG ) THEN
!
!       * ALLOCATE WORK ARRAYS.
!
        ALLOCATE (FINSD(ILGA,LEVA,ISFINTB))
        ALLOCATE (TERM (ILGA,LEVA,ISFINTB))
        ALLOCATE (TERM1(ILGA,LEVA,ISFINTB))
        ALLOCATE (FIL1   (ILGA,LEVA,ISFINT))
        ALLOCATE (FIL2   (ILGA,LEVA,ISFINT))
!
!-----------------------------------------------------------------------
!       * PARAMETERS FOR PIECEWISE LINEAR REPRESENTATION OF SIZE
!       * DISTRIBUTION. THE USE OF LINEAR DISTRIBUTIONS LEADS TO
!       * MORE EFFICIENT CALCULATIONS IN THE FOLLOWING COMPARED
!       * LOGNORMAL DISTRIBUTIONS DESPITE SOMEWHAT REDUCED ACCURACY.
!
!       * NUMBER SIZE DISTRIBUTION [D N/D R**3] FOR BOUNDARIES OF
!       * COAGULATION SIZE SECTIONS FROM PLA SIZE DISTRIBUTION.
!
        TERM=0.
        ISC=0
        DO IS=1,ISAINT
        DO IST=1,COAGP%ISEC
          ISC=ISC+1
          WHERE ( ABS(PIPSI(:,:,IS)-YNA) > YTINY &
                                        .AND. PIPSI(:,:,IS) < YLARGES )
            TERM(:,:,ISC)=-PIPSI(:,:,IS) &
                                     *(FIPHI(ISC)%VL-PIPHI0(:,:,IS))**2
          ENDWHERE
        ENDDO
        ENDDO
        IS=ISAINT
        ISC=ISFINTB
        WHERE ( ABS(PIPSI(:,:,IS)-YNA) > YTINY &
                                        .AND. PIPSI(:,:,IS) < YLARGES )
          TERM(:,:,ISC)=-PIPSI(:,:,IS) &
                                   *(FIPHI(ISC-1)%VR-PIPHI0(:,:,IS))**2
        ENDWHERE
        TERM1=EXP(TERM)
        FINSD=0._R8
        ISC=0
        DO IS=1,ISAINT
        DO IST=1,COAGP%ISEC
          ISC=ISC+1
          WHERE ( ABS(PIPSI(:,:,IS)-YNA) > YTINY &
              .AND. PINUM(:,:,IS) > YTINY .AND. PIMAS(:,:,IS) > YTINY )
            FINSD(:,:,ISC)=PIN0(:,:,IS)*TERM1(:,:,ISC) &
                                                  /(3._R8*FIDRYVB(ISC))
          ENDWHERE
        ENDDO
        ENDDO
        IS=ISAINT
        ISC=ISFINTB
        WHERE ( ABS(PIPSI(:,:,IS)-YNA) > YTINY &
              .AND. PINUM(:,:,IS) > YTINY .AND. PIMAS(:,:,IS) > YTINY )
          FINSD(:,:,ISC)=PIN0(:,:,IS)*TERM1(:,:,ISC) &
                                                  /(3._R8*FIDRYVB(ISC))
        ENDWHERE
!
!       * PARAMETERS FOR PIECE-WISE LINEAR REPRESENTATION OF NUMBER
!       * SIZE DISTRIBUTION [D N/D R**3 = FIL1 + FIL2 * R**3] FROM
!       * LINEAR INTERPOLATION OF NUMBER SIZE DISTRIBUTION [FINSD]
!       * BETWEEN SECTION BOUNDARIES.
!
        FIL1=0._R8
        FIL2=0._R8
        ISC=0
        DO IS=1,ISAINT
        DO IST=1,COAGP%ISEC
          ISC=ISC+1
          WHERE ( ABS(PIPSI(:,:,IS)-YNA) > YTINY &
              .AND. PINUM(:,:,IS) > YTINY .AND. PIMAS(:,:,IS) > YTINY )
            FIL2(:,:,ISC)=(FINSD(:,:,ISC+1)-FINSD(:,:,ISC)) &
                         /(FIDRYVB(ISC+1)-FIDRYVB(ISC))
            FIL1(:,:,ISC)=FINSD(:,:,ISC)-FIL2(:,:,ISC)*FIDRYVB(ISC)
          ENDWHERE
        ENDDO
        ENDDO
!
!       * DEALLOCATE AND ALLOCATE WORK ARRAYS.
!
        DEALLOCATE (FINSD)
        DEALLOCATE (TERM )
        DEALLOCATE (TERM1)
        ALLOCATE (PAR1N (ILGA,LEVA,ISAINT))
        ALLOCATE (PAR2N (ILGA,LEVA,ISAINT))
        ALLOCATE (PAR1M (ILGA,LEVA,ISAINT))
        ALLOCATE (PAR2M (ILGA,LEVA,ISAINT))
        ALLOCATE (PISCL1(ILGA,LEVA,ISAINT))
        ALLOCATE (PISCL2(ILGA,LEVA,ISAINT))
        ALLOCATE (PITST (ILGA,LEVA,ISAINT))
!
!       * ADJUST PARAMETERS FOR INDIVIDUAL PIECES OF LINEAR DISTRIBUTIONS
!       * SO THAT THE TOTAL NUMBER AND MASS IS CONSERVED WITHIN THE
!       * SIZE RANGES GIVEN BY THE ORIGINAL PLA DISTRIBUTIONS. THIS IS
!       * ACCOMPLISHED BY SCALING THE PARAMETERS FIL1 AND FIL2 LINEARLY.
!       * PAR1N+PAR2N IS THE NUMBER AND PARM1+PARM2 IS THE MASS
!       * CORESSPONDING TO THE LINEAR FIT.
!
        PAR1N=0._R8
        PAR2N=0._R8
        PAR1M=0._R8
        PAR2M=0._R8
        ISC=0
        DO IS=1,ISAINT
        DO IST=1,COAGP%ISEC
          ISC=ISC+1
          AT1=               FIDRYVB(ISC+1)   -FIDRYVB(ISC)
          AT2=.5_R8        *(FIDRYVB(ISC+1)**2-FIDRYVB(ISC)**2)
          AT3=(1._R8/3._R8)*(FIDRYVB(ISC+1)**3-FIDRYVB(ISC)**3)
          PAR1N(:,:,IS)=PAR1N(:,:,IS)+AT1*FIL1(:,:,ISC)
          PAR2N(:,:,IS)=PAR2N(:,:,IS)+AT2*FIL2(:,:,ISC)
          PAR1M(:,:,IS)=PAR1M(:,:,IS)+AT2*FIL1(:,:,ISC)
          PAR2M(:,:,IS)=PAR2M(:,:,IS)+AT3*FIL2(:,:,ISC)
        ENDDO
        ENDDO
        IF ( IMETH == 1 ) THEN
          PITST=.FALSE.
          DO IS=1,ISAINT
          DO L=1,LEVA
          DO IL=1,ILGA
            IF ( PINUM(IL,L,IS) > YTINY .AND. PIMAS(IL,L,IS) > YTINY &
                                                               ) THEN
              PISCL2(IL,L,IS)=PIMAS(IL,L,IS)/(YCNST*PIDDN(IL,L,IS))
              IF ( ABS(PAR1M(IL,L,IS)) &
                > MAX(ABS(PISCL2(IL,L,IS))/YLARGE8,YTINY8) ) THEN
                PISCL2(IL,L,IS)=PISCL2(IL,L,IS)/PAR1M(IL,L,IS)
                PITST(IL,L,IS)=.TRUE.
              ENDIF
              PISCL1(IL,L,IS)=PINUM(IL,L,IS)
              IF ( ABS(PAR1N(IL,L,IS)) > MAX(ABS(PISCL1(IL,L,IS)) &
                 /YLARGE8,YTINY8) .AND. PITST(IL,L,IS) ) THEN
                PISCL1(IL,L,IS)=PISCL1(IL,L,IS)/PAR1N(IL,L,IS)
              ELSE
                PITST(IL,L,IS)=.FALSE.
              ENDIF
              IF ( PITST(IL,L,IS) ) PISCL2(IL,L,IS) &
                 =PISCL2(IL,L,IS)-PISCL1(IL,L,IS)
              IF ( ABS(PAR1M(IL,L,IS)) > MAX(ABS(PAR2M(IL,L,IS))/YLARGE8, &
               YTINY8) .AND. ABS(PAR1N(IL,L,IS)) > MAX(ABS(PAR2N(IL,L,IS)) &
                /YLARGE8,YTINY8) .AND. PITST(IL,L,IS) ) THEN
                PISCL1(IL,L,IS)=PAR2M(IL,L,IS)/PAR1M(IL,L,IS) &
                               -PAR2N(IL,L,IS)/PAR1N(IL,L,IS)
              ELSE
                PITST(IL,L,IS)=.FALSE.
              ENDIF
              IF ( ABS(PISCL1(IL,L,IS)) > MAX(ABS(PISCL2(IL,L,IS)) &
                 /YLARGE8,YTINY8) .AND. PITST(IL,L,IS) ) THEN
                PISCL2(IL,L,IS)=PISCL2(IL,L,IS)/PISCL1(IL,L,IS)
              ELSE
                PITST(IL,L,IS)=.FALSE.
              ENDIF
              PISCL1(IL,L,IS)=PINUM(IL,L,IS) &
                             -PAR2N(IL,L,IS)*PISCL2(IL,L,IS)
              IF ( ABS(PAR1N(IL,L,IS)) > MAX(ABS(PISCL1(IL,L,IS)) &
                 /YLARGE8,YTINY8) .AND. PITST(IL,L,IS) ) THEN
                PISCL1(IL,L,IS)=PISCL1(IL,L,IS)/PAR1N(IL,L,IS)
              ELSE
                PITST(IL,L,IS)=.FALSE.
              ENDIF
            ENDIF
          ENDDO
          ENDDO
          ENDDO
          WHERE ( .NOT.PITST )
            PISCL1=1._R8
            PISCL2=1._R8
          ENDWHERE
        ELSE
          WHERE ( PAR1N+PAR2N > MAX(ABS(PINUM)/YLARGE8,YTINY8) )
            PISCL1=PINUM/(PAR1N+PAR2N)
          ELSEWHERE
            PISCL1=1._R8
          ENDWHERE
          WHERE ( PAR1M+PAR2M > MAX(ABS(PIMAS/(YCNST*PIDDN)) &
                                    /YLARGE8,YTINY8) )
            PISCL2=(PIMAS/(YCNST*PIDDN))/(PAR1M+PAR2M)
          ELSEWHERE
            PISCL2=1._R8
          ENDWHERE
          PISCL1=.5*(PISCL1+PISCL2)
          PISCL2=PISCL1
        ENDIF
        ISC=0
        DO IS=1,ISAINT
        DO IST=1,COAGP%ISEC
          ISC=ISC+1
          FIL1(:,:,ISC)=FIL1(:,:,ISC)*PISCL1(:,:,IS)
          FIL2(:,:,ISC)=FIL2(:,:,ISC)*PISCL2(:,:,IS)
        ENDDO
        ENDDO
!
!       * DEALLOCATE AND ALLOCATE WORK ARRAYS.
!
        DEALLOCATE (PAR1N )
        DEALLOCATE (PAR2N )
        DEALLOCATE (PAR1M )
        DEALLOCATE (PAR2M )
        DEALLOCATE (PISCL1)
        DEALLOCATE (PISCL2)
        DEALLOCATE (PITST )
        ALLOCATE (FIPHIC (ILGA,LEVA,ISFINT))
        ALLOCATE (FMOM0S (ILGA,LEVA,ISFTRI))
        ALLOCATE (FMOM3S (ILGA,LEVA,ISFTRI))
        ALLOCATE (FIVOLA (ILGA,LEVA,ISFTRI))
        ALLOCATE (FIVOLW (ILGA,LEVA,ISFINTB))
        ALLOCATE (FIVOLL (ILGA,LEVA,ISFINT))
        ALLOCATE (FIVOLR (ILGA,LEVA,ISFINT))
!
!       * MAP PLA SIZE DISTRIBUTION ONTO APPROPRIATE SIZE DISTRIBUTION
!       * FOR COAGULATION CALCULATIONS. IN THE FOLLOWING, THE SIZE
!       * DISTRIBUTION FOR COAGULATION CALCULATIONS IS FOR A DOUBLING
!       * OF PARTICLE VOLUME BETWEEN SECTION BOUNDARIES AND FOR WET
!       * PARTICLE SIZES.
!
        CALL SDMAP(FIL1,FIL2,FIDDN,FIFRC,FIVOLW,FIVOLA,FIPHIC,VOLGF, &
                   PIDDN,PIFRC,PIWETRC,PIMAS,PIPSI,ILGA,LEVA)
!
!       * TEMPERATURE, PRESSURE, AND PARTICLE SIZE CORRESPONDING
!       * TO TABULATED COAGULATION DATA.
!
        ADIFF=YLARGE
        DO ITT=1,IDA
          ADT=ABS(TA-CTEM(ITT))
          WHERE ( ADT < ADIFF )
            IDT=ITT
            ADIFF=ADT
          ENDWHERE
        ENDDO
        ADIFF=YLARGE
        DO ITT=1,IDD
          ADT=ABS(PA-CPRE(ITT))
          WHERE ( ADT < ADIFF )
            IDP=ITT
            ADIFF=ADT
          ENDWHERE
        ENDDO
        DO IS=1,ISFINT
          ADIFF=YLARGE
          DO ITT=1,IDB
            ADT=ABS(FIPHIC(:,:,IS)-CPHI(ITT))
            WHERE ( ADT < ADIFF )
              IDR(:,:,IS)=ITT
              ADIFF=ADT
            ENDWHERE
          ENDDO
        ENDDO
!
!       * COAGULATION COEFFICIENT (KERNEL).
!
        ITRI=0
        DO IS=1,ISFINT
        DO IK=1,IS
          ITRI=ITRI+1
          DO L=1,LEVA
          DO IL=1,ILGA
            KERN(IL,L,ITRI)=TUNE &
                 *CKER(IDR(IL,L,IS),IDR(IL,L,IK),IDT(IL,L),IDP(IL,L))
          ENDDO
          ENDDO
        ENDDO
        ENDDO
!
!       * MOMENTS OF SIZE DISTRIBUTION FOR COMPLETE SECTIONS.
!
        FIVOLL=FIVOLW(:,:,1:ISFINT )
        FIVOLR=FIVOLW(:,:,2:ISFINTB)
        FMOM0=SDINTL0(FIL1,FIL2,      FIVOLL,FIVOLR,ILGA,LEVA,ISFINT)
        FMOM3=SDINTL (FIL1,FIL2,IMOM3,FIVOLL,FIVOLR,ILGA,LEVA,ISFINT)
        FMOM6=SDINTL (FIL1,FIL2,IMOM6,FIVOLL,FIVOLR,ILGA,LEVA,ISFINT)
!
!       * MOMENTS OVER PARTS OF SECTIONS FOR CALCULATION OF MAIN TERMS.
!
        ITRI=0
        DO IS=1,ISFINT
        DO IK=1,IS
          ITRI=ITRI+1
          VOLRS(:,:,ISS)=FIVOLW(:,:,IS+1)
          VOLLS(:,:,ISS)=MIN(FIVOLA(:,:,ITRI),FIVOLW(:,:,IS+1))
          FIL1S(:,:,ISS)=FIL1(:,:,IS)
          FIL2S(:,:,ISS)=FIL2(:,:,IS)
          SDTMP=SDINTL0(FIL1S,FIL2S,VOLLS,VOLRS,ILGA,LEVA,ISS)
          FMOM0S(:,:,ITRI)=SDTMP(:,:,ISS)
          SDTMP=SDINTL (FIL1S,FIL2S,IMOM3,VOLLS,VOLRS,ILGA,LEVA,ISS)
          FMOM3S(:,:,ITRI)=SDTMP(:,:,ISS)
        ENDDO
        ENDDO

!
!       * DEALLOCATE WORK ARRAYS.
!
        DEALLOCATE (FIL1   )
        DEALLOCATE (FIL2   )
        DEALLOCATE (FIPHIC )
        DEALLOCATE (FIVOLL )
        DEALLOCATE (FIVOLR )
        DEALLOCATE (FIVOLA )
        ALLOCATE (FMOM0I(ILGA,LEVA,ISFTRIM))
        ALLOCATE (FMOM0D(ILGA,LEVA,ISFTRIM))
        ALLOCATE (FMOM3I(ILGA,LEVA,ISFTRIM))
        ALLOCATE (FMOM3D(ILGA,LEVA,ISFTRIM))
!
!       * COEFFICIENTS FOR LINEAR VOLUME INTERPOLATION OF PARTIAL MOMENTS
!       * BETWEEN ADJACENT SECTION BOUNDARIES.
!
        ITRI=0
        ITRIM=0
        DO IS=1,ISFINT
          DO IK=1,IS-1
            ITRI=ITRI+1
            ITRIM=ITRIM+1
            FMOM0I(:,:,ITRIM)=FMOM0S(:,:,ITRI)-FIVOLW(:,:,IK) &
                                     /(FIVOLW(:,:,IK+1)-FIVOLW(:,:,IK)) &
                            *(FMOM0S(:,:,ITRI+1)-FMOM0S(:,:,ITRI))
            FMOM3I(:,:,ITRIM)=FMOM3S(:,:,ITRI)-FIVOLW(:,:,IK) &
                                     /(FIVOLW(:,:,IK+1)-FIVOLW(:,:,IK)) &
                            *(FMOM3S(:,:,ITRI+1)-FMOM3S(:,:,ITRI))
            FMOM0D(:,:,ITRIM)=(FMOM0S(:,:,ITRI+1) &
                   -FMOM0S(:,:,ITRI))/(FIVOLW(:,:,IK+1)-FIVOLW(:,:,IK))
            FMOM3D(:,:,ITRIM)=(FMOM3S(:,:,ITRI+1) &
                   -FMOM3S(:,:,ITRI))/(FIVOLW(:,:,IK+1)-FIVOLW(:,:,IK))
          ENDDO
          ITRI=ITRI+1
        ENDDO
!
!       * DEALLOCATE WORK ARRAYS.
!
        DEALLOCATE (FMOM0S)
        DEALLOCATE (FMOM3S)
        DEALLOCATE (FIVOLW)
        ALLOCATE (DNDTS(ILGA,LEVA,ISFINT,ISAINT))
        ALLOCATE (DNDT (ILGA,LEVA,ISAINT,ISAINT))
!
!-----------------------------------------------------------------------
!       * NUMBER TENDENCY EQUATION.
!
        CALL COAGN(DNDTS,KERN,FMOM0,FMOM3,FMOM0I,FMOM0D,ILGA,LEVA)
!
!-----------------------------------------------------------------------
!       * SUM UP TENDENCIES TO MATCH PLA SECTIONS.
!
        DNDT=0.
        ISC=0
        DO IS=1,ISAINT
        DO IST=1,COAGP%ISEC
          ISC=ISC+1
          WHERE ( ABS(DNDTS(:,:,ISC,:)) > YTINY )
            DNDT(:,:,IS,:)=DNDT(:,:,IS,:)+DNDTS(:,:,ISC,:)
          ENDWHERE
        ENDDO
        ENDDO
!
!       * DEALLOCATE AND ALLOCATE WORK ARRAYS.
!
        DEALLOCATE (DNDTS)
        ALLOCATE (DMDT (ILGA,LEVA,ISAINT,ISAINT,KINT))
        ALLOCATE (DMDTS(ILGA,LEVA,ISFINT,ISAINT))
!
!-----------------------------------------------------------------------
!       * MASS TENDENCY EQUATION. SUM UP TENDENCIES TO MATCH
!       * PLA SECTIONS.
!
        DO IR1=1,ISFINT
        DO IR2=1,IR1
          ITS=ITR(IR1,IR2)
          CTM(:,:,ITS)=KERN(:,:,ITS)*FMOM0(:,:,IR1)*FMOM3(:,:,IR2)
        ENDDO
        ENDDO
        DMDT=0.
        DO IT=1,KINT
          FIFRCS(:,:,:)=FIFRC(:,:,:,IT)
          CALL COAGM(DMDTS,CTM,KERN,FIFRCS,FIDDN,FMOM0,FMOM3,FMOM6, &
                     FMOM0I,FMOM0D,FMOM3I,FMOM3D,VOLGF,ILGA,LEVA)
          ISC=0
          DO IS=1,ISAINT
          DO IST=1,COAGP%ISEC
            ISC=ISC+1
            WHERE ( ABS(DMDTS(:,:,ISC,:)) > YTINY )
              DMDT(:,:,IS,:,IT)=DMDT(:,:,IS,:,IT)+DMDTS(:,:,ISC,:)
            ENDWHERE
          ENDDO
          ENDDO
        ENDDO
!
!       * DEALLOCATE WORK ARRAYS.
!
        DEALLOCATE (DMDTS)
        DEALLOCATE (FMOM0I)
        DEALLOCATE (FMOM0D)
        DEALLOCATE (FMOM3I)
        DEALLOCATE (FMOM3D)
!
!-----------------------------------------------------------------------
!       * INITIAL MASS SIZE DISTRIBUTION.
!
        DO IT=1,KINT
          TIMAS(:,:,:,IT)=PIMAS(:,:,:)*PIFRC(:,:,:,IT)
        ENDDO
!        DO IT=1,KINT
!          TIMAS(:,:,:,IT)=0.
!          ISC=0
!          DO IS=1,ISAINT
!          DO IST=1,COAGP%ISEC
!            ISC=ISC+1
!            WHERE ( ABS(PIPSI(:,:,IS)-YNA) > YTINY )
!              TIMAS(:,:,IS,IT)=TIMAS(:,:,IS,IT)+FIFRC(:,:,ISC,IT)
!     1                  *YCNST*FIDDN(:,:,ISC)*FMOM3(:,:,ISC)/VOLGF(:,:)
!            ENDWHERE
!          ENDDO
!          ENDDO
!        ENDDO
!
!-----------------------------------------------------------------------
!       * ADJUST TIME STEP, IF NECESSARY. THE ADJUSTMENT IS PERFORMED
!       * FOR EACH INDIVIDUAL SECTION (OUTER LOOP) BASED ON THE GATHERED
!       * TENDENCIES FOR COAGULATION OF PARTICLES IN THAT SECTION WITH
!       * PARTICLES WITH EQUAL OR GREATER SIZES TO ACCOUNT FOR DIFFERENCES
!       * IN COAGULATION EFFICIENCIES DEPENDING ON PARTICLE SIZES.
!       * GENERALLY, COAGULATION WILL BE MOST EFFICIENT FOR THE SMALLEST
!       * PARTICLES RESULTING IN POTENTIALLY NECESSARY REDUCTIONS IN
!       * TIME STEP TO AVOID NEGATIVE CONCENTRATIONS FROM SIMPLE EULER
!       * FORWARD INTEGRATION. HOWEVER, LITTLE OR NO ADJUSTMENT MAY BE
!       * REQUIRED FOR LARGER PARTICLES SO THAT THERE MAY ONLY BE MINOR
!       * ADJUSTMENTS TO TIME STEPS FOR COAGULATION OF THE LARGER
!       * PARTICLES.
!
        TINUM=PINUM
        DO ISC=ISAINT,1,-1
          DTC=DT/YSEC
          DO IS=1,ISAINT
          DO L=1,LEVA
          DO IL=1,ILGA
            IF ( DNDT(IL,L,IS,ISC) < -MAX(TINUM(IL,L,IS)/YLARGE,YTINY) &
                                        .AND. DTC(IL,L) > YSMALL ) THEN
              IF ( -TINUM(IL,L,IS)/DTC(IL,L) > DNDT(IL,L,IS,ISC) ) THEN
                DTC(IL,L)=-TINUM(IL,L,IS)/DNDT(IL,L,IS,ISC)
              ENDIF
!              DTC(IL,L)=MAX(MIN(DTC(IL,L),
!     1                     -TINUM(IL,L,IS)/DNDT(IL,L,IS,ISC)),0.)
            ENDIF
            DO IT=1,KINT
              IF ( DMDT(IL,L,IS,ISC,IT) &
                                 < -MAX(TIMAS(IL,L,IS,IT)/YLARGE,YTINY) &
                                        .AND. DTC(IL,L) > YSMALL ) THEN
                IF ( -TIMAS(IL,L,IS,IT)/DTC(IL,L) &
                                          > DMDT(IL,L,IS,ISC,IT) ) THEN
                  DTC(IL,L)=-TIMAS(IL,L,IS,IT)/DMDT(IL,L,IS,ISC,IT)
                ENDIF
!                DTC(IL,L)=MAX(MIN(DTC(IL,L),
!     1                     -TIMAS(IL,L,IS,IT)/DMDT(IL,L,IS,ISC,IT)),0.)
              ENDIF
            ENDDO
          ENDDO
          ENDDO
          ENDDO
          DTC=MAX(MIN(DT,DTC*YSEC),0.)
          DO IS=1,ISAINT
            WHERE ( ABS(PIPSI(:,:,IS)-YNA) > YTINY )
              DNDT(:,:,IS,ISC)=DNDT(:,:,IS,ISC)*DTC(:,:)/DT
              TINUM(:,:,IS)=TINUM(:,:,IS)+DNDT(:,:,IS,ISC)*DT
            ENDWHERE
            DO IT=1,KINT
              WHERE ( ABS(PIPSI(:,:,IS)-YNA) > YTINY )
                DMDT(:,:,IS,ISC,IT)=DMDT(:,:,IS,ISC,IT)*DTC(:,:)/DT
                TIMAS(:,:,IS,IT)=TIMAS(:,:,IS,IT)+DMDT(:,:,IS,ISC,IT)*DT
              ENDWHERE
            ENDDO
          ENDDO
        ENDDO
!
!       * REJECT NEGATIVE SOLUTIONS FROM NUMERICAL TRUNCATION.
!
        TINUM=MAX(TINUM,0.)
        TIMAS=MAX(TIMAS,0.)
!
!       * MASS FIXER TO COMPENSATE FOR TRUNCATION ERRORS.
!
        DO IT=1,KINT
          TOTM1=SUM(PIMAS(:,:,:)*PIFRC(:,:,:,IT),DIM=3)
          TOTM2=SUM(TIMAS(:,:,:,IT),DIM=3)
          WHERE ( TOTM2 > MAX(TOTM1/YLARGE,YTINY) )
            RATM=TOTM1/TOTM2
          ELSEWHERE
            RATM=1.
          ENDWHERE
          DO IS=1,ISAINT
            TIMAS(:,:,IS,IT)=TIMAS(:,:,IS,IT)*RATM(:,:)
          ENDDO
        ENDDO
!
!       * DEALLOCATE WORK ARRAYS.
!
        DEALLOCATE (DNDT )
        DEALLOCATE (DMDT )
!
!       * PARTICLE NUMBER COAGULATION TENDECY FROM DIFFERENCE
!       * OF ADJUSTED NEW AND ORIGINAL PARTICLE NUMBER CONCENTRATION.
!
        WHERE ( ABS(PIPSI-YNA) > YTINY )
          PIDNDT=(TINUM-PINUM)/DT
        ENDWHERE
!
!       * PARTICLE MASS COAGULATION TENDECY FROM DIFFERENCE
!       * OF ADJUSTED NEW AND ORIGINAL PARTICLE MASS CONCENTRATION.
!
        TIMAST=0.
        DO IT=1,KINT
          TIMAST(:,:,:)=TIMAST(:,:,:)+TIMAS(:,:,:,IT)
        ENDDO
        DO IT=1,KINT
          WHERE ( TIMAST > MAX(TIMAS(:,:,:,IT)/YLARGE,YTINY) &
                                         .AND. ABS(PIPSI-YNA) > YTINY )
            PIDFDT(:,:,:,IT)=(TIMAS(:,:,:,IT)/TIMAST(:,:,:) &
                                                   -PIFRC(:,:,:,IT))/DT
          ENDWHERE
        ENDDO
        WHERE ( ABS(PIPSI-YNA) > YTINY )
          PIDMDT=(TIMAST-PIMAS)/DT
        ENDWHERE
      ENDIF
!
      END SUBROUTINE COAG
