C234567890123456789012345678901234567890123456789012345678901234567890
C
C     ONE-DIMENSIONAL HETEROGENEOUS CHEMISTRY
C
      SUBROUTINE HET_1D(TSTEP,TD1,YD1,PSC1,PSC2,AR,HNO3S,H2OS,H2SO4S)
      INCLUDE 'basp.h'
      INCLUDE 'din.h'
      INCLUDE 'garr.h'
      REAL TD1(L),YD1(L,NC),PSC1(L),PSC2(L),AR(L)
      REAL HNO3S(L),H2OS(L),H2SO4S(L)
      REAL H2O(L),HNO3(L),H2SO4(L),gH2O(L),gHNO3(L),TM(L)
      IF(NOEB.NE.0)GOTO3
      DO 2 I=1,L
         TM(I)=TD1(L-I+1)
         HNO3(I)=YD1(I,NHNO3)+HNO3S(I)
         HNO3S(I)=0.
         H2O(I)=YD1(I,KH2O)+H2OS(I)
         H2OS(I)=0.
         YD1(I,NNOY)=YD1(I,NNOY)-YD1(I,NHNO3)
         H2SO4(I)=0.5E-9*YD1(I,NC)
!      write(*,'(4f10.0,5(1pe10.2))')zkm(i),zst(i),zst(i+1),TM(i),
!     $HNO3(i),h2o(i),AR(i),h2so4(i)
  2      CONTINUE
!      pause 
      CALL PSC(HNO3,H2O,H2SO4,TM,ZKM,ZST,L,LL,TSTEP,IPSC1,IPSC2,
     $AR,PSC1,PSC2,gH2O,gHNO3,HNO3S,H2OS)
      DO 1 I=1,L
!	write(*,'(2f6.1,7(1pe9.2))')zkm(i),tm(i),ar(i),psc1(i),psc2(i),
!     $gH2O(i),gHNO3(i),hno3s(i),h2os(i)
         YD1(I,NHNO3)=gHNO3(I)
         YD1(I,NNOY)=YD1(I,NNOY)+YD1(I,NHNO3)!gHNO3(I)
         YD1(I,KH2O)=gH2O(I)
   1     CONTINUE
!      pause
   3  RETURN
      END

C234567890123456789012345678901234567890123456789012345678901234567890

      SUBROUTINE PSC(HNO3,H2O,H2SO4,T,ZKM,ZST,L,LL,TSTEP,IPSC1,IPSC2,
     $AER,PSC1,PSC2,gH2O,gHNO3,HNO3s,H2Os)
C
C
C     INPUT:
C
c       HNO3    - total HNO3, cm^-3
c       H2O     - numb. dens., cm^-3
c       AER     - AEROSOL surface area
C       T       - Temperature, K
C
C     OUTPUT:
C
C       HNO3s   - solid HNO3, cm^-3
C       gHNO3   - gaseous HNO3, cm^-3
C       H2Os    - solid H2O, cm^-3
C       gH2O    - gaseous H2O, cm^-3
C       PSC1    - PSC1 surface area
C       PSC2    - PSC2 surface area
C
      PARAMETER(N_B=40,LHET=100)
      REAL Na,ro(2),AMUR(2)  ! 'kBoltz' not used GM 1/00 
      REAL HNO3(L),HNO3s(L),H2O(L),H2Os(L),H2SO4(L),gH2O(L),gHNO3(L),
     $T(L),AER(L),PSC1(L),PSC2(L),ZKM(L),ZST(LL),HNO3I(LHET),
     $RNO(LHET),RNOI(LHET)
      DIMENSION F_SAD(N_B),F_NAT(N_B),F_ICE(N_B),R_A(N_B),DR(N_B)      
      DIMENSION PAR(LHET),RALT(LHET)
C
C     VOLUME OF ICE,SAD AND STS
C
      DIMENSION VICE(LHET),VOL_SA(LHET),VOLA(LHET)
C
C     DENSITY OF SAD AND STS
C
      DIMENSION RHO0(LHET),RHO(LHET)
C
C     SEDIMENTATION PARAMETERS
C
      DIMENSION W(LHET),SR(LHET),SRS(LHET),SRT(LHET)
      REAL*8 AK_CARS
      COMMON/BINSOL/AK_CARS(7,2),Q_CARS(10,2)
      DATA ssatH2O/1./,Na/6.0225E23/,PI/3.141592/
