      SUBROUTINE XTEAERO4(XEMIS,PRESSG,DSHJ,ZF,
     1                    EBBTROW,EOBTROW,EBFTROW,EOFTROW,ESPTROW,
     2                    ESITROW,ESDTROW,ESSTROW,ESOTROW,ESRTROW,
     3                    IPAM,DOCEM1,DBCEM1,DOCEM2,DBCEM2,DSUEM2,
     4                    DSUEM3,PSUEF2,PSUEF3,
     5                    IOCO,IOCY,IBCO,IBCY,ISO2,
     6                    ILG,IL1,IL2,ILEV,LEV,NTRAC,ZTMST)
C
C     * JUN 17/2013 - K.VONSALZEN/ NEW VERSION FOR GCM17+:
C     *               M.LAZARE.    - ADD CONDITION TO SET WGT=1. FOR
C     *                              SPECIAL CASE WHERE HSO2=LSO2.
C     *                            - ADD SUPPORT FOR PLA OPTION.
C     * APR 28/2012 - K.VONSALZEN/ PREVIOUS VERSION XTEAERO3 FOR GCM16:
C     *               M.LAZARE.    - INITIALIZE LSO2 AND HSO2 TO
C     *                              ILEV INSTEAD OF ZERO, TO
C     *                              AVOID POSSIBLE INVALID 
C     *                              MEMORY REFERENCES.
C     *                            - BUGFIX FOR INCORRECT MEMORY
C     *                              REFERENCE: LSO2(L)->LSO2(IL).
C     * APR 26/2010 - K.VONSALZEN. PREVIOUS VERSION XTEAERO2 FOR GCM15I:
C     *                            - ZF AND PF CALCULATED ONE-TIME AT
C     *                              BEGINNING OF PHYSICS AND ZF
C     *                              PASSED IN, WITH INTERNAL CALCULATION
C     *                              THUS REMOVED AND NO NEED TO PASS
C     *                              IN SHBJ,T.
C     *                            - EMISSIONS ACCUMULATED INTO GENERAL
C     *                              "XEMIS" ARRAY RATHER THAN UPDATING
C     *                              XROW (XROW UPDATED FROM XEMIS IN
C     *                              THE PHYSICS). 
C     *                            - REMOVAL OF DIAGNOSTIC FIELDS
C     *                              CALCULATION, SO NO NEED TO PASS
C     *                              IN "SAVERAD" OR "ISVCHEM".
C     *                            - WGT=0 (IE NO EMISSIONS IN LAYER)
C     *                              IF IT IS ABOVE TOP OR BELOW BASE
C     *                              OF EMISSIONS.
C     * FEB 05/09 - K.VONSALZEN. PREVIOUS VERSION XTEAERO1 FOR GCM15H:
C                                - REMOVE NATURAL EMISSIONS.
C     * JAN 18/08 - X.MA/     PREVIOUS VERSION XTEAERO FOR GCM15G.
C                   M.LAZARE. 
C
      IMPLICIT REAL (A-H,O-Z),
     +INTEGER (I-N)
C
C     * I/O FIELDS.
C
      REAL XEMIS(ILG,ILEV,NTRAC)
      REAL DSHJ(ILG,ILEV),PRESSG(ILG)
      REAL ZF(ILG,ILEV)
      REAL, DIMENSION(ILG,ILEV) :: DOCEM1, DBCEM1, DOCEM2, DBCEM2, 
     1                             DSUEM2, DSUEM3
C
      REAL EBBTROW(ILG),EOBTROW(ILG),EBFTROW(ILG),EOFTROW(ILG)
      REAL ESPTROW(ILG),ESITROW(ILG)
      REAL ESDTROW(ILG),ESSTROW(ILG),ESOTROW(ILG),ESRTROW(ILG)
C
C     * INTERNAL WORK FIELDS.
C
      INTEGER, DIMENSION(ILG) ::  LSO2,HSO2 
C
      COMMON /PARAMS/ WW,     TW,    A,     ASQ,  GRAV, RGAS,  RGOCP,
     1                RGOASQ, CPRES, RGASV, CPRESV
C
      REAL, PARAMETER :: WSO2    =64.059
      REAL, PARAMETER :: WH2SO4  =98.073
      REAL, PARAMETER :: WNH3    =17.030
      REAL, PARAMETER :: WAMSUL  =WH2SO4+2.*WNH3 
      REAL, PARAMETER :: WRAT    =WAMSUL/WSO2
