Skip to content
Snippets Groups Projects
ares.f 51 KiB
Newer Older
  • Learn to ignore specific revisions
  • !MNH_LIC Copyright 1994-2014 CNRS, Meteo-France and Universite Paul Sabatier
    
    !MNH_LIC This is part of the Meso-NH software governed by the CeCILL-C licence
    
    !MNH_LIC version 1. See LICENSE, CeCILL-C_V1-en.txt and CeCILL-C_V1-fr.txt  
    
    !MNH_LIC for details. version 1.
    
    5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300 301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318 319 320 321 322 323 324 325 326 327 328 329 330 331 332 333 334 335 336 337 338 339 340 341 342 343 344 345 346 347 348 349 350 351 352 353 354 355 356 357 358 359 360 361 362 363 364 365 366 367 368 369 370 371 372 373 374 375 376 377 378 379 380 381 382 383 384 385 386 387 388 389 390 391 392 393 394 395 396 397 398 399 400 401 402 403 404 405 406 407 408 409 410 411 412 413 414 415 416 417 418 419 420 421 422 423 424 425 426 427 428 429 430 431 432 433 434 435 436 437 438 439 440 441 442 443 444 445 446 447 448 449 450 451 452 453 454 455 456 457 458 459 460 461 462 463 464 465 466 467 468 469 470 471 472 473 474 475 476 477 478 479 480 481 482 483 484 485 486 487 488 489 490 491 492 493 494 495 496 497 498 499 500 501 502 503 504 505 506 507 508 509 510 511 512 513 514 515 516 517 518 519 520 521 522 523 524 525 526 527 528 529 530 531 532 533 534 535 536 537 538 539 540 541 542 543 544 545 546 547 548 549 550 551 552 553 554 555 556 557 558 559 560 561 562 563 564 565 566 567 568 569 570 571 572 573 574 575 576 577 578 579 580 581 582 583 584 585 586 587 588 589 590 591 592 593 594 595 596 597 598 599 600 601 602 603 604 605 606 607 608 609 610 611 612 613 614 615 616 617 618 619 620 621 622 623 624 625 626 627 628 629 630 631 632 633 634 635 636 637 638 639 640 641 642 643 644 645 646 647 648 649 650 651 652 653 654 655 656 657 658 659 660 661 662 663 664 665 666 667 668 669 670 671 672 673 674 675 676 677 678 679 680 681 682 683 684 685 686 687 688 689 690 691 692 693 694 695 696 697 698 699 700 701 702 703 704 705 706 707 708 709 710 711 712 713 714 715 716 717 718 719 720 721 722 723 724 725 726 727 728 729 730 731 732 733 734 735 736 737 738 739 740 741 742 743 744 745 746 747 748 749 750 751 752 753 754 755 756 757 758 759 760 761 762 763 764 765 766 767 768 769 770 771 772 773 774 775 776 777 778 779 780 781 782 783 784 785 786 787 788 789 790 791 792 793 794 795 796 797 798 799 800 801 802 803 804 805 806 807 808 809 810 811 812 813 814 815 816 817 818 819 820 821 822 823 824 825 826 827 828 829 830 831 832 833 834 835 836 837 838 839 840 841 842 843 844 845 846 847 848 849 850 851 852 853 854 855 856 857 858 859 860 861 862 863 864 865 866 867 868 869 870 871 872 873 874 875 876 877 878 879 880 881 882 883 884 885 886 887 888 889 890 891 892 893 894 895 896 897 898 899 900 901 902 903 904 905 906 907 908 909 910 911 912 913 914 915 916 917 918 919 920 921 922 923 924 925 926 927 928 929 930 931 932 933 934 935 936 937 938 939 940 941 942 943 944 945 946 947 948 949 950 951 952 953 954 955 956 957 958 959 960 961 962 963 964 965 966 967 968 969 970 971 972 973 974 975 976 977 978 979 980 981 982 983 984 985 986 987 988 989 990 991 992 993 994 995 996 997 998 999 1000
    !-----------------------------------------------------------------
    !--------------- special set of characters for RCS information
    !-----------------------------------------------------------------
    ! $Source$ $Revision$
    ! MASDEV4_7 aerosol 2006/05/18 13:07:25
    !-----------------------------------------------------------------
    c/////////////////////////////////////////////////////////////////////////////
                
    C    Calculate the aerosol chemical speciation and water content.
    
          SUBROUTINE RPMARES ( SO4, HNO3, NO3, NH3, NH4, RH, TEMP,
         &                     ASO4, ANO3, AH2O, ANH4, GNH3, GNO3) 
    
    C-----------------------------------------------------------------------
    C
    C Description:
    C
    C   ARES calculates the chemical composition of a sulfate/nitrate/
    C   ammonium/water aerosol based on equilibrium thermodynamics.
    C
    C   This code considers two regimes depending upon the molar ratio 
    C   of ammonium to sulfate. 
    C
    C   For values of this ratio less than 2,the code solves a cubic for 
    C   hydrogen ion molality, HPLUS,  and if enough ammonium and liquid
    C   water are present calculates the dissolved nitric acid. For molal
    C   ionic strengths greater than 50, nitrate is assumed not to be present. 
    C   
    C   For values of the molar ratio of 2 or greater, all sulfate is assumed
    C   to be ammonium sulfate and a calculation is made for the presence of
    C   ammonium nitrate.
    C
    C   The Pitzer multicomponent approach is used in subroutine ACTCOF to
    C   obtain the activity coefficients. Abandoned -7/30/97 FSB 
    
    c   The Bromley method of calculating the activity coefficients is s used
    c    in this version
    
    c   The calculation of liquid water
    C   is done in subroutine water. Details for both calculations are given
    C   in the respective subroutines.
    C
    C   Based upon MARS due to 
    C   P. Saxena, A.B. Hudischewskyj, C. Seigneur, and J.H. Seinfeld, 
    C   Atmos. Environ., vol. 20, Number 7, Pages 1471-1483, 1986.
    C
    C   and SCAPE due to 
    C   Kim, Seinfeld, and Saxeena, Aerosol Ceience and Technology,
    C   Vol 19, number 2, pages 157-181 and pages 182-198, 1993.
    C
    C NOTE: All concentrations supplied to this subroutine are TOTAL
    C       over gas and aerosol phases
    C
    C Parameters:
    C 
    C  SO4   : Total sulfate in MICROGRAMS/M**3 as sulfate (IN)
    C  HNO3  : Nitric Acid in MICROGRAMS/M**3 as nitric acid (IN)
    C  NO3   : Total nitrate in MICROGRAMS/M**3 as nitric acid (IN)
    C  NH3   : Total ammonia in MICROGRAMS/M**3 as ammonia (IN)
    C  NH4   : Ammonium in MICROGRAMS/M**3 as ammonium (IN)
    C  RH    : Fractional relative humidity (IN)
    C  TEMP  : Temperature in Kelvin (IN)
    C  GNO3  : Gas phase nitric acid in MICROGRAMS/M**3 (OUT)
    C  GNH3  : Gas phase ammonia in MICROGRAMS/M**3 (OUT)
    C  ASO4  : Aerosol phase sulfate in MICROGRAMS/M**3 (OUT) 
    C  ANO3  : Aerosol phase nitrate in MICROGRAMS/M**3 (OUT)
    C  ANH4  : Aerosol phase ammonium in MICROGRAMS/M**3 (OUT)
    C  AH2O  : Aerosol phase water in MICROGRAMS/M**3 (OUT)
    C  NITR  : Number of iterations for obtaining activity coefficients  (OUT) 
    C  NR    : Number of real roots to the cubic in the low ammonia case (OUT)
    C 
    C Revision History:
    C      Who       When        Detailed description of changes
    C   ---------   --------  -------------------------------------------
    C   S.Roselle   11/10/87  Received the first version of the MARS code
    C   S.Roselle   12/30/87  Restructured code
    C   S.Roselle   2/12/88   Made correction to compute liquid-phase 
    C                         concentration of H2O2.  
    C   S.Roselle   5/26/88   Made correction as advised by SAI, for 
    C                         computing H+ concentration.
    C   S.Roselle   3/1/89    Modified to operate with EM2
    C   S.Roselle   5/19/89   Changed the maximum ionic strength from 
    C                         100 to 20, for numerical stability.
    C   F.Binkowski 3/3/91    Incorporate new method for ammonia rich case
    C                         using equations for nitrate budget.
    C   F.Binkowski 6/18/91   New ammonia poor case which 
    C                         omits letovicite.
    C   F.Binkowski 7/25/91   Rearranged entire code, restructured
    C                         ammonia poor case.
    C   F.Binkowski 9/9/91    Reconciled all cases of ASO4 to be output
    C                         as SO4--
    C   F.Binkowski 12/6/91   Changed the ammonia defficient case so that 
    C                         there is only neutralized sulfate (ammonium
    C                         sulfate) and sulfuric acid.
    C   F.Binkowski 3/5/92    Set RH bound on AWAS to 37 % to be in agreement 
    C                          with the Cohen et al. (1987)  maximum molality
    C                          of 36.2 in Table III.( J. Phys Chem (91) page
    C                          4569, and Table IV p 4587.)
    C   F.Binkowski 3/9/92    Redid logic for ammonia defficient case to remove
    C                         possibility for denomenator becoming zero; 
    C                         this involved solving for HPLUS first.
    C                         Note that for a relative humidity
    C                          less than 50%, the model assumes that there is no 
    C                          aerosol nitrate.
    C   F.Binkowski 4/17/95   Code renamed  ARES (AeRosol Equilibrium System)  
    C                          Redid logic as follows
    C                         1. Water algorithm now follows Spann & Richardson
    C                         2. Pitzer Multicomponent method used
    C                         3. Multicomponent practical osmotic coefficient 
    C                            use to close iterations.
    C                         4. The model now assumes that for a water
    C                            mass fraction WFRAC less than 50% there is
    C                            no aerosol nitrate.
    C   F.Binkowski 7/20/95   Changed how nitrate is calculated in ammonia poor 
    C                         case, and changed the WFRAC criterion to 40%.
    C                         For ammonium to sulfate ratio less than 1.0 
    C                         all ammonium is aerosol and no nitrate aerosol
    C                         exists.
    C   F.Binkowski 7/21/95   Changed ammonia-ammonium in ammonia poor case to
    C                         allow gas-phase ammonia to exist. 
    C   F.Binkowski 7/26/95   Changed equilibrium constants to values from 
    C                         Kim et al. (1993) 
    C   F.Binkowski 6/27/96   Changed to new water format
    c   F.Binkowski 7/30/97   Changed to Bromley method for multicomponent 
    c                         activity coefficients. The binary activity coefficients 
    c                         are the same as the previous version
    c   F.Binkowski 8/1/97    Chenged minimum sulfate from 0.0 to 1.0e-6 i.e.
    c                         1 picogram per cubic meter
    C
    C-----------------------------------------------------------------------
    *   K. Suhre    6/6/98    
    *
    *   note different values of the following constants,
    *   the latter ones are used in subroutine awater ..
    *
    *           MWSO4   = 96.0576
    *           mwso4   = 96.0636
    *
    *           MWNH4   = 18.03858
    *           mwnh4   = 18.0985
    *
    *           MWNO3   = 62.0049
    *           mwno3   = 62.0649
    *
    * --
    *
    *   increased precision by
    *   changing  TOLER2 = 0.001  to   TOLER2 = 0.00001
    *
    C-----------------------------------------------------------------------
     
          IMPLICIT NONE
    
    C...........INCLUDES and their descriptions
    
    ccc      INCLUDE SUBST_CONST          ! constants
    
    C...........PARAMETERS and their descriptions:
    
          REAL        MWNACL           ! molecular weight for NaCl
          PARAMETER ( MWNACL = 58.44277 )
    
          REAL        MWNO3            ! molecular weight for NO3
          PARAMETER ( MWNO3  = 62.0049 ) 
    
          REAL        MWHNO3           ! molecular weight for HNO3
          PARAMETER ( MWHNO3 = 63.01287 )       
    
          REAL        MWSO4            ! molecular weight for SO4
          PARAMETER ( MWSO4 = 96.0576 )
    
          REAL        MWHSO4           ! molecular weight for HSO4
          PARAMETER ( MWHSO4 = MWSO4 + 1.0080 ) 
    
          REAL        MH2SO4           ! molecular weight for H2SO4
          PARAMETER ( MH2SO4 = 98.07354 ) 
    
          REAL        MWNH3            ! molecular weight for NH3
          PARAMETER ( MWNH3 = 17.03061 ) 
    
          REAL        MWNH4            ! molecular weight for NH4
          PARAMETER ( MWNH4 = 18.03858 )
    
          REAL        MWORG            ! molecular weight for Organic Species
          PARAMETER ( MWORG = 16.0 )
    
          REAL        MWCL             ! molecular weight for Chloride  
          PARAMETER ( MWCL = 35.453 )
    
          REAL        MWAIR            ! molecular weight for AIR
          PARAMETER ( MWAIR = 28.964 )
    
          REAL        MWLCT            ! molecular weight for Letovicite
          PARAMETER ( MWLCT = 3.0 * MWNH4 + 2.0 * MWSO4 + 1.0080 )
    
          REAL        MWAS             ! molecular weight for Ammonium Sulfate
          PARAMETER ( MWAS = 2.0 * MWNH4 + MWSO4 )
    
          REAL        MWABS            ! molecular weight for Ammonium Bisulfate 
          PARAMETER ( MWABS = MWNH4 + MWSO4 + 1.0080 )
    
    C...........ARGUMENTS and their descriptions
    
          REAL        SO4              ! Total sulfate in micrograms / m**3 
    ciamodels3
          REAL        HNO3             ! Total nitric acid in micrograms / m**3
          REAL        NO3              ! Total nitrate in micrograms / m**3
          REAL        NH3              ! Total ammonia in micrograms / m**3
          REAL        NH4              ! Total ammonium in micrograms / m**3
          REAL        RH               ! Fractional relative humidity 
          REAL        TEMP             ! Temperature in Kelvin 
          REAL        ASO4             ! Aerosol sulfate in micrograms / m**3 
          REAL        ANO3             ! Aerosol nitrate in micrograms / m**3
          REAL        AH2O             ! Aerosol liquid water content water in micrograms / m**3
          REAL        ANH4             ! Aerosol ammonium in micrograms / m**3
          REAL        GNO3             ! Gas-phase nitric acid in micrograms / m**3
          REAL        GNH3             ! Gas-phase ammonia in micrograms / m**3
    
    C...........SCRATCH LOCAL VARIABLES and their descriptions:
           
          INTEGER     IRH              ! Index set to percent relative humidity  
          INTEGER     NITR             ! Number of iterations for activity coefficients
          INTEGER     NNN              ! Loop index for iterations 
          INTEGER     NR               ! Number of roots to cubic equation for HPLUS
    
          REAL        A0               ! Coefficients and roots of 
          REAL        A1               ! Coefficients and roots of 
          REAL        A2               ! Coefficients and roots of 
          REAL        AA               ! Coefficients and discriminant for quadratic equation for ammonium nitrate
          REAL        BAL              ! internal variables ( high ammonia case)
          REAL        BB               ! Coefficients and discriminant for quadratic equation for ammonium nitrate
          REAL        BHAT             ! Variables used for ammonia solubility 
          REAL        CC               ! Coefficients and discriminant for quadratic equation for ammonium nitrate
          REAL        CONVT            ! Factor for conversion of units  
          REAL        DD               ! Coefficients and discriminant for quadratic equation for ammonium nitrate
          REAL        DISC             ! Coefficients and discriminant for quadratic equation for ammonium nitrate
          REAL        EROR             ! Relative error used for convergence test  
          REAL        FNH3             ! "Free ammonia concentration", that which exceeds TWOSO4       
          REAL        GAMAAB           ! Activity Coefficient for (NH4+, HSO4-)GAMS( 2,3 )
          REAL        GAMAAN           ! Activity coefficient for (NH4+, NO3-) GAMS( 2,2 )
          REAL        GAMAHAT          ! Variables used for ammonia solubility 
          REAL        GAMANA           ! Activity coefficient for (H+ ,NO3-)   GAMS( 1,2 )
          REAL        GAMAS1           ! Activity coefficient for (2H+, SO4--) GAMS( 1,1 )
          REAL        GAMAS2           ! Activity coefficient for (H+, HSO4-)  GAMS( 1,3 )
          REAL        GAMOLD           ! used for convergence of iteration
          REAL        GASQD            ! internal variables ( high ammonia case)
          REAL        HPLUS            ! Hydrogen ion (low ammonia case) (moles / kg water)
          REAL        K1A              ! Equilibrium constant for ammoniua to ammonium
          REAL        K2SA             ! Equilibrium constant for sulfate-bisulfate (aqueous)
          REAL        K3               ! Dissociation constant for ammonium nitrate 
          REAL        KAN              ! Equilibrium constant for ammonium nitrate (aqueous)
          REAL        KHAT             ! Variables used for ammonia solubility 
          REAL        KNA              ! Equilibrium constant for nitric acid (aqueous)   
          REAL        KPH              ! Henry's Law Constant for ammonia       
          REAL        KW               ! Equilibrium constant for water dissociation             
          REAL        KW2              ! Internal variable using KAN 
          REAL        MAN              ! Nitrate (high ammonia case) (moles / kg water)
          REAL        MAS              ! Sulfate (high ammonia case) (moles / kg water)
          REAL        MHSO4            ! Bisulfate (low ammonia case) (moles / kg water)
          REAL        MNA              ! Nitrate (low ammonia case) (moles / kg water)
          REAL        MNH4             ! Ammonium (moles / kg water)
          REAL        MOLNU            ! Total number of moles of all ions
          REAL        MSO4             ! Sulfate (low ammonia case) (moles / kg water)
          REAL        PHIBAR           ! Practical osmotic coefficient      
          REAL        PHIOLD           ! Previous value of practical osmotic coefficient used for convergence of iteration
          REAL        RATIO            ! Molar ratio of ammonium to sulfate
          REAL        RK2SA            ! Internal variable using K2SA
          REAL        RKNA             ! Internal variables using KNA
          REAL        RKNWET           ! Internal variables using KNA
          REAL        RR1
          REAL        RR2
          REAL        STION            ! Ionic strength
          REAL        T1               ! Internal variables for temperature corrections
          REAL        T2               ! Internal variables for temperature corrections
          REAL        T21              ! Internal variables of convenience (low ammonia case)
          REAL        T221             ! Internal variables of convenience (low ammonia case)
          REAL        T3               ! Internal variables for temperature corrections
          REAL        T4               ! Internal variables for temperature corrections
          REAL        T6               ! Internal variables for temperature corrections
          REAL        TNH4             ! Total ammonia and ammonium in micromoles / meter ** 3
          REAL        TNO3             ! Total nitrate in micromoles / meter ** 3
          REAL        TOLER1           ! Tolerances for convergence test 
          REAL        TOLER2           ! Tolerances for convergence test 
          REAL        TSO4             ! Total sulfate in micromoles / meter ** 3
          REAL        TWOSO4           ! 2.0 * TSO4  (high ammonia case) (moles / kg water)
          REAL        WFRAC            ! Water mass fraction 
          REAL        WH2O             ! Aerosol liquid water content (internally)
                                       ! micrograms / meter **3 on output
                                       ! internally it is 10 ** (-6) kg (water) / meter ** 3
                                       ! the conversion factor (1000 g = 1 kg) is applied 
                                       ! for AH2O output
          REAL        WSQD             ! internal variables ( high ammonia case)
          REAL        XNO3             ! Nitrate aerosol concentration in micromoles / meter ** 3
          REAL        XXQ              ! Variable used in quadratic solution
          REAL        YNH4             ! Ammonium aerosol concentration in micromoles / meter** 3
          REAL        ZH2O             ! Water variable saved in case ionic strength too high.
          REAL        ZSO4             ! Total sulfate molality - mso4 + mhso4 (low ammonia case) (moles / kg water)
    
          REAL        CAT( 2 )         ! Array for cations (1, H+); (2, NH4+) (moles / kg water)
          REAL        AN ( 3 )         ! Array for anions (1, SO4--); (2, NO3-); (3, HSO4-)  (moles / kg water) 
          REAL        CRUTES( 3 )      ! Coefficients and roots of 
          REAL        GAMS( 2, 3 )     ! Array of activity coefficients 
          REAL        MINSO4           ! Minimum value of sulfate laerosol concentration
           PARAMETER( MINSO4 = 1.0E-6 / MWSO4 ) 
          REAL        FLOOR
           PARAMETER( FLOOR = 1.0E-30) ! minimum concentration       
    
    C-----------------------------------------------------------------------
    C  begin body of subroutine RPMARES
                                                                             
    C...convert into micromoles/m**3
     
    ccc      WRITE( 10, * ) 'SO4, NO3, NH3 ', SO4, NO3, NH3
    Ciamodels3 merge NH3/NH4 , HNO3,NO3 here
          TSO4 = MAX( 0.0, SO4 / MWSO4  )
          TNO3 = MAX( 0.0, (NO3 / MWNO3 + HNO3 / MWHNO3) )
          TNH4 = MAX( 0.0, (NH3 / MWNH3 + NH4 / MWNH4)  )
    ccc      WRITE( 10, * ) 'TSO4, TNO3, TNH3, RH ', TSO4, TNO3, TNH3, RH
     
    C...now set humidity index IRH as a percent
    
          IRH = NINT( 100.0 * RH )
    
    C...Check for valid IRH
    
          IRH = MAX(  1, IRH )
          IRH = MIN( 99, IRH )
    ccc      WRITE(10,*)'RH,IRH ',RH,IRH
    
    C...Specify the equilibrium constants at  correct
    C...  temperature.  Also change units from ATM to MICROMOLE/M**3 (for KAN,
    C...  KPH, and K3 )
    C...  Values from Kim et al. (1993) except as noted.
     
          CONVT = 1.0 / ( 0.082 * TEMP ) 
          T6 = 0.082E-9 *  TEMP
          T1 = 298.0 / TEMP
          T2 = ALOG( T1 )
          T3 = T1 - 1.0
          T4 = 1.0 + T2 - T1
          KNA  = 2.511E+06 *  EXP(  29.17 * T3 + 16.83 * T4 ) * T6
          K1A  = 1.805E-05 *  EXP(  -1.50 * T3 + 26.92 * T4 )
          K2SA = 1.015E-02 *  EXP(   8.85 * T3 + 25.14 * T4 )
          KW   = 1.010E-14 *  EXP( -22.52 * T3 + 26.92 * T4 )
          KPH  = 57.639    *  EXP(  13.79 * T3 - 5.39  * T4 ) * T6
    ccc      K3   =  5.746E-17 * EXP( -74.38 * T3 + 6.12  * T4 ) * T6 * T6  
          KHAT =  KPH * K1A / KW  
          KAN  =  KNA * KHAT  
    
    C...Compute temperature dependent equilibrium constant for NH4NO3
    C...  ( from Mozurkewich, 1993)
    
          K3 = EXP( 118.87  - 24084.0 / TEMP -  6.025  * ALOG( TEMP ) )
    
    C...Convert to (micromoles/m**3) **2
    
          K3 = K3 * CONVT * CONVT
    
          WH2O   = 0.0
          STION  = 0.0
          AH2O   = 0.0
          MAS    = 0.0
          MAN    = 0.0
          HPLUS  = 0.0
          TOLER1 = 0.00001
          TOLER2 = 0.00001
          NITR   = 0
          NR     = 0
          RATIO  = 0.0
          GAMAAN = 1.0
          GAMOLD = 1.0
    
    C...set the ratio according to the amount of sulfate and nitrate
    
          IF ( TSO4 .GT. MINSO4 ) THEN
            RATIO = TNH4 / TSO4
    
    C...If there is no sulfate and no nitrate, there can be no ammonium
    C...  under the current paradigm. Organics are ignored in this version.
    
          ELSE 
          
           IF ( TNO3 .EQ. 0.0 ) THEN
    
    C *** If there is very little sulfate and no nitrate set concentrations
    c      to a very small value and return.        
              ASO4 = MAX(FLOOR, ASO4)
              ANO3 = MAX(FLOOR, ANO3 )          
              WH2O = 0.0
              AH2O = 0.0
              GNH3 = MAX(FLOOR,GNH3)
              GNO3 = MAX(FLOOR,GNO3)
              RETURN
           END IF
           
    C...For the case of no sulfate and nonzero nitrate, set ratio to 5
    C...  to send the code to the high ammonia case
    
            RATIO = 5.0
          END IF 
    
     
    C....................................
    C......... High Ammonia Case ........
    C....................................
     
          IF ( RATIO .GT. 2.0 ) THEN
     
            GAMAAN = 0.1
    
    C...Set up twice the sulfate for future use.
    
            TWOSO4 = 2.0 * TSO4
            XNO3 = 0.0            
            YNH4 = TWOSO4
    
    C...Treat different regimes of relative humidity 
    
    C...ZSR relationship is used to set water levels. Units are
    C...  10**(-6) kg water/ (cubic meter of air)
    C...  start with ammomium sulfate solution without nitrate
    
          CALL awater(IRH,TSO4,YNH4,TNO3,AH2O) !**** note TNO3
            WH2O = 1.0E-3 * AH2O  
            ASO4 = TSO4 * MWSO4
            ANO3 = 0.0
            ANH4 = YNH4 * MWNH4
            WFRAC = AH2O / ( ASO4 + ANH4 +  AH2O )
    ccc        IF ( WFRAC .EQ. 0.0 )  RETURN   ! No water       
            IF ( WFRAC .LT. 0.2 ) THEN
     
    C..."dry" ammonium sulfate and ammonium nitrate
    C...  compute free ammonia
    
              FNH3 = TNH4 - TWOSO4
              CC = TNO3 * FNH3 - K3
    
    C...check for not enough to support aerosol      
    
              IF ( CC .LE. 0.0 ) THEN
                XNO3 = 0.0
              ELSE
                AA = 1.0
                BB = -( TNO3 + FNH3 ) 
                DISC = BB * BB - 4.0 * CC
    
    C...Check for complex roots of the quadratic
    C...  set nitrate to zero and RETURN if complex roots are found
    
                IF ( DISC .LT. 0.0 ) THEN
                  XNO3 = 0.0
                  AH2O = 1000.0 * WH2O
                  YNH4 = TWOSO4
                  GNO3 = TNO3 * MWHNO3
                  GNH3 = ( TNH4 - YNH4 ) * MWNH3
                  ASO4 = TSO4 * MWSO4
                  ANO3 = 0.0
                  ANH4 = YNH4 * MWNH4
                  RETURN
                END IF
    
    C...to get here, BB .lt. 0.0, CC .gt. 0.0 always      
    
                DD = SQRT( DISC )
                XXQ = -0.5 * ( BB + SIGN ( 1.0, BB ) * DD )
    
    C...Since both roots are positive, select smaller root.      
    
                XNO3 = MIN( XXQ / AA, CC / XXQ )
              
              END IF
              AH2O = 1000.0 * WH2O
              YNH4 = 2.0 * TSO4 + XNO3
              GNO3 = ( TNO3 - XNO3 ) * MWHNO3
              GNH3 = ( TNH4 - YNH4 ) * MWNH3
              ASO4 = TSO4 * MWSO4
              ANO3 = XNO3 * MWNO3
              ANH4 = YNH4 * MWNH4
              RETURN
    
            END IF
    
    C...liquid phase containing completely neutralized sulfate and
    C...  some nitrate.  Solve for composition and quantity.
     
            MAS = TSO4 / WH2O
            MAN = 0.0
            XNO3 = 0.0
            YNH4 = TWOSO4
            PHIOLD = 1.0
    
    C...Start loop for iteration
     
    C...The assumption here is that all sulfate is ammonium sulfate,
    C...  and is supersaturated at lower relative humidities.
     
            DO 1501 NNN = 1, 150 
              NITR = NNN
              GASQD = GAMAAN * GAMAAN
              WSQD = WH2O * WH2O
              KW2 = KAN * WSQD / GASQD
              AA = 1.0 - KW2
              BB = TWOSO4 + KW2 * ( TNO3 + TNH4 - TWOSO4 )
              CC = -KW2 * TNO3 * ( TNH4 - TWOSO4 )
    
    C...This is a quadratic for XNO3 [MICROMOLES / M**3] of nitrate in solution.
    
              DISC = BB * BB - 4.0 * AA * CC
    
    C...Check for complex roots, if so set nitrate to zero and RETURN 
    
              IF ( DISC .LT. 0.0 ) THEN
                XNO3 = 0.0
                AH2O = 1000.0 * WH2O
                YNH4 = TWOSO4
                GNO3 = TNO3 * MWHNO3
                GNH3 = ( TNH4 - YNH4 ) * MWNH3
                ASO4 = TSO4 * MWSO4
                ANO3 = 0.0
                ANH4 = YNH4 * MWNH4
    ccc            WRITE( 10, * ) ' COMPLEX ROOTS '
                RETURN
              END IF
    
              DD = SQRT( DISC )
              XXQ = -0.5 * ( BB + SIGN ( 1.0, BB ) * DD )
              RR1 = XXQ / AA
              RR2 = CC / XXQ
    
    C...choose minimum positve root         
    
              IF ( ( RR1 * RR2 ) .LT. 0.0 ) THEN
                XNO3 = MAX( RR1, RR2 )
              ELSE 
                XNO3 = MIN( RR1, RR2 )
              END IF
    
              XNO3 = MIN( XNO3, TNO3 )
    
    C...This version assumes no solid sulfate forms (supersaturated ) 
    C...  Now update water
    
              CALL AWATER ( IRH, TSO4, YNH4, XNO3, AH2O )
    
    C...ZSR relationship is used to set water levels. Units are
    C...  10**(-6) kg water/ (cubic meter of air)
    C...  The conversion from micromoles to moles is done by the units of WH2O.
    
              WH2O = 1.0E-3 * AH2O
    
    C...Ionic balance determines the ammonium in solution.
    
              MAN = XNO3 / WH2O
              MAS = TSO4 / WH2O
              MNH4 = 2.0 * MAS + MAN
              YNH4 = MNH4 * WH2O
    
    C...MAS, MAN and MNH4 are the aqueous concentrations of sulfate, nitrate,
    C...  and ammonium in molal units (moles/(kg water) ).
    
              STION = 3.0 * MAS + MAN
              CAT( 1 ) = 0.0
              CAT( 2 ) = MNH4 
              AN ( 1 ) = MAS
              AN ( 2 ) = MAN
              AN ( 3 ) = 0.0
              CALL ACTCOF ( CAT, AN, GAMS, MOLNU, PHIBAR )
              GAMAAN = GAMS( 2, 2 )
    
    C...Use GAMAAN for convergence control
    
              EROR = ABS( GAMOLD - GAMAAN ) / GAMOLD
              GAMOLD = GAMAAN
    
    C...Check to see if we have a solution
    
              IF ( EROR .LE. TOLER1 ) THEN 
    ccc            WRITE( 11, * ) RH, STION, GAMS( 1, 1 ),GAMS( 1, 2 ), GAMS( 1, 3 ),
    ccc     &      GAMS( 2, 1 ), GAMS( 2, 2 ), GAMS( 2, 3 ), PHIBAR
    
                ASO4 = TSO4 * MWSO4
                ANO3 = XNO3 * MWNO3
                ANH4 = YNH4 * MWNH4
                GNO3 = ( TNO3 - XNO3 ) * MWHNO3
                GNH3 = ( TNH4 - YNH4 ) * MWNH3
                AH2O = 1000.0 * WH2O
                RETURN
              END IF
    
    1501    CONTINUE
    
    C...If after NITR iterations no solution is found, then:
    
            ASO4 = TSO4 * MWSO4
            ANO3 = 0.0
            YNH4 = TWOSO4
            ANH4 = YNH4 * MWNH4
            CALL AWATER ( IRH, TSO4, YNH4, XNO3, AH2O )
            GNO3 = TNO3 * MWHNO3
            GNH3 = ( TNH4 - YNH4 ) * MWNH3
            RETURN
    
          ELSE
           
    C......................................
    C......... Low Ammonia Case ...........
    C......................................
          
    C...coded by Dr. Francis S. Binkowski 12/8/91.(4/26/95)
    
    C...All cases covered by this logic
     
            WH2O = 0.0
            CALL AWATER ( IRH, TSO4, TNH4, TNO3, AH2O )
            WH2O = 1.0E-3 * AH2O
            ZH2O = AH2O
    
    C...convert 10**(-6) kg water/(cubic meter of air) to micrograms of water
    C...  per cubic meter of air (1000 g = 1 kg)
          
            ASO4 = TSO4 * MWSO4
            ANH4 = TNH4 * MWNH4
            ANO3 = 0.0
            GNO3 = TNO3 * MWHNO3
            GNH3 = 0.0 
    
    C...Check for zero water.      
    
            IF ( WH2O .EQ. 0.0 ) RETURN
            ZSO4 = TSO4 / WH2O 
    
    C...ZSO4 is the molality of total sulfate i.e. MSO4 + MHSO4      
    
    ccc         IF ( ZSO4 .GT. 11.0 ) THEN
    
    C...do not solve for aerosol nitrate for total sulfate molality
    C...  greater than 11.0 because the model parameters break down
    C...  greater than  9.0 because the model parameters break down
    
            IF ( ZSO4 .GT. 9.0 ) THEN   ! 18 June 97
              RETURN
            END IF
    
    C...First solve with activity coeffs of 1.0, then iterate.
    
            PHIOLD = 1.0
            GAMANA = 1.0
            GAMAS1 = 1.0
            GAMAS2 = 1.0
            GAMAAB = 1.0
            GAMOLD = 1.0
    
    C...All ammonia is considered to be aerosol ammonium. 
    
            MNH4 = TNH4 / WH2O
    
    C...MNH4 is the molality of ammonium ion.
    
            YNH4 = TNH4
          
    C...loop for iteration
     
            DO 1601 NNN = 1, 150
              NITR = NNN
    
    C...set up equilibrium constants including activities
    C...  solve the system for hplus first then sulfate & nitrate
    
              RK2SA = K2SA * GAMAS2 * GAMAS2 / ( GAMAS1 * GAMAS1 * GAMAS1 )
              RKNA = KNA / ( GAMANA * GAMANA )
              RKNWET = RKNA * WH2O       
              T21  = ZSO4 - MNH4
              T221 = ZSO4 + T21
    
    C...set up coefficients for cubic       
    
              A2 = RK2SA + RKNWET - T21
              A1 = RK2SA * RKNWET - T21 * ( RK2SA + RKNWET )
         &       - RK2SA * ZSO4 - RKNA * TNO3
              A0 = - (T21 * RK2SA * RKNWET
         &       + RK2SA * RKNWET * ZSO4 + RK2SA * RKNA * TNO3 )
              CALL CUBIC ( A2, A1, A0, NR, CRUTES )
           
    C...Code assumes the smallest positive root is in CRUTES(1)
     
              HPLUS = CRUTES( 1 )
              BAL = HPLUS **3 + A2 * HPLUS**2 + A1 * HPLUS + A0
              MSO4 = RK2SA * ZSO4 / ( HPLUS + RK2SA )   ! molality of sulfate ion
              MHSO4 = ZSO4 - MSO4                       ! molality of bisulfate ion
              MNA = RKNA * TNO3 / ( HPLUS + RKNWET )    ! molality of nitrate ion
              MNA = MAX( 0.0, MNA )
              MNA = MIN( MNA, TNO3 / WH2O )
              XNO3 = MNA * WH2O
              ANO3 = MNA * WH2O * MWNO3
              GNO3 = ( TNO3 - XNO3 ) * MWHNO3
            
    C...Calculate ionic strength      
    
              STION = 0.5 * ( HPLUS + MNA + MNH4 + MHSO4 + 4.0 * MSO4 )
              
    C...Update water
    
              CALL AWATER ( IRH, TSO4, YNH4, XNO3, AH2O )
    
    C...Convert 10**(-6) kg water/(cubic meter of air) to micrograms of water
    C...  per cubic meter of air (1000 g = 1 kg)                       
    
              WH2O = 1.0E-3 * AH2O 
              CAT( 1 ) = HPLUS
              CAT( 2 ) = MNH4
              AN ( 1 ) = MSO4
              AN ( 2 ) = MNA
              AN ( 3 ) = MHSO4
    
              CALL ACTCOF ( CAT, AN, GAMS, MOLNU, PHIBAR )
    
              GAMANA = GAMS( 1, 2 )
              GAMAS1 = GAMS( 1, 1 )
              GAMAS2 = GAMS( 1, 3 )
              GAMAAN = GAMS( 2, 2 )
    
              GAMAHAT = ( GAMAS2 * GAMAS2 / ( GAMAAB * GAMAAB ) )
              BHAT = KHAT * GAMAHAT 
    ccc          EROR = ABS ( ( PHIOLD - PHIBAR ) / PHIOLD )
    ccc          PHIOLD = PHIBAR
              EROR = ABS ( GAMOLD - GAMAHAT ) / GAMOLD 
              GAMOLD = GAMAHAT   
    
    C...write out molalities and activity coefficient
    C...  and return with good solution
    
              IF ( EROR .LE. TOLER2 ) THEN
    ccc            WRITE(12,*) RH, STION,HPLUS,ZSO4,MSO4,MHSO4,MNH4,MNA
    ccc            WRITE(11,*) RH, STION, GAMS(1,1),GAMS(1,2),GAMS(1,3),
    ccc     &                  GAMS(2,1),GAMS(2,2),GAMS(2,3), PHIBAR
                RETURN
              END IF
    
    1601    CONTINUE     
    
    C...after NITR iterations, failure to solve the system, no ANO3
    
            GNO3 = TNO3 * MWHNO3
            ANO3 = 0.0
            CALL AWATER ( IRH, TSO4, TNH4, TNO3, AH2O )      
            RETURN
                
          END IF   ! ratio .gt. 2.0
    
          END ! end RPMares
    C ///////////////////////////////////////////////////
    C ///////////////////////////////////////////////////
           SUBROUTINE ACTCOF ( CAT, AN, GAMA, MOLNU, PHIMULT )
    
    C-----------------------------------------------------------------------
    C
    C DESCRIPTION:
    C
    C  This subroutine computes the activity coefficients of (2NH4+,SO4--),
    C  (NH4+,NO3-),(2H+,SO4--),(H+,NO3-),AND (H+,HSO4-) in aqueous
    C  multicomponent solution, using Bromley's model and Pitzer's method.
    C
    C REFERENCES:
    C
    C   Bromley, L.A. (1973) Thermodynamic properties of strong electrolytes
    C     in aqueous solutions.  AIChE J. 19, 313-320.
    C
    C   Chan, C.K. R.C. Flagen, & J.H.  Seinfeld (1992) Water Activities of 
    C     NH4NO3 / (NH4)2SO4 solutions, Atmos. Environ. (26A): 1661-1673.
    C
    C   Clegg, S.L. & P. Brimblecombe (1988) Equilibrium partial pressures 
    C     of strong acids over saline solutions - I HNO3, 
    C     Atmos. Environ. (22): 91-100
    C
    C   Clegg, S.L. & P. Brimblecombe (1990) Equilibrium partial pressures
    C     and mean activity and osmotic coefficients of 0-100% nitric acid 
    C     as a function of temperature,   J. Phys. Chem (94): 5369 - 5380
    C
    C   Pilinis, C. and J.H. Seinfeld (1987) Continued development of a
    C     general equilibrium model for inorganic multicomponent atmospheric
    C     aerosols.  Atmos. Environ. 21(11), 2453-2466.
    C
    
    
    C
    C ARGUMENT DESCRIPTION:
    C
    C     CAT(1) : conc. of H+    (moles/kg)
    C     CAT(2) : conc. of NH4+  (moles/kg)
    C     AN(1)  : conc. of SO4-- (moles/kg)
    C     AN(2)  : conc. of NO3-  (moles/kg)
    C     AN(3)  : conc. of HSO4- (moles/kg)
    C     GAMA(2,1)    : mean molal ionic activity coeff for (2NH4+,SO4--)
    C     GAMA(2,2)    :  "    "     "       "       "    "  (NH4+,NO3-)
    C     GAMA(2,3)    :  "    "     "       "       "    "  (NH4+. HSO4-)
    C     GAMA(1,1)    :  "    "     "       "       "    "  (2H+,SO4--)
    C     GAMA(1,2)    :  "    "     "       "       "    "  (H+,NO3-)
    C     GAMA(1,3)    :  "    "     "       "       "    "  (H+,HSO4-)
    C     MOLNU   : the total number of moles of all ions.
    C     PHIMULT : the multicomponent paractical osmotic coefficient.
    C
    C REVISION HISTORY:
    C      Who       When        Detailed description of changes
    C   ---------   --------  -------------------------------------------
    C   S.Roselle   7/26/89   Copied parts of routine BROMLY, and began this
    C                         new routine using a method described by Pilinis
    C                         and Seinfeld 1987, Atmos. Envirn. 21 pp2453-2466.
    C   S.Roselle   7/30/97   Modified for use in Models-3
    C   F.Binkowski 8/7/97    Modified coefficients BETA0, BETA1, CGAMA
    C
    C-----------------------------------------------------------------------
    
          IMPLICIT NONE
    
    C...........INCLUDES and their descriptions
    
    c      INCLUDE SUBST_XSTAT     ! M3EXIT status codes
    C....................................................................
    
          INTEGER    XSTAT0       ! Normal, successful completion
          PARAMETER (XSTAT0 = 0)
          INTEGER    XSTAT1       ! File I/O error
          PARAMETER (XSTAT1 = 1)
          INTEGER    XSTAT2       ! Execution error
          PARAMETER (XSTAT2 = 2)
          INTEGER    XSTAT3       ! Special  error
          PARAMETER (XSTAT3 = 3)
    
          CHARACTER*120 XMSG
    
    C...........PARAMETERS and their descriptions:
    
          INTEGER      NCAT                 ! number of cations
          PARAMETER  ( NCAT = 2 )
    
          INTEGER      NAN                  ! number of anions
          PARAMETER  ( NAN = 3 )
    
    C...........ARGUMENTS and their descriptions
    
          REAL         MOLNU                ! tot # moles of all ions
          REAL         PHIMULT              ! multicomponent paractical osmotic coef
          REAL         CAT( NCAT )          ! cation conc in moles/kg (input)
          REAL         AN ( NAN )           ! anion conc in moles/kg (input)
          REAL         GAMA( NCAT, NAN )    ! mean molal ionic activity coefs
    
    C...........SCRATCH LOCAL VARIABLES and their descriptions:
    
          CHARACTER*16 PNAME            ! driver program name
          SAVE         PNAME
    
          INTEGER      IAN                  ! anion indX
          INTEGER      ICAT                 ! cation indX
    
          REAL         FGAMA                ! 
          REAL         I                    ! ionic strength 
          REAL         R                    ! 
          REAL         S                    ! 
          REAL         TA                   ! 
          REAL         TB                   ! 
          REAL         TC                   ! 
          REAL         TEXPV                ! 
          REAL         TRM                  ! 
          REAL         TWOI                 ! 2*ionic strength
          REAL         TWOSRI               ! 2*sqrt of ionic strength
          REAL         ZBAR                 ! 
          REAL         ZBAR2                ! 
          REAL         ZOT1                 ! 
          REAL         SRI                  ! square root of ionic strength 
          REAL         F2( NCAT )           ! 
          REAL         F1( NAN )            ! 
          REAL         ZP( NCAT )           ! absolute value of charges of cation
          REAL         ZM( NAN )            ! absolute value of charges of anion
          REAL         BGAMA ( NCAT, NAN )  ! 
          REAL         X     ( NCAT, NAN )  ! 
          REAL         M     ( NCAT, NAN )  ! molality of each electrolyte
          REAL         LGAMA0( NCAT, NAN )  ! binary activity coefficients
          REAL         Y     ( NAN, NCAT )  ! 
          REAL         BETA0 ( NCAT, NAN )  ! binary activity coefficient parameter
          REAL         BETA1 ( NCAT, NAN )  ! binary activity coefficient parameter
          REAL         CGAMA ( NCAT, NAN )  ! binary activity coefficient parameter
          REAL         V1    ( NCAT, NAN )  ! number of cations in electrolyte formula
          REAL         V2    ( NCAT, NAN )  ! number of anions in electrolyte formula
    
          DATA         ZP / 1.0, 1.0 /
          DATA         ZM / 2.0, 1.0, 1.0 /
          DATA          XMSG / ' ' /
          DATA         PNAME / 'ACTCOF' /
    
    C *** Sources for the coefficients BETA0, BETA1, CGAMA:
     
    C *** (1,1);(1,3)  - Clegg & Brimblecombe (1988)
    C *** (2,3)        - Pilinis & Seinfeld (1987), cgama different 
    C *** (1,2)        - Clegg & Brimblecombe (1990)
    C *** (2,1);(2,2)  - Chan, Flagen & Seinfeld (1992)                                 
          
    c *** now set the basic constants, BETA0, BETA1, CGAMA 
    
          DATA BETA0(1,1) /2.98E-2/,    BETA1(1,1) / 0.0/, 
         &     CGAMA(1,1) / 4.38E-2/                                ! 2H+SO4-
         
          DATA BETA0(1,2) / 1.2556E-1/,   BETA1(1,2) / 2.8778E-1/, 
         &     CGAMA(1,2) / -5.59E-3/                               ! HNO3
         
          DATA BETA0(1,3) / 2.0651E-1/,   BETA1(1,3) / 5.556E-1/,
         &     CGAMA(1,3) /0.0/                                     ! H+HSO4-
         
          DATA BETA0(2,1) /4.6465E-2/,   BETA1(2,1) /-0.54196/, 
         &     CGAMA(2,1) /-1.2683E-3/                              ! (NH4)2SO4
         
          DATA BETA0(2,2) /-7.26224E-3/, BETA1(2,2) /-1.168858/, 
         &     CGAMA(2,2) /3.51217E-5/                              ! NH4NO3
         
          DATA BETA0(2,3) / 4.494E-2/,    BETA1(2,3) / 2.3594E-1/,
         &     CGAMA(2,3) /-2.962E-3/                               ! NH4HSO4
    
          DATA V1(1,1), V2(1,1) / 2.0, 1.0 /     ! 2H+SO4-
          DATA V1(2,1), V2(2,1) / 2.0, 1.0 /     ! (NH4)2SO4
          DATA V1(1,2), V2(1,2) / 1.0, 1.0 /     ! HNO3 
          DATA V1(2,2), V2(2,2) / 1.0, 1.0 /     ! NH4NO3
          DATA V1(1,3), V2(1,3) / 1.0, 1.0 /     ! H+HSO4-
          DATA V1(2,3), V2(2,3) / 1.0, 1.0 /     ! NH4HSO4
    
    C-----------------------------------------------------------------------
    C  begin body of subroutine ACTCOF
    
    C...compute ionic strength
    
          I = 0.0
    
          DO ICAT = 1, NCAT
            I = I + CAT( ICAT ) * ZP( ICAT ) * ZP( ICAT )
          END DO
    
          DO IAN = 1, NAN
            I = I + AN( IAN ) * ZM( IAN ) * ZM( IAN )
          END DO
    
          I = 0.5 * I
    
    C...check for problems in the ionic strength
    
          IF ( I .EQ. 0.0 ) THEN
    
            DO IAN = 1, NAN
              DO ICAT = 1, NCAT
                GAMA( ICAT, IAN ) = 0.0
              END DO
            END DO
            
    c        XMSG = 'Ionic strength is zero...returning zero activities'
    c        CALL M3WARN ( PNAME, 0, 0, XMSG )
            write(*,*)  'Ionic strength is zero...returning zero activities'
            RETURN
    
          ELSE IF ( I .LT. 0.0 ) THEN
    c        XMSG = 'Ionic strength below zero...negative concentrations'
    c        CALL M3EXIT ( PNAME, 0, 0, XMSG, XSTAT1 )
    CALL ABORT
           STOP 'Ionic strength below zero...negative concentrations'
          END IF
    
    C...compute some essential expressions
    
          SRI    = SQRT( I )
          TWOSRI = 2.0 * SRI
          TWOI   = 2.0 * I
          TEXPV  = 1.0 - EXP( -TWOSRI ) * ( 1.0 + TWOSRI - TWOI )
          R      = 1.0 + 0.75 * I
          S      = 1.0 + 1.5  * I
          ZOT1   = 0.511 * SRI / ( 1.0 + SRI )
    
    C...Compute binary activity coeffs
    
          FGAMA = -0.392 * ( ( SRI / ( 1.0 + 1.2 * SRI )
         &      + ( 2.0 / 1.2 ) * ALOG( 1.0 + 1.2 * SRI ) ) )
    
          DO ICAT = 1, NCAT
            DO IAN = 1, NAN
    
            BGAMA( ICAT, IAN ) = 2.0 * BETA0( ICAT, IAN )
         &                    + ( 2.0 * BETA1( ICAT, IAN ) / ( 4.0 * I ) )
         &                    * TEXPV
    
    C...compute the molality of each electrolyte for given ionic strength
    
              M( ICAT, IAN ) = ( CAT( ICAT )**V1( ICAT, IAN )
         &                   *   AN( IAN )**V2( ICAT, IAN ) )**( 1.0
         &                   / ( V1( ICAT, IAN ) + V2( ICAT, IAN ) ) )
    
    C...calculate the binary activity coefficients
    
              LGAMA0( ICAT, IAN ) = ( ZP( ICAT ) * ZM( IAN ) * FGAMA
         &                      + M( ICAT, IAN )
         &                      * ( 2.0 * V1( ICAT, IAN ) * V2( ICAT, IAN )
         &                      / ( V1( ICAT, IAN ) + V2( ICAT, IAN ) )