C not used GM 1/00   DATA ssatHNO3/1./, R/8.205E-5/,kBoltz/1.3806e-19/  
      DATA ro/ 1.62, 0.928/, AMUR/117.,18./
C not used GM 1/00   DATA AN_SA/8./, AN_NAT/0.2/, AN_ICE/0.01/
      DATA R0_SA/0.07E-4/, R0_NAT/0.4E-4/, R0_ICE/5.0E-4/
      DATA S_SA/1.7/, S_NAT/1.6/, S_ICE/1.8/
C234567890123456789012345678901234567890123456789012345678901234567890
      DO 35 I=1,N_B
  35     R_A(I)=1.E-4*10**(0.08*(I-24))         !  PARTICLE RADIUS, cm
      DR(1)=R_A(2)-R_A(1)
      DO 41 I=2,N_B-1
  41     DR(I)=0.5*(R_A(I+1)-R_A(I-1))
      DR(N_B)=R_A(N_B)-R_A(N_B-1)
C
C     MOMENTS OF LOGNORMAL DISTRIBUTION
C
      CALL DIS_LN(N_B,R_A,F_SAD,R0_SA, S_SA, RA,  S2_SA, S3_SA )
      CALL DIS_LN(N_B,R_A,F_NAT,R0_NAT,S_NAT,RNAT,S2_NAT,S3_NAT)
      CALL DIS_LN(N_B,R_A,F_ICE,R0_ICE,S_ICE,RICE,S2_ICE,S3_ICE)
      DO 1 I=1,L
         PAR(I)=AER(I)/(4.*PI*S2_SA)          ! Number of SAD particles
!         if(par(i).gt.0)then
!            print *,zkm(i),par(i),aer(i)
!            pause
!         endif
         PSC1(I)=0.                           ! PSC1 surface area, cm(-1)
         PSC2(I)=0.                           ! PSC2 surface area, cm(-1)
   1     CONTINUE
      AVRNAT=0.
      CAVRNAT=0.
      DO 2 I=1,L
         Temp=T(I)
C
C        WATER VAPOR CAN NOT EXCEED ITS SATURATED VALUE
C
         H2Osat=H2O_SAT(Temp)                  ! SATURATED H2O
         H2Oss=ssatH2O*H2Osat                  ! OVERCOOLED VALUE OF H2O
         H2Os(I)=AMAX1(0.,H2O(I)-H2Oss)        ! SOLID H2O (ice)
         gH2O(I)=H2O(I)-H2Os(I)                ! WATER VAPOR
         VICE(I)=AMUR(2)*H2Os(I)/(Na*ro(2))    ! VOLUME OF ICE PARTICLES
C
C        VOLUME OF BINARY SOLUTION H2SO4/H2O (cm(-3)/cm(-3))
C
!         VOL_SA(I)=4.*PI*PAR(I)*S3_SA/3.!+H2SO4s(I)*98.12/(Na*RHO0(I))
C
C        TERNARY SOLUTION H2SO4/HNO3/H2O
C
C                VOLA  - VOLUME OF TERNARY SOLUTION
C                HNO3s - SOLID HNO3
C
         CALL CARS(Temp,H2O(I),HNO3(I),H2SO4(I),VOLA(I),
     *   WS,WN,RHO(I),HNO3s(I),HNO3I(I))
         gHNO3(I)=HNO3(I)-HNO3s(I)-HNO3I(I)    ! GASEOUS HNO3, cm(-3)
         RALT(I)=R0_NAT
         if(VOLA(i).gt.1.e-15)then
           IF(par(i).eq.0.)PAR(I)=10.
           RALT(I)=(3.*VOLA(I)/(4.*PI*PAR(I)))**(1./3.) !*0.595238
           AVRNAT=AVRNAT+RALT(I)
           CAVRNAT=CAVRNAT+1.
!            CALL VILKA(VOLA(I),S_SA,R_A,F_SAD,DR,N_B,PAR(I),RALT(I))
!            print *,zkm(i),par(i),vola(i),vol_sa(i),aer(i)
!            pause
         endif
   2     CONTINUE