C-------------------------------------------------------------------
C
C     * EBBT, EBFT, EOBT, AND EOFT
C
      POM2OC=1./1.4
      DO 30 IL=IL1,IL2
        FACT                 = ZTMST*GRAV/(DSHJ(IL,ILEV)*PRESSG(IL))
        EBC=EBBTROW(IL)+EBFTROW(IL)
        EOC=POM2OC*(EOBTROW(IL)+EOFTROW(IL))
        IF ( IPAM.EQ.0 ) THEN
          XEMIS(IL,ILEV,IBCO)=XEMIS(IL,ILEV,IBCO)+0.8*EBC*FACT
          XEMIS(IL,ILEV,IBCY)=XEMIS(IL,ILEV,IBCY)+0.2*EBC*FACT
          XEMIS(IL,ILEV,IOCO)=XEMIS(IL,ILEV,IOCO)+0.5*EOC*FACT
          XEMIS(IL,ILEV,IOCY)=XEMIS(IL,ILEV,IOCY)+0.5*EOC*FACT
        ELSE
          DOCEM1(IL,ILEV)=DOCEM1(IL,ILEV)+EOBTROW(IL)*FACT/ZTMST
          DOCEM2(IL,ILEV)=DOCEM2(IL,ILEV)+EOFTROW(IL)*FACT/ZTMST
          DBCEM1(IL,ILEV)=DBCEM1(IL,ILEV)+EBBTROW(IL)*FACT/ZTMST
          DBCEM2(IL,ILEV)=DBCEM2(IL,ILEV)+EBFTROW(IL)*FACT/ZTMST
        ENDIF
   30 CONTINUE      
C
      JT=ISO2
C
C     * ESDT, ESST, ESOT AND ESRT
C
      DO 50 IL=IL1,IL2
        FACT    = ZTMST*GRAV/(DSHJ(IL,ILEV)*PRESSG(IL))
        ESDT=ESDTROW(IL)*(1.-PSUEF2)*32.064/64.059
        ESO2=ESSTROW(IL)*(1.-PSUEF2)*32.064/64.059
        ESO2=ESO2+( ESOTROW(IL)+ESRTROW(IL) )*(1.-PSUEF2)*32.064/64.059
        XEMIS(IL,ILEV,JT)=XEMIS(IL,ILEV,JT)+(ESDT+ESO2)*FACT
        IF ( IPAM.EQ.1 ) THEN
          DSUEM2(IL,ILEV)=DSUEM2(IL,ILEV)
     1          +(ESDTROW(IL)+ESOTROW(IL)+ESRTROW(IL)+ESSTROW(IL))*FACT
     2                    *WRAT*PSUEF2/ZTMST
        ENDIF
   50 CONTINUE
C
C     * DETERMINE LAYER CORRESPONDING TO POWER PLANT (ESPTROW) AND
C     * INDUSTRIAL (ESITROW) SO2 EMISSIONS (100-300m).
C
      DO 100 IL = IL1, IL2
        LSO2(IL) = ILEV
        HSO2(IL) = ILEV
  100 CONTINUE
C
      DO 105 L = ILEV, 1, -1
      DO 105 IL = IL1, IL2
          IF(ZF(IL,L) .LE. 100.)   LSO2(IL)=L
          IF(ZF(IL,L) .LE. 300.)   HSO2(IL)=L
  105 CONTINUE
C
      DO 110 L=1, ILEV
      DO 110 IL = IL1, IL2
        ESPI = ESPTROW(IL)+ESITROW(IL)
        FACT                 = ZTMST*GRAV/(DSHJ(IL,L)*PRESSG(IL))
        IF(L.EQ.HSO2(IL) .AND. L.EQ.LSO2(IL)) THEN
         WGT=1.
        ELSE IF(L.EQ.HSO2(IL)) THEN
         WGT=(300.-ZF(IL,L))/200.
        ELSE IF(L.GT.HSO2(IL).AND.L.LT.LSO2(IL)) THEN
         WGT=(ZF(IL,L-1)-ZF(IL,L))/200.
        ELSE IF(L.EQ.LSO2(IL)) THEN
         WGT=(ZF(IL,L-1)-100.)/200.
        ELSE
         WGT=0.
        ENDIF
        XEMIS(IL,L,JT)=XEMIS(IL,L,JT)+WGT*ESPI*FACT*32.064/64.059
     1                 *(1.-PSUEF2)
        IF ( IPAM.EQ.1 ) THEN
          DSUEM2(IL,L)=DSUEM2(IL,L)+WGT*ESPI*FACT*WRAT*PSUEF2/ZTMST
        ENDIF
  110 CONTINUE
C
      RETURN
      END