C
C     SEDIMENTATION
C
      IF(CAVRNAT.EQ.0.)GOTO8
      AVRNAT=AVRNAT/CAVRNAT
      CALL DIS_LN(N_B,R_A,F_NAT,AVRNAT,S_NAT,RNAT,S2_NAT,S3_NAT)
      VU_NAT=4.*PI*S3_NAT/3.
      VU_ICE=4.*PI*S3_ICE/3. 
      DO 3 JI=1,N_B
         VOL_PAR=4.*PI*R_A(JI)**3.*DR(JI)/3.
C
C        ICE AND SOLID H2O AND HNO3 SEDIMENTATION
C
         VUF_ICE=F_ICE(JI)*VOL_PAR/VU_ICE
         DO 5 I=1,L
            W(I)=WSED(R_A(JI),ZKM(I),ro(2))
            SRS(I)=H2Os(I)*VUF_ICE
            SRT(I)=HNO3I(I)*VUF_ICE
   5        SR(I)=VICE(I)*VUF_ICE
         CALL SED(VICE,SR,L,ZST,LL,w,TSTEP)
         CALL SED(H2Os,SRS,L,ZST,LL,w,TSTEP)
         CALL SED(HNO3I,SRT,L,ZST,LL,w,TSTEP)
C
C        STS AND LIQUID HNO3 SEDIMENTATION
C
         VUF_NAT=F_NAT(JI)*VOL_PAR/VU_NAT
         DO 6 I=1,L
            W(I)=WSED(R_A(JI),ZKM(I),RHO(I))
            SRS(I)=HNO3s(I)*VUF_NAT
   6        SR(I)=VOLA(I)*VUF_NAT
!         CALL SED(VOLA,SR,L,ZST,LL,w,TSTEP)
!         CALL SED(HNO3s,SRS,L,ZST,LL,w,TSTEP)
   3     CONTINUE
C
C     SURFACE AREA 
C
      DO 7 I=1,L
         Temp=T(I)
!         CALL DIS_LN(N_B,R_A,F_NAT,RALT(I),S_NAT,RNAT,S2_NAT,S3_NAT)
CSPS         IF(Temp.LT.205.)THEN
CSPS            AER(I)=0.
CSPS         ELSE IF(Temp.LT.212.)THEN
CSPS            AER(I)=3.*VOLA(I)*S2_SA/S3_SA
CSPS         ENDIF
C
C        PSCI
C
!         if(VOLA(i).gt.1.e-15.and.par(i).gt.0.)then
!            RALT(I)=(3.*VOLA(I)/(4.*PI*PAR(I)))**(1./3.)*0.595238
!         ENDIF
         IF(IPSC1.EQ.0.AND.Temp.LE.220..AND.VOLA(i).gt.1.e-15)THEN
            RALT(I)=(3.*VOLA(I)/(4.*PI*PAR(I)))**(1./3.)*0.5 !*0.595238
            CALL DIS_LN(N_B,R_A,F_NAT,RALT(I),S_NAT,RNAT,S2_NAT,S3_NAT)
            SURFACE=3.*VOLA(I)/RNAT !*S2_NAT/S3_NAT ! cm(2)/cm(3    
            IF(Temp.LE.205.)THEN
               PSC1(I)=SURFACE
!               PSC1(I)=3.*VOLA(I)/RALT(I) 
!               AER(I)=0.
            ELSE
!               AER(I)=AMAX1(AER(I),SURFACE) !3.*VOLA(I)/RALT(I))
               PSC1(I)=0.
            ENDIF
         ENDIF
C
C        PSCII
C
         IF(IPSC2.EQ.0.AND.Temp.LE.188.)
     *   PSC2(I)=3.*VICE(I)*S2_ICE/S3_ICE   ! cm(2)/cm(3)
         HNO3S(I)=HNO3s(I)+HNO3I(I)
   7     CONTINUE
   8  RETURN
      END
C234567890123456789012345678901234567890123456789012345678901234567890
C
      SUBROUTINE H2O_S(H2O,T,H2OSAT)
C
C     H2O CAN'T EXCEED IT'S SATURATION VALUE
C
      H2OSAT=H2O_SAT(T)
      IF(H2O.GT.H2OSAT)H2O=H2OSAT
      RETURN
      END
C234567890123456789012345678901234567890123456789012345678901234567890

      FUNCTION H2O_SAT(T)
      REAL NA,KBOLTZ
      DATA
     &         KBOLTZ / 1.3806e-19 /,   ! Boltzman's const
     &         NA / 6.0225E23 /         ! Avogadro's const, mol^-1
      AMP0=KBOLTZ*T
      IF(T.LT.273.15)THEN
         H2O_SAT=10.**(12.537-2663.5/T)/(AMP0*100.)
      ELSE
         TC=T-273.15
         E=6.1*10.**(7.45*TC/(235.+TC))
         H2O_SAT=217.E-6*E/T/18.*NA
      ENDIF
      RETURN
      END
C234567890123456789012345678901234567890123456789012345678901234567890

      FUNCTION WSED(Rcm,Z,RO)
C
C     FALLING SPEED OF AEROSOL PARTICLES
C     F.KASTEN, J.APPL.MET.,7,944-947,1968.
C
C     INPUT:
C            Rcm  - RADIUS OF PARTICLE, cm
C            Z    - ALTITUDE, km
C            RO   - MASS DENSITY OF AEROSOL PARTICLES, gm/cm(3)
C     OUTPUT:
C            WSED - FALLING SPEED, cm/sec
C      
      PARAMETER(LL=31,NN=8)
      DIMENSION WF(LL,NN),ZF(LL),RF(NN),AMN(NN)

      DATA WF/
     *        4.1150,5.1727,6.5809,8.4837,11.096,14.747,20.005,27.337,
     *        37.362,51.065,69.792,95.383,130.00,176.66,239.37,323.44,
     *        435.86,586.05,782.70,1038.3,1368.8,1793.3,2335.9,3025.5,
     *        3893.3,4988.9,6391.9,8177.8,10492.,13510.,17463.,
     *        1.4283,1.7825,2.2537,2.8901,3.7630,4.9831,6.7369,9.1809,
     *        12.522,17.090,23.332,31.862,43.399,58.953,79.856,107.88,
     *        145.35,195.41,260.96,346.17,456.33,597.82,778.68,1008.5,
     *        1297.8,1663.0,2130.7,2726.0,3497.3,4503.5,5821.0,
     *        4.8096,5.8811,7.3058,9.2287,11.864,15.544,20.814,28.140,
     *        38.159,51.858,70.582,96.165,130.77,177.43,240.13,324.20,
     *        436.61,586.79,783.42,1039.0,1369.5,1794.0,2336.5,3026.1,
     *        3893.9,4989.6,6392.5,8178.4,10492.,13511.,17463.,
     *        2.3182,2.6796,3.1600,3.8083,4.6970,5.9377,7.6918,10.115,
     *        13.440,17.995,24.227,32.743,44.268,59.812,80.706,108.72,
     *        146.19,196.23,261.76,346.95,457.10,598.57,779.42,1009.3,
     *        1298.5,1663.8,2131.4,2726.7,3498.1,4504.3,5821.7,
     *        1.4008,1.5224,1.6778,1.8818,2.1565,2.5363,3.0490,3.7427,
     *        4.7117,6.0547,7.9057,10.441,13.882,18.533,24.790,33.186,
     *        44.418,59.418,79.065,104.61,137.64,180.08,234.32,303.27,
     *        390.04,499.61,639.90,818.50,1049.9,1351.8,1747.0,
     *        1.3188,1.3890,1.4725,1.5736,1.6985,1.8575,2.0299,2.2259,
     *        2.5077,2.9122,3.4878,4.2890,5.3987,6.9179,8.9785,11.757,
     *        15.484,20.463,26.993,35.491,46.488,60.619,78.688,101.66,
     *        130.58,167.10,213.86,273.40,350.54,451.17,582.93,
     *        1.1264,1.1738,1.2279,1.2903,1.3635,1.4510,1.5235,1.5776,
     *        1.6521,1.7559,1.9027,2.1034,2.3936,2.8074,3.3872,4.1873,
     *        5.2775,6.7402,8.6733,11.201,14.482,18.706,24.114,30.993,
     *        39.660,50.615,64.642,82.506,105.65,135.85,175.38,
     *        1.2280,1.2745,1.3264,1.3848,1.4509,1.5269,1.5769,1.5942,
     *        1.6183,1.6515,1.6972,1.7489,1.8237,1.9311,2.0855,2.3077,
     *        2.6247,3.0564,3.6503,4.4497,5.5071,6.8847,8.6620,10.934,
     *        13.808,17.453,22.124,28.080,35.799,45.868,59.050/
      DATA RF/0.003,0.01,0.03,0.1,0.3,1.,3.,10./, 
     *     AMN/1.E-6,1.E-5,1.E-5,1.E-4,1.E-3,1.E-2,1.E-1,1./
      R=Rcm*1.E4
      DO 1 I=1,LL
   1     ZF(I)=2.*REAL(I-1)
      II=1
      DO WHILE(ZF(II)+0.1.LT.Z)
         II=II+1
         IF(II.GT.LL)GOTO2
      ENDDO
      IF(II.EQ.1)II=II+1
      ZXP=ZF(II)-Z
      ZXM=Z-ZF(II-1)
      KK=1
      DO WHILE(RF(KK)+0.0001.LT.R)
         KK=KK+1
         IF(KK.GE.NN)GOTO4
      ENDDO
      IF(KK.EQ.1)KK=KK+1
      YXP=(RF(KK)-R)*AMN(KK-1)
      YXM=(R-RF(KK-1))*AMN(KK)
      ZN=1./((ZF(II)-ZF(II-1))*(RF(KK)-RF(KK-1)))
      WSED=((WF(II-1,KK-1)*ZXP+WF(II,KK-1)*ZXM)*YXP+
     *      (WF(II-1,KK)*ZXP+WF(II,KK)*ZXM)*YXM)*ZN
      GOTO3
   4  YXP=ALOG10(RF(KK)/R)
      YXM=ALOG10(R/RF(KK-1))
      ZN=1./((ZF(II)-ZF(II-1))*(ALOG10(RF(KK)/RF(KK-1))))
      WSED=10.**(((ALOG10(WF(II-1,KK-1)*AMN(KK-1))*ZXP+
     *     ALOG10(WF(II,KK-1)*AMN(KK-1))*ZXM)*YXP+
     *     (ALOG10(WF(II-1,KK)*AMN(KK))*ZXP+
     *     ALOG10(WF(II,KK)*AMN(KK))*ZXM)*YXM)*ZN)
      GOTO3
   2  WSED=0.
   3  WSED=RO*WSED
      RETURN
      END
C234567890123456789012345678901234567890123456789012345678901234567890

      SUBROUTINE SED(Y,S,L,Z,LL,W,DT)
C     
C     SEDIMENTATION OF ATMOSPHERIC PARTICLE Y IN SOLID STATE WITH SPEED W
C
C     INPUT:
C            L       - NUMBER OF ALTITUDE LEVELS
C            Z(LL)   - STRAGGLED ALTITUDE VALUES, CM
C            LL=L+1  - STRAGGLED GRID
C            W       - FALLING SPEED, CM/SEC
C            DT      - TIME STEP, SEC
C            S       - SEDIMENTED CONCENTRATION OF Y, cm(-3)
C
C     OUTPUT:
C            Y(L)    - IMPLICIT CONCENTRATION OF Y, cm(-3)
C
      DIMENSION Y(L),Z(LL),W(L),S(L)   
      
      DO 1 I=1,L
         II=I
         E=W(I)*DT/(Z(I+1)-Z(I))   
   3     IF(E.LE.1.)GOTO2
         Y(II)=Y(II)-S(I)
         IF(II.EQ.1)GOTO1
         Y(II-1)=Y(II-1)+S(I)
         E=E-1.
         II=II-1
         GOTO3
   2     HSED=S(I)*E
         Y(II)=Y(II)-HSED
         IF(II.EQ.1)GOTO1
         Y(II-1)=Y(II-1)+HSED
   1     CONTINUE
      RETURN
      END
C234567890123456789012345678901234567890123456789012345678901234567890
      SUBROUTINE CARS(T,H2O,HNO3,H2SO4,V,WS,WN,RHO,HNO3A,HNO3I)
C
C     CARSLAW ET AL., GRL,22,1877-1880,1995.
C
C     INPUT:
C    
C     T     - TEMPERATURE (K)
C     H2O   - TOTAL H2O,   CM(-3)
C     HNO3  - TOTAL HNO3,  CM(-3)
C     H2SO4 - TOTAL H2SO4, CM(-3)
C
C     OUTPUT:
C
C     V     - VOLUME OF AEROSOL PER UNIT VOLUME OF AIR, 
C     WS    - WEIGHT FRACTION OF H2SO4 IN THE AEROSOL
C     WN    - WEIGHT FRACTION OF HNO3 IN THE AEROSOL
C     RHO   - TERNARY SOLUTION MASS DENSITY, KG/M(3)
C     PN    - EQUILIBRIUM VAPOUR PRESSURE OF HNO3 OVER THE DROPLET, ATM
C     HNO3A - SOLID HNO3, CM(-3)
C
      EXTERNAL AM_BIN,RHO_S_CAR,RHO_N_CAR
      REAL*8 AK_CARS,AM_BIN,RHO_S_CAR,RHO_N_CAR
      REAL*8 AMNB,AMSB,ANB,ASB,TT,ANBT,ASBT,AMSB2,AMNB2,AMNB_AMSB,
     $A2B3,A2C,A,B,C,PN0,PN,PHI,PI,AMS,AMN,RHOS,RHON,DISCR
      COMMON/BINSOL/AK_CARS(7,2),Q_CARS(10,2)
      DATA R/8.205E-5/,           ! M(3)*ATM*MOL(-1)*K(-1)
     *     AVOG/6.0225E23/,       ! MOL(-1)
     *     AKBOLZ/ 1.3806E-19 /   ! Boltzman's const
C
C     WATER PARTIAL PRESSURE
C
      PW=H2O*R/AVOG*T*1.E6        ! ATM
      PW_MB=PW*1.013E3            ! mb
C
C     HNO3 PARTIAL PRESSURE
C
      PN0=HNO3*R/AVOG*T*1.E6      ! ATM
C
C     H2SO4 MOL/M(3)
C
      PH=H2SO4*1.E6/AVOG
C
C     VALID ONLY FOR 2.E-5MB < PW < 2.E-3MB
C                    MAX(TICE-3K,185K) < T < 215K
C
      IF(PW_MB.GT.2.E-5.AND.PW_MB.LT.2.E-3.AND.T.LE.240.)GOTO3
 8    V=0.
      HNO3A=0.
      HNO3I=0.
      GOTO2
 3    TICE=2663.5/(12.537-ALOG10(AKBOLZ*T*H2O*100.))
      TMIN=AMAX1(TICE-3.,187.)
      IF(T.LE.220.)GOTO9
      HNO3I=0.                                        ! NO FROZEN HNO3
 10   AMS=AM_BIN(T,PW,AK_CARS(1,2)) ! (KG H2O)(-1)
      AMN=0.
      PN=PN0
      GOTO7      
 9    IF(T.GE.TMIN)GOTO1
 6    AM = -2.7836 - 0.00088*T
      B = 36.1047 - 11397./T + 0.009179*T
      AMP0=AKBOLZ*T
      AMP=AMP0*H2O
      pH2O = AMP/1013.25
      pHNO3s=AM*(ALOG10(pH2O)+2.8808)+B
      HNO3s=AMIN1(HNO3,10.**(pHNO3S)/AMP0*1013.25)    ! SATURATED HNO3
c      print *,' T=',T,' HNO3=',hno3,' HNO3s=',hno3s
c      pause
      IF(T.LT.TMIN)THEN
         HNO3A=0.                                     ! ALL HNO3 IS FROZEN
         HNO3I = HNO3-hno3s                                 ! SOLID HNO3
         V=117.*HNO3I/(AVOG*0.928)                      ! VOLUME OF SOLID HNO3
      ELSE
         HNO3I=0.                                     ! NO FROZEN HNO3
         HNO3A = HNO3-HNO3s                           ! SOLID HNO3
         RHO=1.62                                     ! DENSITY OF SOLID HNO3
         V=117.*HNO3A/(AVOG*RHO)                      ! VOLUME OF SOLID HNO3
      ENDIF
      GOTO2
 1    HNO3I=0.
C
C     CONCENTRATION OF H2SO4 IN THE BINARY SOLUTION WITH H2O
C
      AMSB=AM_BIN(T,PW,AK_CARS(1,2)) ! (KG H2O)(-1)
      IF(PH)6,6,5
  5   CONTINUE
C
C     H2SO4 PARTIAL PRESSURE
C
      TT=R*T*PH
C
C     HENRY LAW CONSTANTS
C
      HNNB=HENRIC(T,PW,Q_CARS(1,1))  ! MOL*KG(-1)*ATM(-1)
      HNSB=HENRIC(T,PW,Q_CARS(1,2))  ! MOL*KG(-1)*ATM(-1)
C
C     CONCENTRATION OF HNO3 IN THE BINARY SOLUTION WITH H2O
C
      AMNB=AM_BIN(T,PW,AK_CARS(1,1)) ! (KG H2O)(-1)
      IF(AMNB.LE.0.)GOTO10
      ANB=HNNB*AMNB
      ASB=HNSB*AMSB
      ANBT=TT*ANB
      ASBT=TT*ASB
      AMSB2=AMSB*AMSB
      AMNB2=AMNB*AMNB
      AMNB_AMSB=AMNB-AMSB
      A=((ANBT-ASBT)*AMNB-2.*AMNB2*AMSB+AMNB*AMSB2+
     *(ANB-ASB)*AMSB*PN0)/(AMNB2-AMNB*AMSB)
      B=AMSB*(-2.*ANBT+ASBT+AMSB*AMNB-HNNB*AMSB*PN0)/AMNB_AMSB
      C=ANBT*AMSB2/AMNB_AMSB
      A2B3=A*A-3.*B
      A2C=-2.*A*A*A+9.*A*B-27.*C
      DISCR=4.*A2B3*A2B3*A2B3-A2C*A2C
      IF(DISCR.LE.0.D0)GOTO10
      PHI=DATAN(DSQRT(DISCR)/A2C)
      PI=2.*DASIN(1.D0)
      IF(PHI.LT.0.D0)PHI=PHI+PI
      AMS=-(A+2.*DSQRT(A2B3)*DCOS((PI+PHI)/3.))/3.
      IF(AMS.LE.0.)GOTO6
      AMN=DMAX1(0.D0,AMNB*(1.-AMS/AMSB))
      PN=AMN*(AMN+AMS)/(HNNB*AMN+HNSB*AMS)
 7    ZNAM=1./(1.+AMS*0.09812+AMN*0.06303)
      WS=AMS*0.09812*ZNAM
      WN=AMN*0.06303*ZNAM
      RHOS=RHO_S_CAR(T,AMS)
      RHON=RHO_N_CAR(T,AMN)
      RHO=1.E-3/((1./RHOS)*AMS/(AMN+AMS)+(1./RHON)*AMN/(AMN+AMS))
      V=PH*0.09812E-3/(WS*RHO)
      HNO3A=DMAX1(0.D0,(PN0-PN)*AVOG/(R*T*1.D6))
 2    RETURN
      END

      REAL*8 FUNCTION AM_BIN(T,P,AK)
      EXTERNAL X_CARS
      REAL*8 AK(7),X,X_CARS
      X=X_CARS(T,P,AK)
      AM_BIN=DMAX1(55.51D0*X/(1.D0-X),0.D0)
      RETURN
      END

      REAL*8 FUNCTION X_CARS(T,P,AK)
      REAL*8 AK,ASUM1,ASUM2
      DIMENSION AK(7)
      ASUM1=AK(1)+AK(2)/T
      ASUM2=2.*(AK(3)+AK(4)/T)
      X_CARS=(-ASUM1-DSQRT(ASUM1*ASUM1-2.*ASUM2*
     *  (AK(5)+AK(6)/T+AK(7)*ALOG(T)-ALOG(P))))/
     *  ASUM2
      RETURN
      END

      FUNCTION HENRIC(T,P,Q)
      DIMENSION Q(10)
      PR=ALOG(P)+18.4
      TR=1.E4*(1./T-1./230.)
      TR2=TR*TR
      PR2=PR*PR
      AHEN=Q(1)+Q(2)*TR2+(Q(3)+Q(4)*TR+Q(5)*TR2+Q(6)*TR2*TR)*PR+
     *     (Q(7)+Q(8)*TR+Q(9)*TR2)*PR2+Q(10)*TR*PR2*PR
      HENRIC=EXP(AHEN)
      RETURN
      END

      REAL*8 FUNCTION RHO_S_CAR(T,AM)
C
C     DENSITY OF H2SO4 IN A BINARY SOLUTION, KG/M(3)
C      
      REAL*8 AM,T2
      T2=T*T
C      print *,am
      RHO_S_CAR=1000.+(123.64-5.6E-4*T2)*AM-(29.54-1.81E-4*T2)*AM**1.5
     +         +(2.343-1.487E-3*T-1.324E-5*T2)*AM*AM
      RETURN
      END

      REAL*8 FUNCTION RHO_N_CAR(T,AM)
C
C     DENSITY OF HNO3 IN A BINARY SOLUTION, KG/M(3)
C      
      REAL*8 AM,T2
      T2=T*T
      RHO_N_CAR=1000.+(85.11-5.04E-4*T2)*AM-(18.96-1.427E-4*T2)*AM**1.5+
     +     (1.458-1.198E-3*T-9.703E-6*T2)*AM*AM
      RETURN
      END

C234567890123456789012345678901234567890123456789012345678901234567890
      
      FUNCTION FLOG(R,S,R0)
      PI=3.141592
      SL2=ALOG(S)**2.
      R0R=1./R0
      FLOG=1./(SQRT(2.*PI*SL2)*R)*
     *EXP(-(ALOG(R*R0R))**2./(2.*SL2))
      RETURN
      END
C234567890123456789012345678901234567890123456789012345678901234567890

      BLOCK DATA SPS
      REAL*8 AK_CARS
      COMMON/BINSOL/AK_CARS(7,2),Q_CARS(10,2)
C
C     COEFFICIENTS FOR THE CALCULATION OF THE CONCENTRATION OF A
C     BINARY SOLUTION OF HNO3 AND H2SO4
C
      DATA  AK_CARS/
     *-3.9136D1,6.3584D3,8.329D1,-1.7650D4,1.9853D2,-1.1948D4,-2.8469D1,
     *-2.1661D1,2.7242D3,5.181D1,-1.5732D4,4.7004D1,-6.9690D3,-4.6183D0/
C
C     COEFFICIENTS FOR THE CALCULATION OF HENRY'S LAW CONSTANT
C     OF HNO3 AND H2SO4; CARSLAW ET AL.,GRL,22,1877-1880,1995
C
      DATA Q_CARS/14.5734,0.0615994,-1.14895,0.691693,-0.098863,
     *0.0051579,0.123472,-0.115574,0.0110113,0.0097914,
     *14.47,0.0638795,-3.29597,1.778224,-0.223244,0.0086486,
     *0.536695,-0.335164,0.0265153,0.0157550/
      END
C234567890123456789012345678901234567890123456789012345678901234567890

      SUBROUTINE LOGN(R,N,S,R0,F)
      DIMENSION R(N),F(N)
      PI=3.141592
      SL2=ALOG(S)**2.
      R0R=1./R0
      DO 1 I=1,N
   1     F(I)=1./(SQRT(2.*PI*SL2)*R(I))*
     $   EXP(-(ALOG(R(I)*R0R))**2./(2.*SL2))
      RETURN
      END
C234567890123456789012345678901234567890123456789012345678901234567890

      SUBROUTINE MOM(R,N,F,S1,S2,S3)
      DIMENSION R(N),F(N)
      S1=0.
      S2=0.
      S3=0.
      DO 1 I=1,N
         IF(I.EQ.1)THEN
            DR=R(2)-R(1)
         ELSE IF(I.EQ.N)THEN
            DR=R(N)-R(N-1)
         ELSE
            DR=0.5*(R(I+1)-R(I-1))
         ENDIF
         R2=R(I)*R(I)
         S1=S1+F(I)*R(I)*DR
         S2=S2+F(I)*R2*DR
         S3=S3+F(I)*R2*R(I)*DR
   1     CONTINUE
      RETURN
      END
C234567890123456789012345678901234567890123456789012345678901234567890

      SUBROUTINE DIS_LN(N,R,F,R0,S,S1,S2,S3)
      DIMENSION R(N),F(N)
      CALL LOGN(R,N,S,R0,F)
      CALL MOM(R,N,F,S1,S2,S3)
      RETURN
      END
