chem/isocom.F

   1 !=======================================================================
   2 !
   3 ! *** ISORROPIA CODE
   4 ! *** SUBROUTINE ISOROPIA
   5 ! *** THIS SUBROUTINE IS THE MASTER ROUTINE FOR THE ISORROPIA
   6 !     THERMODYNAMIC EQUILIBRIUM AEROSOL MODEL (VERSION 1.1 and above)
   7 !
   8 ! ======================== ARGUMENTS / USAGE ===========================
   9 !
  10 !  INPUT:
  11 !  1. [WI]
  12 !     DOUBLE PRECISION array of length [5].
  13 !     Concentrations, expressed in moles/m3. Depending on the type of
  14 !     problem solved (specified in CNTRL(1)), WI contains either
  15 !     GAS+AEROSOL or AEROSOL only concentratios.
  16 !     WI(1) - sodium
  17 !     WI(2) - sulfate
  18 !     WI(3) - ammonium
  19 !     WI(4) - nitrate
  20 !     WI(5) - chloride
  21 !
  22 !  2. [RHI]
  23 !     DOUBLE PRECISION variable.
  24 !     Ambient relative humidity expressed on a (0,1) scale.
  25 !
  26 !  3. [TEMPI]
  27 !     DOUBLE PRECISION variable.
  28 !     Ambient temperature expressed in Kelvins.
  29 !
  30 !  4. [CNTRL]
  31 !     DOUBLE PRECISION array of length [2].
  32 !     Parameters that control the type of problem solved.
  33 !
  34 !     CNTRL(1): Defines the type of problem solved.
  35 !     0 - Forward problem is solved. In this case, array WI contains
  36 !         GAS and AEROSOL concentrations together.
  37 !     1 - Reverse problem is solved. In this case, array WI contains
  38 !         AEROSOL concentrations only.
  39 !
  40 !     CNTRL(2): Defines the state of the aerosol
  41 !     0 - The aerosol can have both solid+liquid phases (deliquescent)
  42 !     1 - The aerosol is in only liquid state (metastable aerosol)
  43 !
  44 !  OUTPUT:
  45 !  1. [WT]
  46 !     DOUBLE PRECISION array of length [5].
  47 !     Total concentrations (GAS+AEROSOL) of species, expressed in moles/m3.
  48 !     If the foreward probelm is solved (CNTRL(1)=0), array WT is
  49 !     identical to array WI.
  50 !     WT(1) - total sodium
  51 !     WT(2) - total sulfate
  52 !     WT(3) - total ammonium
  53 !     WT(4) - total nitrate
  54 !     WT(5) - total chloride
  55 !
  56 !  2. [GAS]
  57 !     DOUBLE PRECISION array of length [03].
  58 !     Gaseous species concentrations, expressed in moles/m3.
  59 !     GAS(1) - NH3
  60 !     GAS(2) - HNO3
  61 !     GAS(3) - HCl
  62 !
  63 !  3. [AERLIQ]
  64 !     DOUBLE PRECISION array of length [11].
  65 !     Liquid aerosol species concentrations, expressed in moles/m3.
  66 !     AERLIQ(01) - H+(aq)
  67 !     AERLIQ(02) - Na+(aq)
  68 !     AERLIQ(03) - NH4+(aq)
  69 !     AERLIQ(04) - Cl-(aq)
  70 !     AERLIQ(05) - SO4--(aq)
  71 !     AERLIQ(06) - HSO4-(aq)
  72 !     AERLIQ(07) - NO3-(aq)
  73 !     AERLIQ(08) - H2O
  74 !     AERLIQ(09) - NH3(aq) (undissociated)
  75 !     AERLIQ(10) - HNCl(aq) (undissociated)
  76 !     AERLIQ(11) - HNO3(aq) (undissociated)
  77 !     AERLIQ(12) - OH-(aq)
  78 !
  79 !  4. [AERSLD]
  80 !     DOUBLE PRECISION array of length [09].
  81 !     Solid aerosol species concentrations, expressed in moles/m3.
  82 !     AERSLD(01) - NaNO3(s)
  83 !     AERSLD(02) - NH4NO3(s)
  84 !     AERSLD(03) - NaCl(s)
  85 !     AERSLD(04) - NH4Cl(s)
  86 !     AERSLD(05) - Na2SO4(s)
  87 !     AERSLD(06) - (NH4)2SO4(s)
  88 !     AERSLD(07) - NaHSO4(s)
  89 !     AERSLD(08) - NH4HSO4(s)
  90 !     AERSLD(09) - (NH4)4H(SO4)2(s)
  91 !
  92 !  5. [SCASI]
  93 !     CHARACTER*15 variable.
  94 !     Returns the subcase which the input corresponds to.
  95 !
  96 !  6. [OTHER]
  97 !     DOUBLE PRECISION array of length [6].
  98 !     Returns solution information.
  99 !
 100 !     OTHER(1): Shows if aerosol water exists.
 101 !     0 - Aerosol is WET
 102 !     1 - Aerosol is DRY
 103 !
 104 !     OTHER(2): Aerosol Sulfate ratio, defined as (in moles/m3) :
 105 !               (total ammonia + total Na) / (total sulfate)
 106 !
 107 !     OTHER(3): Sulfate ratio based on aerosol properties that defines
 108 !               a sulfate poor system:
 109 !               (aerosol ammonia + aerosol Na) / (aerosol sulfate)
 110 !
 111 !     OTHER(4): Aerosol sodium ratio, defined as (in moles/m3) :
 112 !               (total Na) / (total sulfate)
 113 !
 114 !     OTHER(5): Ionic strength of the aqueous aerosol (if it exists).
 115 !
 116 !     OTHER(6): Total number of calls to the activity coefficient
 117 !               calculation subroutine.
 118 !
 119 ! *** COPYRIGHT 1996-2006, UNIVERSITY OF MIAMI, CARNEGIE MELLON UNIVERSITY,
 120 ! *** GEORGIA INSTITUTE OF TECHNOLOGY
 121 ! *** WRITTEN BY ATHANASIOS NENES
 122 ! *** UPDATED BY CHRISTOS FOUNTOUKIS
 123 !
 124 !=======================================================================
 125 !
 126       SUBROUTINE ISOROPIA (WI, RHI, TEMPI,  CNTRL,    &
 127                            WT, GAS, AERLIQ, AERSLD, SCASI, OTHER)
 128       USE ISRPIA
 129       IMPLICIT DOUBLE PRECISION (A-H, O-Z)
 130 !      INCLUDE 'ISRPIA.INC'
 131       PARAMETER (NCTRL=2,NOTHER=6)
 132       CHARACTER SCASI*15
 133       DIMENSION WI(NCOMP), WT(NCOMP),   GAS(NGASAQ),  AERSLD(NSLDS),    &
 134                 AERLIQ(NIONS+NGASAQ+2), CNTRL(NCTRL), OTHER(NOTHER)
 135 !
 136 ! *** PROBLEM TYPE (0=FOREWARD, 1=REVERSE) ******************************
 137 !
 138       IPROB   = NINT(CNTRL(1))
 139 !
 140 ! *** AEROSOL STATE (0=SOLID+LIQUID, 1=METASTABLE) **********************
 141 !
 142       METSTBL = NINT(CNTRL(2))
 143 !
 144 ! *** SOLVE FOREWARD PROBLEM ********************************************
 145 !
 146 50    IF (IPROB.EQ.0) THEN
 147          IF (WI(1)+WI(2)+WI(3)+WI(4)+WI(5) .LE. TINY) THEN ! Everything=0
 148             CALL INIT1 (WI, RHI, TEMPI)
 149          ELSE IF (WI(1)+WI(4)+WI(5) .LE. TINY) THEN        ! Na,Cl,NO3=0
 150             CALL ISRP1F (WI, RHI, TEMPI)
 151          ELSE IF (WI(1)+WI(5) .LE. TINY) THEN              ! Na,Cl=0
 152             CALL ISRP2F (WI, RHI, TEMPI)
 153          ELSE
 154             CALL ISRP3F (WI, RHI, TEMPI)
 155          ENDIF
 156 !
 157 ! *** SOLVE REVERSE PROBLEM *********************************************
 158 !
 159       ELSE
 160          IF (WI(1)+WI(2)+WI(3)+WI(4)+WI(5) .LE. TINY) THEN ! Everything=0
 161             CALL INIT1 (WI, RHI, TEMPI)
 162          ELSE IF (WI(1)+WI(4)+WI(5) .LE. TINY) THEN        ! Na,Cl,NO3=0
 163             CALL ISRP1R (WI, RHI, TEMPI)
 164          ELSE IF (WI(1)+WI(5) .LE. TINY) THEN              ! Na,Cl=0
 165             CALL ISRP2R (WI, RHI, TEMPI)
 166          ELSE
 167             CALL ISRP3R (WI, RHI, TEMPI)
 168          ENDIF
 169       ENDIF
 170 !
 171 ! *** ADJUST MASS BALANCE ***********************************************
 172 !
 173       IF (NADJ.EQ.1) CALL ADJUST (WI)
 174 !cC
 175 !cC *** IF METASTABLE AND NO WATER - RESOLVE AS NORMAL ********************
 176 !cC
 177 !c      IF (WATER.LE.TINY .AND. METSTBL.EQ.1) THEN
 178 !c         METSTBL = 0
 179 !c         GOTO 50
 180 !c      ENDIF
 181 !
 182 ! *** SAVE RESULTS TO ARRAYS (units = mole/m3) ****************************
 183 !
 184       GAS(1) = GNH3                ! Gaseous aerosol species
 185       GAS(2) = GHNO3
 186       GAS(3) = GHCL
 187 !
 188       DO 10 I=1,NIONS              ! Liquid aerosol species
 189          AERLIQ(I) = MOLAL(I)
 190   10  CONTINUE
 191       DO 20 I=1,NGASAQ
 192          AERLIQ(NIONS+1+I) = GASAQ(I)
 193   20  CONTINUE
 194       AERLIQ(NIONS+1)        = WATER*1.0D3/18.0D0
 195       AERLIQ(NIONS+NGASAQ+2) = COH
 196 !
 197       AERSLD(1) = CNANO3           ! Solid aerosol species
 198       AERSLD(2) = CNH4NO3
 199       AERSLD(3) = CNACL
 200       AERSLD(4) = CNH4CL
 201       AERSLD(5) = CNA2SO4
 202       AERSLD(6) = CNH42S4
 203       AERSLD(7) = CNAHSO4
 204       AERSLD(8) = CNH4HS4
 205       AERSLD(9) = CLC
 206 !
 207       IF(WATER.LE.TINY) THEN       ! Dry flag
 208         OTHER(1) = 1.d0
 209       ELSE
 210         OTHER(1) = 0.d0
 211       ENDIF
 212 !
 213       OTHER(2) = SULRAT            ! Other stuff
 214       OTHER(3) = SULRATW
 215       OTHER(4) = SODRAT
 216       OTHER(5) = IONIC
 217       OTHER(6) = ICLACT
 218 !
 219       SCASI = SCASE
 220 !
 221       WT(1) = WI(1)                ! Total gas+aerosol phase
 222       WT(2) = WI(2)
 223       WT(3) = WI(3)
 224       WT(4) = WI(4)
 225       WT(5) = WI(5)
 226       IF (IPROB.GT.0 .AND. WATER.GT.TINY) THEN
 227          WT(3) = WT(3) + GNH3
 228          WT(4) = WT(4) + GHNO3
 229          WT(5) = WT(5) + GHCL
 230       ENDIF
 231 !
 232       RETURN
 233 !
 234 ! *** END OF SUBROUTINE ISOROPIA ******************************************
 235 !
 236       END
 237 !=======================================================================
 238 !
 239 ! *** ISORROPIA CODE
 240 ! *** SUBROUTINE SETPARM
 241 ! *** THIS SUBROUTINE REDEFINES THE SOLUTION PARAMETERS OF ISORROPIA
 242 !
 243 ! ======================== ARGUMENTS / USAGE ===========================
 244 !
 245 ! *** NOTE: IF NEGATIVE VALUES ARE GIVEN FOR A PARAMETER, IT IS
 246 !     IGNORED AND THE CURRENT VALUE IS USED INSTEAD.
 247 !
 248 !  INPUT:
 249 !  1. [WFTYPI]
 250 !     INTEGER variable.
 251 !     Defines the type of weighting algorithm for the solution in Mutual
 252 !     Deliquescence Regions (MDR's):
 253 !     0 - MDR's are assumed dry. This is equivalent to the approach
 254 !         used by SEQUILIB.
 255 !     1 - The solution is assumed "half" dry and "half" wet throughout
 256 !         the MDR.
 257 !     2 - The solution is a relative-humidity weighted mean of the
 258 !         dry and wet solutions (as defined in Nenes et al., 1998)
 259 !
 260 !  2. [IACALCI]
 261 !     INTEGER variable.
 262 !     Method of activity coefficient calculation:
 263 !     0 - Calculate coefficients during runtime
 264 !     1 - Use precalculated tables
 265 !
 266 !  3. [EPSI]
 267 !     DOUBLE PRECITION variable.
 268 !     Defines the convergence criterion for all iterative processes
 269 !     in ISORROPIA, except those for activity coefficient calculations
 270 !     (EPSACTI controls that).
 271 !
 272 !  4. [MAXITI]
 273 !     INTEGER variable.
 274 !     Defines the maximum number of iterations for all iterative
 275 !     processes in ISORROPIA, except for activity coefficient calculations
 276 !     (NSWEEPI controls that).
 277 !
 278 !  5. [NSWEEPI]
 279 !     INTEGER variable.
 280 !     Defines the maximum number of iterations for activity coefficient
 281 !     calculations.
 282 !
 283 !  6. [EPSACTI]
 284 !     DOUBLE PRECISION variable.
 285 !     Defines the convergence criterion for activity coefficient
 286 !     calculations.
 287 !
 288 !  7. [NDIV]
 289 !     INTEGER variable.
 290 !     Defines the number of subdivisions needed for the initial root
 291 !     tracking for the bisection method. Usually this parameter should
 292 !     not be altered, but is included for completeness.
 293 !
 294 !  8. [NADJ]
 295 !     INTEGER variable.
 296 !     Forces the solution obtained to satisfy total mass balance
 297 !     to machine precision
 298 !     0 - No adjustment done (default)
 299 !     1 - Do adjustment
 300 !
 301 ! *** COPYRIGHT 1996-2006, UNIVERSITY OF MIAMI, CARNEGIE MELLON UNIVERSITY,
 302 ! *** GEORGIA INSTITUTE OF TECHNOLOGY
 303 ! *** WRITTEN BY ATHANASIOS NENES
 304 ! *** UPDATED BY CHRISTOS FOUNTOUKIS
 305 !
 306 !=======================================================================
 307 !
 308       SUBROUTINE SETPARM (WFTYPI,  IACALCI, EPSI, MAXITI, NSWEEPI,    &
 309                           EPSACTI, NDIVI, NADJI)
 310       USE ISRPIA
 311       IMPLICIT DOUBLE PRECISION (A-H, O-Z)
 312 !      INCLUDE 'ISRPIA.INC'
 313       INTEGER  WFTYPI
 314 !
 315 ! *** SETUP SOLUTION PARAMETERS *****************************************
 316 !
 317       IF (WFTYPI .GE. 0)   WFTYP  = WFTYPI
 318       IF (IACALCI.GE. 0)   IACALC = IACALCI
 319       IF (EPSI   .GE.ZERO) EPS    = EPSI
 320       IF (MAXITI .GT. 0)   MAXIT  = MAXITI
 321       IF (NSWEEPI.GT. 0)   NSWEEP = NSWEEPI
 322       IF (EPSACTI.GE.ZERO) EPSACT = EPSACTI
 323       IF (NDIVI  .GT. 0)   NDIV   = NDIVI
 324       IF (NADJI  .GE. 0)   NADJ   = NADJI
 325 !
 326 ! *** END OF SUBROUTINE SETPARM *****************************************
 327 !
 328       RETURN
 329       END
 330
 331
 332
 333
 334 !=======================================================================
 335 !
 336 ! *** ISORROPIA CODE
 337 ! *** SUBROUTINE GETPARM
 338 ! *** THIS SUBROUTINE OBTAINS THE CURRENT VAULES OF THE SOLUTION
 339 !     PARAMETERS OF ISORROPIA
 340 !
 341 ! ======================== ARGUMENTS / USAGE ===========================
 342 !
 343 ! *** THE PARAMETERS ARE THOSE OF SUBROUTINE SETPARM
 344 !
 345 ! *** COPYRIGHT 1996-2006, UNIVERSITY OF MIAMI, CARNEGIE MELLON UNIVERSITY,
 346 ! *** GEORGIA INSTITUTE OF TECHNOLOGY
 347 ! *** WRITTEN BY ATHANASIOS NENES
 348 ! *** UPDATED BY CHRISTOS FOUNTOUKIS
 349 !
 350 !=======================================================================
 351 !
 352       SUBROUTINE GETPARM (WFTYPI,  IACALCI, EPSI, MAXITI, NSWEEPI,    &
 353                           EPSACTI, NDIVI, NADJI)
 354       USE ISRPIA
 355       IMPLICIT DOUBLE PRECISION (A-H, O-Z)
 356 !      INCLUDE 'ISRPIA.INC'
 357       INTEGER  WFTYPI
 358 !
 359 ! *** GET SOLUTION PARAMETERS *******************************************
 360 !
 361       WFTYPI  = WFTYP
 362       IACALCI = IACALC
 363       EPSI    = EPS
 364       MAXITI  = MAXIT
 365       NSWEEPI = NSWEEP
 366       EPSACTI = EPSACT
 367       NDIVI   = NDIV
 368       NADJI   = NADJ
 369 !
 370 ! *** END OF SUBROUTINE GETPARM *****************************************
 371 !
 372       RETURN
 373       END
 374
 375 !=======================================================================
 376 !
 377 ! *** ISORROPIA CODE
 378 ! *** BLOCK DATA BLKISO
 379 ! *** THIS SUBROUTINE PROVIDES INITIAL (DEFAULT) VALUES TO PROGRAM
 380 !     PARAMETERS VIA DATA STATEMENTS
 381 !
 382 ! *** COPYRIGHT 1996-2006, UNIVERSITY OF MIAMI, CARNEGIE MELLON UNIVERSITY,
 383 ! *** GEORGIA INSTITUTE OF TECHNOLOGY
 384 ! *** WRITTEN BY ATHANASIOS NENES
 385 ! *** UPDATED BY CHRISTOS FOUNTOUKIS
 386 !
 387 ! *** ZSR RELATIONSHIP PARAMETERS MODIFIED BY DOUGLAS WALDRON
 388 ! *** OCTOBER 2003
 389 ! *** BASED ON AIM MODEL III (http://mae.ucdavis.edu/wexler/aim)
 390 !
 391 !=======================================================================
 392 !
 393       BLOCK DATA BLKISO
 394       USE ISRPIA
 395       IMPLICIT DOUBLE PRECISION (A-H, O-Z)
 396 !      INCLUDE 'ISRPIA.INC'
 397 !
 398 ! *** DEFAULT VALUES *************************************************
 399 ! *** OTHER PARAMETERS ***********************************************
 400 !
 401 ! *** ZSR RELATIONSHIP PARAMETERS **************************************
 402 ! *** END OF BLOCK DATA SUBPROGRAM *************************************
 403 !
 404       END
 405 !=======================================================================
 406 !
 407 ! *** ISORROPIA CODE
 408 ! *** SUBROUTINE INIT1
 409 ! *** THIS SUBROUTINE INITIALIZES ALL GLOBAL VARIABLES FOR AMMONIUM
 410 !     SULFATE AEROSOL SYSTEMS (SUBROUTINE ISRP1)
 411 !
 412 ! *** COPYRIGHT 1996-2006, UNIVERSITY OF MIAMI, CARNEGIE MELLON UNIVERSITY,
 413 ! *** GEORGIA INSTITUTE OF TECHNOLOGY
 414 ! *** WRITTEN BY ATHANASIOS NENES
 415 ! *** UPDATED BY CHRISTOS FOUNTOUKIS
 416 !
 417 !=======================================================================
 418 !
 419       SUBROUTINE INIT1 (WI, RHI, TEMPI)
 420       USE ISRPIA
 421       IMPLICIT DOUBLE PRECISION (A-H, O-Z)
 422 !      INCLUDE 'ISRPIA.INC'
 423       DIMENSION WI(NCOMP)
 424       REAL      IC,GII,GI0,XX,LN10
 425       PARAMETER (LN10=2.3025851)
 426 !
 427 ! *** SAVE INPUT VARIABLES IN COMMON BLOCK ******************************
 428 !
 429       IF (IPROB.EQ.0) THEN                 ! FORWARD CALCULATION
 430          DO 10 I=1,NCOMP
 431             W(I) = MAX(WI(I), TINY)
 432 10       CONTINUE
 433       ELSE
 434          DO 15 I=1,NCOMP                   ! REVERSE CALCULATION
 435             WAER(I) = MAX(WI(I), TINY)
 436             W(I)    = ZERO
 437 15       CONTINUE
 438       ENDIF
 439       RH      = RHI
 440       TEMP    = TEMPI
 441 !
 442 ! *** CALCULATE EQUILIBRIUM CONSTANTS ***********************************
 443 !
 444       XK1  = 1.015e-2  ! HSO4(aq)         <==> H(aq)     + SO4(aq)
 445       XK21 = 57.639    ! NH3(g)           <==> NH3(aq)
 446       XK22 = 1.805e-5  ! NH3(aq)          <==> NH4(aq)   + OH(aq)
 447       XK7  = 1.817     ! (NH4)2SO4(s)     <==> 2*NH4(aq) + SO4(aq)
 448       XK12 = 1.382e2   ! NH4HSO4(s)       <==> NH4(aq)   + HSO4(aq)
 449       XK13 = 29.268    ! (NH4)3H(SO4)2(s) <==> 3*NH4(aq) + HSO4(aq) + SO4(aq)
 450       XKW  = 1.010e-14 ! H2O              <==> H(aq)     + OH(aq)
 451 !
 452       IF (INT(TEMP) .NE. 298) THEN   ! FOR T != 298K or 298.15K
 453          T0  = 298.15
 454          T0T = T0/TEMP
 455          COEF= 1.0+LOG(T0T)-T0T
 456          XK1 = XK1 *EXP(  8.85*(T0T-1.0) + 25.140*COEF)
 457          XK21= XK21*EXP( 13.79*(T0T-1.0) -  5.393*COEF)
 458          XK22= XK22*EXP( -1.50*(T0T-1.0) + 26.920*COEF)
 459          XK7 = XK7 *EXP( -2.65*(T0T-1.0) + 38.570*COEF)
 460          XK12= XK12*EXP( -2.87*(T0T-1.0) + 15.830*COEF)
 461          XK13= XK13*EXP( -5.19*(T0T-1.0) + 54.400*COEF)
 462          XKW = XKW *EXP(-22.52*(T0T-1.0) + 26.920*COEF)
 463       ENDIF
 464       XK2 = XK21*XK22
 465 !
 466 ! *** CALCULATE DELIQUESCENCE RELATIVE HUMIDITIES (UNICOMPONENT) ********
 467 !
 468       DRH2SO4  = 0.0000D0
 469       DRNH42S4 = 0.7997D0
 470       DRNH4HS4 = 0.4000D0
 471       DRLC     = 0.6900D0
 472       IF (INT(TEMP) .NE. 298) THEN
 473          T0       = 298.15d0
 474          TCF      = 1.0/TEMP - 1.0/T0
 475          DRNH42S4 = DRNH42S4*EXP( 80.*TCF)
 476          DRNH4HS4 = DRNH4HS4*EXP(384.*TCF)
 477          DRLC     = DRLC    *EXP(186.*TCF)
 478       ENDIF
 479 !
 480 ! *** CALCULATE MUTUAL DELIQUESCENCE RELATIVE HUMIDITIES ****************
 481 !
 482       DRMLCAB = 0.3780D0              ! (NH4)3H(SO4)2 & NH4HSO4
 483       DRMLCAS = 0.6900D0              ! (NH4)3H(SO4)2 & (NH4)2SO4
 484 !CC      IF (INT(TEMP) .NE. 298) THEN      ! For the time being.
 485 !CC         T0       = 298.15d0
 486 !CC         TCF      = 1.0/TEMP - 1.0/T0
 487 !CC         DRMLCAB  = DRMLCAB*EXP(507.506*TCF)
 488 !CC         DRMLCAS  = DRMLCAS*EXP(133.865*TCF)
 489 !CC      ENDIF
 490 !
 491 ! *** LIQUID PHASE ******************************************************
 492 !
 493       CHNO3  = ZERO
 494       CHCL   = ZERO
 495       CH2SO4 = ZERO
 496       COH    = ZERO
 497       WATER  = TINY
 498 !
 499       DO 20 I=1,NPAIR
 500          MOLALR(I)=ZERO
 501          GAMA(I)  =0.1
 502          GAMIN(I) =GREAT
 503          GAMOU(I) =GREAT
 504          M0(I)    =1d5
 505  20   CONTINUE
 506 !
 507       DO 30 I=1,NPAIR
 508          GAMA(I) = 0.1d0
 509  30   CONTINUE
 510 !
 511       DO 40 I=1,NIONS
 512          MOLAL(I)=ZERO
 513 40    CONTINUE
 514       COH = ZERO
 515 !
 516       DO 50 I=1,NGASAQ
 517          GASAQ(I)=ZERO
 518 50    CONTINUE
 519 !
 520 ! *** SOLID PHASE *******************************************************
 521 !
 522       CNH42S4= ZERO
 523       CNH4HS4= ZERO
 524       CNACL  = ZERO
 525       CNA2SO4= ZERO
 526       CNANO3 = ZERO
 527       CNH4NO3= ZERO
 528       CNH4CL = ZERO
 529       CNAHSO4= ZERO
 530       CLC    = ZERO
 531 !
 532 ! *** GAS PHASE *********************************************************
 533 !
 534       GNH3   = ZERO
 535       GHNO3  = ZERO
 536       GHCL   = ZERO
 537 !
 538 ! *** CALCULATE ZSR PARAMETERS ******************************************
 539 !
 540       IRH    = MIN (INT(RH*NZSR+0.5),NZSR)  ! Position in ZSR arrays
 541       IRH    = MAX (IRH, 1)
 542 !
 543 !      M0(01) = AWSC(IRH)      ! NACl
 544 !      IF (M0(01) .LT. 100.0) THEN
 545 !         IC = M0(01)
 546 !         CALL KMTAB(IC,298.0,     GI0,XX,XX,XX,XX,XX,XX,XX,XX,XX,XX,XX)
 547 !         CALL KMTAB(IC,SNGL(TEMP),GII,XX,XX,XX,XX,XX,XX,XX,XX,XX,XX,XX)
 548 !         M0(01) = M0(01)*EXP(LN10*(GI0-GII))
 549 !      ENDIF
 550 !
 551 !      M0(02) = AWSS(IRH)      ! (NA)2SO4
 552 !      IF (M0(02) .LT. 100.0) THEN
 553 !         IC = 3.0*M0(02)
 554 !         CALL KMTAB(IC,298.0,     XX,GI0,XX,XX,XX,XX,XX,XX,XX,XX,XX,XX)
 555 !         CALL KMTAB(IC,SNGL(TEMP),XX,GII,XX,XX,XX,XX,XX,XX,XX,XX,XX,XX)
 556 !         M0(02) = M0(02)*EXP(LN10*(GI0-GII))
 557 !      ENDIF
 558 !
 559 !      M0(03) = AWSN(IRH)      ! NANO3
 560 !      IF (M0(03) .LT. 100.0) THEN
 561 !         IC = M0(03)
 562 !         CALL KMTAB(IC,298.0,     XX,XX,GI0,XX,XX,XX,XX,XX,XX,XX,XX,XX)
 563 !         CALL KMTAB(IC,SNGL(TEMP),XX,XX,GII,XX,XX,XX,XX,XX,XX,XX,XX,XX)
 564 !         M0(03) = M0(03)*EXP(LN10*(GI0-GII))
 565 !      ENDIF
 566 !
 567       M0(04) = AWAS(IRH)      ! (NH4)2SO4
 568 !      IF (M0(04) .LT. 100.0) THEN
 569 !         IC = 3.0*M0(04)
 570 !         CALL KMTAB(IC,298.0,     XX,XX,XX,GI0,XX,XX,XX,XX,XX,XX,XX,XX)
 571 !         CALL KMTAB(IC,SNGL(TEMP),XX,XX,XX,GII,XX,XX,XX,XX,XX,XX,XX,XX)
 572 !         M0(04) = M0(04)*EXP(LN10*(GI0-GII))
 573 !      ENDIF
 574 !
 575 !      M0(05) = AWAN(IRH)      ! NH4NO3
 576 !      IF (M0(05) .LT. 100.0) THEN
 577 !         IC     = M0(05)
 578 !         CALL KMTAB(IC,298.0,     XX,XX,XX,XX,GI0,XX,XX,XX,XX,XX,XX,XX)
 579 !         CALL KMTAB(IC,SNGL(TEMP),XX,XX,XX,XX,GII,XX,XX,XX,XX,XX,XX,XX)
 580 !         M0(05) = M0(05)*EXP(LN10*(GI0-GII))
 581 !      ENDIF
 582 !
 583 !      M0(06) = AWAC(IRH)      ! NH4CL
 584 !      IF (M0(06) .LT. 100.0) THEN
 585 !         IC = M0(06)
 586 !         CALL KMTAB(IC,298.0,     XX,XX,XX,XX,XX,GI0,XX,XX,XX,XX,XX,XX)
 587 !         CALL KMTAB(IC,SNGL(TEMP),XX,XX,XX,XX,XX,GII,XX,XX,XX,XX,XX,XX)
 588 !         M0(06) = M0(06)*EXP(LN10*(GI0-GII))
 589 !      ENDIF
 590 !
 591       M0(07) = AWSA(IRH)      ! 2H-SO4
 592 !      IF (M0(07) .LT. 100.0) THEN
 593 !         IC = 3.0*M0(07)
 594 !         CALL KMTAB(IC,298.0,     XX,XX,XX,XX,XX,XX,GI0,XX,XX,XX,XX,XX)
 595 !         CALL KMTAB(IC,SNGL(TEMP),XX,XX,XX,XX,XX,XX,GII,XX,XX,XX,XX,XX)
 596 !         M0(07) = M0(07)*EXP(LN10*(GI0-GII))
 597 !      ENDIF
 598 !
 599       M0(08) = AWSA(IRH)      ! H-HSO4
 600 !CC      IF (M0(08) .LT. 100.0) THEN     ! These are redundant, because M0(8) is not used
 601 !CC         IC = M0(08)
 602 !CC         CALL KMTAB(IC,298.0,     XX,XX,XX,XX,XX,XX,XX,GI0,XX,XX,XX,XX)
 603 !CC         CALL KMTAB(IC,SNGL(TEMP),XX,XX,XX,XX,XX,XX,XX,GI0,XX,XX,XX,XX)
 604 !CCCCC         CALL KMTAB(IC,SNGL(TEMP),XX,XX,XX,XX,XX,XX,XX,GII,XX,XX,XX,XX)
 605 !CC         M0(08) = M0(08)*EXP(LN10*(GI0-GII))
 606 !CC      ENDIF
 607 !
 608       M0(09) = AWAB(IRH)      ! NH4HSO4
 609 !      IF (M0(09) .LT. 100.0) THEN
 610 !         IC = M0(09)
 611 !         CALL KMTAB(IC,298.0,     XX,XX,XX,XX,XX,XX,XX,XX,GI0,XX,XX,XX)
 612 !         CALL KMTAB(IC,SNGL(TEMP),XX,XX,XX,XX,XX,XX,XX,XX,GII,XX,XX,XX)
 613 !         M0(09) = M0(09)*EXP(LN10*(GI0-GII))
 614 !      ENDIF
 615 !
 616 !      M0(12) = AWSB(IRH)      ! NAHSO4
 617 !      IF (M0(12) .LT. 100.0) THEN
 618 !         IC = M0(12)
 619 !         CALL KMTAB(IC,298.0,     XX,XX,XX,XX,XX,XX,XX,XX,XX,XX,XX,GI0)
 620 !         CALL KMTAB(IC,SNGL(TEMP),XX,XX,XX,XX,XX,XX,XX,XX,XX,XX,XX,GII)
 621 !         M0(12) = M0(12)*EXP(LN10*(GI0-GII))
 622 !      ENDIF
 623 !
 624       M0(13) = AWLC(IRH)      ! (NH4)3H(SO4)2
 625 !      IF (M0(13) .LT. 100.0) THEN
 626 !         IC     = 4.0*M0(13)
 627 !         CALL KMTAB(IC,298.0,     XX,XX,XX,GI0,XX,XX,XX,XX,GII,XX,XX,XX)
 628 !         G130   = 0.2*(3.0*GI0+2.0*GII)
 629 !         CALL KMTAB(IC,SNGL(TEMP),XX,XX,XX,GI0,XX,XX,XX,XX,GII,XX,XX,XX)
 630 !         G13I   = 0.2*(3.0*GI0+2.0*GII)
 631 !         M0(13) = M0(13)*EXP(LN10*SNGL(G130-G13I))
 632 !      ENDIF
 633 !
 634 ! *** OTHER INITIALIZATIONS *********************************************
 635 !
 636       ICLACT  = 0
 637       CALAOU  = .TRUE.
 638       CALAIN  = .TRUE.
 639       FRST    = .TRUE.
 640       SCASE   = '??'
 641       SULRATW = 2.D0
 642       SODRAT  = ZERO
 643       NOFER   = 0
 644       STKOFL  =.FALSE.
 645       DO 60 I=1,NERRMX
 646          ERRSTK(I) =-999
 647          ERRMSG(I) = 'MESSAGE N/A'
 648    60 CONTINUE
 649 !
 650 ! *** END OF SUBROUTINE INIT1 *******************************************
 651 !
 652       END
 653
 654
 655
 656 !=======================================================================
 657 !
 658 ! *** ISORROPIA CODE
 659 ! *** SUBROUTINE INIT2
 660 ! *** THIS SUBROUTINE INITIALIZES ALL GLOBAL VARIABLES FOR AMMONIUM,
 661 !     NITRATE, SULFATE AEROSOL SYSTEMS (SUBROUTINE ISRP2)
 662 !
 663 ! *** COPYRIGHT 1996-2006, UNIVERSITY OF MIAMI, CARNEGIE MELLON UNIVERSITY,
 664 ! *** GEORGIA INSTITUTE OF TECHNOLOGY
 665 ! *** WRITTEN BY ATHANASIOS NENES
 666 ! *** UPDATED BY CHRISTOS FOUNTOUKIS
 667 !
 668 !=======================================================================
 669 !
 670       SUBROUTINE INIT2 (WI, RHI, TEMPI)
 671       USE ISRPIA
 672       IMPLICIT DOUBLE PRECISION (A-H, O-Z)
 673 !      INCLUDE 'ISRPIA.INC'
 674       DIMENSION WI(NCOMP)
 675       REAL      IC,GII,GI0,XX,LN10
 676       PARAMETER (LN10=2.3025851)
 677 !
 678 ! *** SAVE INPUT VARIABLES IN COMMON BLOCK ******************************
 679 !
 680       IF (IPROB.EQ.0) THEN                 ! FORWARD CALCULATION
 681          DO 10 I=1,NCOMP
 682             W(I) = MAX(WI(I), TINY)
 683 10       CONTINUE
 684       ELSE
 685          DO 15 I=1,NCOMP                   ! REVERSE CALCULATION
 686             WAER(I) = MAX(WI(I), TINY)
 687             W(I)    = ZERO
 688 15       CONTINUE
 689       ENDIF
 690       RH      = RHI
 691       TEMP    = TEMPI
 692 !
 693 ! *** CALCULATE EQUILIBRIUM CONSTANTS ***********************************
 694 !
 695       XK1  = 1.015e-2  ! HSO4(aq)         <==> H(aq)     + SO4(aq)
 696       XK21 = 57.639    ! NH3(g)           <==> NH3(aq)
 697       XK22 = 1.805e-5  ! NH3(aq)          <==> NH4(aq)   + OH(aq)
 698       XK4  = 2.511e6   ! HNO3(g)          <==> H(aq)     + NO3(aq) ! ISORR
 699 !CC      XK4  = 3.638e6   ! HNO3(g)          <==> H(aq)     + NO3(aq) ! SEQUIL
 700       XK41 = 2.100e5   ! HNO3(g)          <==> HNO3(aq)
 701       XK7  = 1.817     ! (NH4)2SO4(s)     <==> 2*NH4(aq) + SO4(aq)
 702       XK10 = 5.746e-17 ! NH4NO3(s)        <==> NH3(g)    + HNO3(g) ! ISORR
 703 !CC      XK10 = 2.985e-17 ! NH4NO3(s)        <==> NH3(g)    + HNO3(g) ! SEQUIL
 704       XK12 = 1.382e2   ! NH4HSO4(s)       <==> NH4(aq)   + HSO4(aq)
 705       XK13 = 29.268    ! (NH4)3H(SO4)2(s) <==> 3*NH4(aq) + HSO4(aq) + SO4(aq)
 706       XKW  = 1.010e-14 ! H2O              <==> H(aq)     + OH(aq)
 707 !
 708       IF (INT(TEMP) .NE. 298) THEN   ! FOR T != 298K or 298.15K
 709          T0  = 298.15D0
 710          T0T = T0/TEMP
 711          COEF= 1.0+LOG(T0T)-T0T
 712          XK1 = XK1 *EXP(  8.85*(T0T-1.0) + 25.140*COEF)
 713          XK21= XK21*EXP( 13.79*(T0T-1.0) -  5.393*COEF)
 714          XK22= XK22*EXP( -1.50*(T0T-1.0) + 26.920*COEF)
 715          XK4 = XK4 *EXP( 29.17*(T0T-1.0) + 16.830*COEF) !ISORR
 716 !CC         XK4 = XK4 *EXP( 29.47*(T0T-1.0) + 16.840*COEF) ! SEQUIL
 717          XK41= XK41*EXP( 29.17*(T0T-1.0) + 16.830*COEF)
 718          XK7 = XK7 *EXP( -2.65*(T0T-1.0) + 38.570*COEF)
 719          XK10= XK10*EXP(-74.38*(T0T-1.0) +  6.120*COEF) ! ISORR
 720 !CC         XK10= XK10*EXP(-75.11*(T0T-1.0) + 13.460*COEF) ! SEQUIL
 721          XK12= XK12*EXP( -2.87*(T0T-1.0) + 15.830*COEF)
 722          XK13= XK13*EXP( -5.19*(T0T-1.0) + 54.400*COEF)
 723          XKW = XKW *EXP(-22.52*(T0T-1.0) + 26.920*COEF)
 724       ENDIF
 725       XK2  = XK21*XK22
 726       XK42 = XK4/XK41
 727 !
 728 ! *** CALCULATE DELIQUESCENCE RELATIVE HUMIDITIES (UNICOMPONENT) ********
 729 !
 730       DRH2SO4  = ZERO
 731       DRNH42S4 = 0.7997D0
 732       DRNH4HS4 = 0.4000D0
 733       DRNH4NO3 = 0.6183D0
 734       DRLC     = 0.6900D0
 735       IF (INT(TEMP) .NE. 298) THEN
 736          T0       = 298.15D0
 737          TCF      = 1.0/TEMP - 1.0/T0
 738          DRNH4NO3 = DRNH4NO3*EXP(852.*TCF)
 739          DRNH42S4 = DRNH42S4*EXP( 80.*TCF)
 740          DRNH4HS4 = DRNH4HS4*EXP(384.*TCF)
 741          DRLC     = DRLC    *EXP(186.*TCF)
 742          DRNH4NO3 = MIN (DRNH4NO3,DRNH42S4) ! ADJUST FOR DRH CROSSOVER AT T<271K
 743       ENDIF
 744 !
 745 ! *** CALCULATE MUTUAL DELIQUESCENCE RELATIVE HUMIDITIES ****************
 746 !
 747       DRMLCAB = 0.3780D0              ! (NH4)3H(SO4)2 & NH4HSO4
 748       DRMLCAS = 0.6900D0              ! (NH4)3H(SO4)2 & (NH4)2SO4
 749       DRMASAN = 0.6000D0              ! (NH4)2SO4     & NH4NO3
 750 !CC      IF (INT(TEMP) .NE. 298) THEN    ! For the time being
 751 !CC         T0       = 298.15d0
 752 !CC         TCF      = 1.0/TEMP - 1.0/T0
 753 !CC         DRMLCAB  = DRMLCAB*EXP( 507.506*TCF)
 754 !CC         DRMLCAS  = DRMLCAS*EXP( 133.865*TCF)
 755 !CC         DRMASAN  = DRMASAN*EXP(1269.068*TCF)
 756 !CC      ENDIF
 757 !
 758 ! *** LIQUID PHASE ******************************************************
 759 !
 760       CHNO3  = ZERO
 761       CHCL   = ZERO
 762       CH2SO4 = ZERO
 763       COH    = ZERO
 764       WATER  = TINY
 765 !
 766       DO 20 I=1,NPAIR
 767          MOLALR(I)=ZERO
 768          GAMA(I)  =0.1
 769          GAMIN(I) =GREAT
 770          GAMOU(I) =GREAT
 771          M0(I)    =1d5
 772  20   CONTINUE
 773 !
 774       DO 30 I=1,NPAIR
 775          GAMA(I) = 0.1d0
 776  30   CONTINUE
 777 !
 778       DO 40 I=1,NIONS
 779          MOLAL(I)=ZERO
 780 40    CONTINUE
 781       COH = ZERO
 782 !
 783       DO 50 I=1,NGASAQ
 784          GASAQ(I)=ZERO
 785 50    CONTINUE
 786 !
 787 ! *** SOLID PHASE *******************************************************
 788 !
 789       CNH42S4= ZERO
 790       CNH4HS4= ZERO
 791       CNACL  = ZERO
 792       CNA2SO4= ZERO
 793       CNANO3 = ZERO
 794       CNH4NO3= ZERO
 795       CNH4CL = ZERO
 796       CNAHSO4= ZERO
 797       CLC    = ZERO
 798 !
 799 ! *** GAS PHASE *********************************************************
 800 !
 801       GNH3   = ZERO
 802       GHNO3  = ZERO
 803       GHCL   = ZERO
 804 !
 805 ! *** CALCULATE ZSR PARAMETERS ******************************************
 806 !
 807       IRH    = MIN (INT(RH*NZSR+0.5),NZSR)  ! Position in ZSR arrays
 808       IRH    = MAX (IRH, 1)
 809 !
 810 !      M0(01) = AWSC(IRH)      ! NACl
 811 !      IF (M0(01) .LT. 100.0) THEN
 812 !         IC = M0(01)
 813 !         CALL KMTAB(IC,298.0,     GI0,XX,XX,XX,XX,XX,XX,XX,XX,XX,XX,XX)
 814 !         CALL KMTAB(IC,SNGL(TEMP),GII,XX,XX,XX,XX,XX,XX,XX,XX,XX,XX,XX)
 815 !         M0(01) = M0(01)*EXP(LN10*(GI0-GII))
 816 !      ENDIF
 817 !
 818 !      M0(02) = AWSS(IRH)      ! (NA)2SO4
 819 !      IF (M0(02) .LT. 100.0) THEN
 820 !         IC = 3.0*M0(02)
 821 !         CALL KMTAB(IC,298.0,     XX,GI0,XX,XX,XX,XX,XX,XX,XX,XX,XX,XX)
 822 !         CALL KMTAB(IC,SNGL(TEMP),XX,GII,XX,XX,XX,XX,XX,XX,XX,XX,XX,XX)
 823 !         M0(02) = M0(02)*EXP(LN10*(GI0-GII))
 824 !      ENDIF
 825 !
 826 !      M0(03) = AWSN(IRH)      ! NANO3
 827 !      IF (M0(03) .LT. 100.0) THEN
 828 !         IC = M0(03)
 829 !         CALL KMTAB(IC,298.0,     XX,XX,GI0,XX,XX,XX,XX,XX,XX,XX,XX,XX)
 830 !         CALL KMTAB(IC,SNGL(TEMP),XX,XX,GII,XX,XX,XX,XX,XX,XX,XX,XX,XX)
 831 !         M0(03) = M0(03)*EXP(LN10*(GI0-GII))
 832 !      ENDIF
 833 !
 834       M0(04) = AWAS(IRH)      ! (NH4)2SO4
 835 !      IF (M0(04) .LT. 100.0) THEN
 836 !         IC = 3.0*M0(04)
 837 !         CALL KMTAB(IC,298.0,     XX,XX,XX,GI0,XX,XX,XX,XX,XX,XX,XX,XX)
 838 !         CALL KMTAB(IC,SNGL(TEMP),XX,XX,XX,GII,XX,XX,XX,XX,XX,XX,XX,XX)
 839 !         M0(04) = M0(04)*EXP(LN10*(GI0-GII))
 840 !      ENDIF
 841 !
 842       M0(05) = AWAN(IRH)      ! NH4NO3
 843 !      IF (M0(05) .LT. 100.0) THEN
 844 !         IC     = M0(05)
 845 !         CALL KMTAB(IC,298.0,     XX,XX,XX,XX,GI0,XX,XX,XX,XX,XX,XX,XX)
 846 !         CALL KMTAB(IC,SNGL(TEMP),XX,XX,XX,XX,GII,XX,XX,XX,XX,XX,XX,XX)
 847 !         M0(05) = M0(05)*EXP(LN10*(GI0-GII))
 848 !      ENDIF
 849 !
 850 !      M0(06) = AWAC(IRH)      ! NH4CL
 851 !      IF (M0(06) .LT. 100.0) THEN
 852 !         IC = M0(06)
 853 !         CALL KMTAB(IC,298.0,     XX,XX,XX,XX,XX,GI0,XX,XX,XX,XX,XX,XX)
 854 !         CALL KMTAB(IC,SNGL(TEMP),XX,XX,XX,XX,XX,GII,XX,XX,XX,XX,XX,XX)
 855 !         M0(06) = M0(06)*EXP(LN10*(GI0-GII))
 856 !      ENDIF
 857 !
 858       M0(07) = AWSA(IRH)      ! 2H-SO4
 859 !      IF (M0(07) .LT. 100.0) THEN
 860 !         IC = 3.0*M0(07)
 861 !         CALL KMTAB(IC,298.0,     XX,XX,XX,XX,XX,XX,GI0,XX,XX,XX,XX,XX)
 862 !         CALL KMTAB(IC,SNGL(TEMP),XX,XX,XX,XX,XX,XX,GII,XX,XX,XX,XX,XX)
 863 !         M0(07) = M0(07)*EXP(LN10*(GI0-GII))
 864 !      ENDIF
 865 !
 866       M0(08) = AWSA(IRH)      ! H-HSO4
 867 !CC      IF (M0(08) .LT. 100.0) THEN     ! These are redundant, because M0(8) is not used
 868 !CC         IC = M0(08)
 869 !CC         CALL KMTAB(IC,298.0,     XX,XX,XX,XX,XX,XX,XX,GI0,XX,XX,XX,XX)
 870 !CC         CALL KMTAB(IC,SNGL(TEMP),XX,XX,XX,XX,XX,XX,XX,GI0,XX,XX,XX,XX)
 871 !CCCCC         CALL KMTAB(IC,SNGL(TEMP),XX,XX,XX,XX,XX,XX,XX,GII,XX,XX,XX,XX)
 872 !CC         M0(08) = M0(08)*EXP(LN10*(GI0-GII))
 873 !CC      ENDIF
 874 !
 875       M0(09) = AWAB(IRH)      ! NH4HSO4
 876 !      IF (M0(09) .LT. 100.0) THEN
 877 !         IC = M0(09)
 878 !         CALL KMTAB(IC,298.0,     XX,XX,XX,XX,XX,XX,XX,XX,GI0,XX,XX,XX)
 879 !         CALL KMTAB(IC,SNGL(TEMP),XX,XX,XX,XX,XX,XX,XX,XX,GII,XX,XX,XX)
 880 !         M0(09) = M0(09)*EXP(LN10*(GI0-GII))
 881 !      ENDIF
 882 !
 883 !      M0(12) = AWSB(IRH)      ! NAHSO4
 884 !      IF (M0(12) .LT. 100.0) THEN
 885 !         IC = M0(12)
 886 !         CALL KMTAB(IC,298.0,     XX,XX,XX,XX,XX,XX,XX,XX,XX,XX,XX,GI0)
 887 !         CALL KMTAB(IC,SNGL(TEMP),XX,XX,XX,XX,XX,XX,XX,XX,XX,XX,XX,GII)
 888 !         M0(12) = M0(12)*EXP(LN10*(GI0-GII))
 889 !      ENDIF
 890 !
 891       M0(13) = AWLC(IRH)      ! (NH4)3H(SO4)2
 892 !      IF (M0(13) .LT. 100.0) THEN
 893 !         IC     = 4.0*M0(13)
 894 !         CALL KMTAB(IC,298.0,     XX,XX,XX,GI0,XX,XX,XX,XX,GII,XX,XX,XX)
 895 !         G130   = 0.2*(3.0*GI0+2.0*GII)
 896 !         CALL KMTAB(IC,SNGL(TEMP),XX,XX,XX,GI0,XX,XX,XX,XX,GII,XX,XX,XX)
 897 !         G13I   = 0.2*(3.0*GI0+2.0*GII)
 898 !         M0(13) = M0(13)*EXP(LN10*SNGL(G130-G13I))
 899 !      ENDIF
 900 !
 901 ! *** OTHER INITIALIZATIONS *********************************************
 902 !
 903       ICLACT  = 0
 904       CALAOU  = .TRUE.
 905       CALAIN  = .TRUE.
 906       FRST    = .TRUE.
 907       SCASE   = '??'
 908       SULRATW = 2.D0
 909       SODRAT  = ZERO
 910       NOFER   = 0
 911       STKOFL  =.FALSE.
 912       DO 60 I=1,NERRMX
 913          ERRSTK(I) =-999
 914          ERRMSG(I) = 'MESSAGE N/A'
 915    60 CONTINUE
 916 !
 917 ! *** END OF SUBROUTINE INIT2 *******************************************
 918 !
 919       END
 920
 921
 922
 923
 924
 925 !=======================================================================
 926 !
 927 ! *** ISORROPIA CODE
 928 ! *** SUBROUTINE ISOINIT3
 929 ! *** THIS SUBROUTINE INITIALIZES ALL GLOBAL VARIABLES FOR AMMONIUM,
 930 !     SODIUM, CHLORIDE, NITRATE, SULFATE AEROSOL SYSTEMS (SUBROUTINE
 931 !     ISRP3)
 932 !
 933 ! *** COPYRIGHT 1996-2006, UNIVERSITY OF MIAMI, CARNEGIE MELLON UNIVERSITY,
 934 ! *** GEORGIA INSTITUTE OF TECHNOLOGY
 935 ! *** WRITTEN BY ATHANASIOS NENES
 936 ! *** UPDATED BY CHRISTOS FOUNTOUKIS
 937 !
 938 !=======================================================================
 939 !
 940       SUBROUTINE ISOINIT3 (WI, RHI, TEMPI)
 941       USE ISRPIA
 942       IMPLICIT DOUBLE PRECISION (A-H, O-Z)
 943 !      INCLUDE 'ISRPIA.INC'
 944       DIMENSION WI(NCOMP)
 945       REAL      IC,GII,GI0,XX,LN10
 946       PARAMETER (LN10=2.3025851)
 947 !
 948 ! *** SAVE INPUT VARIABLES IN COMMON BLOCK ******************************
 949 !
 950       IF (IPROB.EQ.0) THEN                 ! FORWARD CALCULATION
 951          DO 10 I=1,NCOMP
 952             W(I) = MAX(WI(I), TINY)
 953 10       CONTINUE
 954       ELSE
 955          DO 15 I=1,NCOMP                   ! REVERSE CALCULATION
 956             WAER(I) = MAX(WI(I), TINY)
 957             W(I)    = ZERO
 958 15       CONTINUE
 959       ENDIF
 960       RH      = RHI
 961       TEMP    = TEMPI
 962 !
 963 ! *** CALCULATE EQUILIBRIUM CONSTANTS ***********************************
 964 !
 965       XK1  = 1.015D-2  ! HSO4(aq)         <==> H(aq)     + SO4(aq)
 966       XK21 = 57.639D0  ! NH3(g)           <==> NH3(aq)
 967       XK22 = 1.805D-5  ! NH3(aq)          <==> NH4(aq)   + OH(aq)
 968       XK3  = 1.971D6   ! HCL(g)           <==> H(aq)     + CL(aq)
 969       XK31 = 2.500e3   ! HCL(g)           <==> HCL(aq)
 970       XK4  = 2.511e6   ! HNO3(g)          <==> H(aq)     + NO3(aq) ! ISORR
 971 !CC      XK4  = 3.638e6   ! HNO3(g)          <==> H(aq)     + NO3(aq) ! SEQUIL
 972       XK41 = 2.100e5   ! HNO3(g)          <==> HNO3(aq)
 973       XK5  = 0.4799D0  ! NA2SO4(s)        <==> 2*NA(aq)  + SO4(aq)
 974       XK6  = 1.086D-16 ! NH4CL(s)         <==> NH3(g)    + HCL(g)
 975       XK7  = 1.817D0   ! (NH4)2SO4(s)     <==> 2*NH4(aq) + SO4(aq)
 976       XK8  = 37.661D0  ! NACL(s)          <==> NA(aq)    + CL(aq)
 977       XK10 = 5.746D-17 ! NH4NO3(s)        <==> NH3(g)    + HNO3(g) ! ISORR
 978 !CC      XK10 = 2.985e-17 ! NH4NO3(s)        <==> NH3(g)    + HNO3(g) ! SEQUIL
 979       XK11 = 2.413D4   ! NAHSO4(s)        <==> NA(aq)    + HSO4(aq)
 980       XK12 = 1.382D2   ! NH4HSO4(s)       <==> NH4(aq)   + HSO4(aq)
 981       XK13 = 29.268D0  ! (NH4)3H(SO4)2(s) <==> 3*NH4(aq) + HSO4(aq) + SO4(aq)
 982       XK14 = 22.05D0   ! NH4CL(s)         <==> NH4(aq)   + CL(aq)
 983       XKW  = 1.010D-14 ! H2O              <==> H(aq)     + OH(aq)
 984       XK9  = 11.977D0  ! NANO3(s)         <==> NA(aq)    + NO3(aq)
 985 !
 986       IF (INT(TEMP) .NE. 298) THEN   ! FOR T != 298K or 298.15K
 987          T0  = 298.15D0
 988          T0T = T0/TEMP
 989          COEF= 1.0+LOG(T0T)-T0T
 990          XK1 = XK1 *EXP(  8.85*(T0T-1.0) + 25.140*COEF)
 991          XK21= XK21*EXP( 13.79*(T0T-1.0) -  5.393*COEF)
 992          XK22= XK22*EXP( -1.50*(T0T-1.0) + 26.920*COEF)
 993          XK3 = XK3 *EXP( 30.20*(T0T-1.0) + 19.910*COEF)
 994          XK31= XK31*EXP( 30.20*(T0T-1.0) + 19.910*COEF)
 995          XK4 = XK4 *EXP( 29.17*(T0T-1.0) + 16.830*COEF) !ISORR
 996 !CC         XK4 = XK4 *EXP( 29.47*(T0T-1.0) + 16.840*COEF) ! SEQUIL
 997          XK41= XK41*EXP( 29.17*(T0T-1.0) + 16.830*COEF)
 998          XK5 = XK5 *EXP(  0.98*(T0T-1.0) + 39.500*COEF)
 999          XK6 = XK6 *EXP(-71.00*(T0T-1.0) +  2.400*COEF)
1000          XK7 = XK7 *EXP( -2.65*(T0T-1.0) + 38.570*COEF)
1001          XK8 = XK8 *EXP( -1.56*(T0T-1.0) + 16.900*COEF)
1002          XK9 = XK9 *EXP( -8.22*(T0T-1.0) + 16.010*COEF)
1003          XK10= XK10*EXP(-74.38*(T0T-1.0) +  6.120*COEF) ! ISORR
1004 !CC         XK10= XK10*EXP(-75.11*(T0T-1.0) + 13.460*COEF) ! SEQUIL
1005          XK11= XK11*EXP(  0.79*(T0T-1.0) + 14.746*COEF)
1006          XK12= XK12*EXP( -2.87*(T0T-1.0) + 15.830*COEF)
1007          XK13= XK13*EXP( -5.19*(T0T-1.0) + 54.400*COEF)
1008          XK14= XK14*EXP( 24.55*(T0T-1.0) + 16.900*COEF)
1009          XKW = XKW *EXP(-22.52*(T0T-1.0) + 26.920*COEF)
1010       ENDIF
1011       XK2  = XK21*XK22
1012       XK42 = XK4/XK41
1013       XK32 = XK3/XK31
1014 !
1015 ! *** CALCULATE DELIQUESCENCE RELATIVE HUMIDITIES (UNICOMPONENT) ********
1016 !
1017       DRH2SO4  = ZERO
1018       DRNH42S4 = 0.7997D0
1019       DRNH4HS4 = 0.4000D0
1020       DRLC     = 0.6900D0
1021       DRNACL   = 0.7528D0
1022       DRNANO3  = 0.7379D0
1023       DRNH4CL  = 0.7710D0
1024       DRNH4NO3 = 0.6183D0
1025       DRNA2SO4 = 0.9300D0
1026       DRNAHSO4 = 0.5200D0
1027       IF (INT(TEMP) .NE. 298) THEN
1028          T0       = 298.15D0
1029          TCF      = 1.0/TEMP - 1.0/T0
1030          DRNACL   = DRNACL  *EXP( 25.*TCF)
1031          DRNANO3  = DRNANO3 *EXP(304.*TCF)
1032          DRNA2SO4 = DRNA2SO4*EXP( 80.*TCF)
1033          DRNH4NO3 = DRNH4NO3*EXP(852.*TCF)
1034          DRNH42S4 = DRNH42S4*EXP( 80.*TCF)
1035          DRNH4HS4 = DRNH4HS4*EXP(384.*TCF)
1036          DRLC     = DRLC    *EXP(186.*TCF)
1037          DRNH4CL  = DRNH4Cl *EXP(239.*TCF)
1038          DRNAHSO4 = DRNAHSO4*EXP(-45.*TCF)
1039 !
1040 ! *** ADJUST FOR DRH "CROSSOVER" AT LOW TEMPERATURES
1041 !
1042          DRNH4NO3  = MIN (DRNH4NO3, DRNH4CL, DRNH42S4, DRNANO3, DRNACL)
1043          DRNANO3   = MIN (DRNANO3, DRNACL)
1044          DRNH4CL   = MIN (DRNH4Cl, DRNH42S4)
1045 !
1046       ENDIF
1047 !
1048 ! *** CALCULATE MUTUAL DELIQUESCENCE RELATIVE HUMIDITIES ****************
1049 !
1050       DRMLCAB = 0.378D0    ! (NH4)3H(SO4)2 & NH4HSO4
1051       DRMLCAS = 0.690D0    ! (NH4)3H(SO4)2 & (NH4)2SO4
1052       DRMASAN = 0.600D0    ! (NH4)2SO4     & NH4NO3
1053       DRMG1   = 0.460D0    ! (NH4)2SO4, NH4NO3, NA2SO4, NH4CL
1054       DRMG2   = 0.691D0    ! (NH4)2SO4, NA2SO4, NH4CL
1055       DRMG3   = 0.697D0    ! (NH4)2SO4, NA2SO4
1056       DRMH1   = 0.240D0    ! NA2SO4, NANO3, NACL, NH4NO3, NH4CL
1057       DRMH2   = 0.596D0    ! NA2SO4, NANO3, NACL, NH4CL
1058       DRMI1   = 0.240D0    ! LC, NAHSO4, NH4HSO4, NA2SO4, (NH4)2SO4
1059       DRMI2   = 0.363D0    ! LC, NAHSO4, NA2SO4, (NH4)2SO4  - NO DATA -
1060       DRMI3   = 0.610D0    ! LC, NA2SO4, (NH4)2SO4
1061       DRMQ1   = 0.494D0    ! (NH4)2SO4, NH4NO3, NA2SO4
1062       DRMR1   = 0.663D0    ! NA2SO4, NANO3, NACL
1063       DRMR2   = 0.735D0    ! NA2SO4, NACL
1064       DRMR3   = 0.673D0    ! NANO3, NACL
1065       DRMR4   = 0.694D0    ! NA2SO4, NACL, NH4CL
1066       DRMR5   = 0.731D0    ! NA2SO4, NH4CL
1067       DRMR6   = 0.596D0    ! NA2SO4, NANO3, NH4CL
1068       DRMR7   = 0.380D0    ! NA2SO4, NANO3, NACL, NH4NO3
1069       DRMR8   = 0.380D0    ! NA2SO4, NACL, NH4NO3
1070       DRMR9   = 0.494D0    ! NA2SO4, NH4NO3
1071       DRMR10  = 0.476D0    ! NA2SO4, NANO3, NH4NO3
1072       DRMR11  = 0.340D0    ! NA2SO4, NACL, NH4NO3, NH4CL
1073       DRMR12  = 0.460D0    ! NA2SO4, NH4NO3, NH4CL
1074       DRMR13  = 0.438D0    ! NA2SO4, NANO3, NH4NO3, NH4CL
1075 !CC      IF (INT(TEMP) .NE. 298) THEN
1076 !CC         T0       = 298.15d0
1077 !CC         TCF      = 1.0/TEMP - 1.0/T0
1078 !CC         DRMLCAB  = DRMLCAB*EXP( 507.506*TCF)
1079 !CC         DRMLCAS  = DRMLCAS*EXP( 133.865*TCF)
1080 !CC         DRMASAN  = DRMASAN*EXP(1269.068*TCF)
1081 !CC         DRMG1    = DRMG1  *EXP( 572.207*TCF)
1082 !CC         DRMG2    = DRMG2  *EXP(  58.166*TCF)
1083 !CC         DRMG3    = DRMG3  *EXP(  22.253*TCF)
1084 !CC         DRMH1    = DRMH1  *EXP(2116.542*TCF)
1085 !CC         DRMH2    = DRMH2  *EXP( 650.549*TCF)
1086 !CC         DRMI1    = DRMI1  *EXP( 565.743*TCF)
1087 !CC         DRMI2    = DRMI2  *EXP(  91.745*TCF)
1088 !CC         DRMI3    = DRMI3  *EXP( 161.272*TCF)
1089 !CC         DRMQ1    = DRMQ1  *EXP(1616.621*TCF)
1090 !CC         DRMR1    = DRMR1  *EXP( 292.564*TCF)
1091 !CC         DRMR2    = DRMR2  *EXP(  14.587*TCF)
1092 !CC         DRMR3    = DRMR3  *EXP( 307.907*TCF)
1093 !CC         DRMR4    = DRMR4  *EXP(  97.605*TCF)
1094 !CC         DRMR5    = DRMR5  *EXP(  98.523*TCF)
1095 !CC         DRMR6    = DRMR6  *EXP( 465.500*TCF)
1096 !CC         DRMR7    = DRMR7  *EXP( 324.425*TCF)
1097 !CC         DRMR8    = DRMR8  *EXP(2660.184*TCF)
1098 !CC         DRMR9    = DRMR9  *EXP(1617.178*TCF)
1099 !CC         DRMR10   = DRMR10 *EXP(1745.226*TCF)
1100 !CC         DRMR11   = DRMR11 *EXP(3691.328*TCF)
1101 !CC         DRMR12   = DRMR12 *EXP(1836.842*TCF)
1102 !CC         DRMR13   = DRMR13 *EXP(1967.938*TCF)
1103 !CC      ENDIF
1104 !
1105 ! *** LIQUID PHASE ******************************************************
1106 !
1107       CHNO3  = ZERO
1108       CHCL   = ZERO
1109       CH2SO4 = ZERO
1110       COH    = ZERO
1111       WATER  = TINY
1112 !
1113       DO 20 I=1,NPAIR
1114          MOLALR(I)=ZERO
1115          GAMA(I)  =0.1
1116          GAMIN(I) =GREAT
1117          GAMOU(I) =GREAT
1118          M0(I)    =1d5
1119  20   CONTINUE
1120 !
1121       DO 30 I=1,NPAIR
1122          GAMA(I) = 0.1d0
1123  30   CONTINUE
1124 !
1125       DO 40 I=1,NIONS
1126          MOLAL(I)=ZERO
1127 40    CONTINUE
1128       COH = ZERO
1129 !
1130       DO 50 I=1,NGASAQ
1131          GASAQ(I)=ZERO
1132 50    CONTINUE
1133 !
1134 ! *** SOLID PHASE *******************************************************
1135 !
1136       CNH42S4= ZERO
1137       CNH4HS4= ZERO
1138       CNACL  = ZERO
1139       CNA2SO4= ZERO
1140       CNANO3 = ZERO
1141       CNH4NO3= ZERO
1142       CNH4CL = ZERO
1143       CNAHSO4= ZERO
1144       CLC    = ZERO
1145 !
1146 ! *** GAS PHASE *********************************************************
1147 !
1148       GNH3   = ZERO
1149       GHNO3  = ZERO
1150       GHCL   = ZERO
1151 !
1152 ! *** CALCULATE ZSR PARAMETERS ******************************************
1153 !
1154       IRH    = MIN (INT(RH*NZSR+0.5),NZSR)  ! Position in ZSR arrays
1155       IRH    = MAX (IRH, 1)
1156 !
1157       M0(01) = AWSC(IRH)      ! NACl
1158 !      IF (M0(01) .LT. 100.0) THEN
1159 !         IC = M0(01)
1160 !         CALL KMTAB(IC,298.0,     GI0,XX,XX,XX,XX,XX,XX,XX,XX,XX,XX,XX)
1161 !         CALL KMTAB(IC,SNGL(TEMP),GII,XX,XX,XX,XX,XX,XX,XX,XX,XX,XX,XX)
1162 !         M0(01) = M0(01)*EXP(LN10*(GI0-GII))
1163 !      ENDIF
1164 !
1165       M0(02) = AWSS(IRH)      ! (NA)2SO4
1166 !      IF (M0(02) .LT. 100.0) THEN
1167 !         IC = 3.0*M0(02)
1168 !         CALL KMTAB(IC,298.0,     XX,GI0,XX,XX,XX,XX,XX,XX,XX,XX,XX,XX)
1169 !         CALL KMTAB(IC,SNGL(TEMP),XX,GII,XX,XX,XX,XX,XX,XX,XX,XX,XX,XX)
1170 !         M0(02) = M0(02)*EXP(LN10*(GI0-GII))
1171 !      ENDIF
1172 !
1173       M0(03) = AWSN(IRH)      ! NANO3
1174 !      IF (M0(03) .LT. 100.0) THEN
1175 !         IC = M0(03)
1176 !         CALL KMTAB(IC,298.0,     XX,XX,GI0,XX,XX,XX,XX,XX,XX,XX,XX,XX)
1177 !         CALL KMTAB(IC,SNGL(TEMP),XX,XX,GII,XX,XX,XX,XX,XX,XX,XX,XX,XX)
1178 !         M0(03) = M0(03)*EXP(LN10*(GI0-GII))
1179 !      ENDIF
1180 !
1181       M0(04) = AWAS(IRH)      ! (NH4)2SO4
1182 !      IF (M0(04) .LT. 100.0) THEN
1183 !         IC = 3.0*M0(04)
1184 !         CALL KMTAB(IC,298.0,     XX,XX,XX,GI0,XX,XX,XX,XX,XX,XX,XX,XX)
1185 !         CALL KMTAB(IC,SNGL(TEMP),XX,XX,XX,GII,XX,XX,XX,XX,XX,XX,XX,XX)
1186 !         M0(04) = M0(04)*EXP(LN10*(GI0-GII))
1187 !      ENDIF
1188 !
1189       M0(05) = AWAN(IRH)      ! NH4NO3
1190 !      IF (M0(05) .LT. 100.0) THEN
1191 !        IC     = M0(05)
1192 !        CALL KMTAB(IC,298.0,     XX,XX,XX,XX,GI0,XX,XX,XX,XX,XX,XX,XX)
1193 !        CALL KMTAB(IC,SNGL(TEMP),XX,XX,XX,XX,GII,XX,XX,XX,XX,XX,XX,XX)
1194 !         M0(05) = M0(05)*EXP(LN10*(GI0-GII))
1195 !      ENDIF
1196 !
1197       M0(06) = AWAC(IRH)      ! NH4CL
1198 !      IF (M0(06) .LT. 100.0) THEN
1199 !         IC = M0(06)
1200 !         CALL KMTAB(IC,298.0,     XX,XX,XX,XX,XX,GI0,XX,XX,XX,XX,XX,XX)
1201 !         CALL KMTAB(IC,SNGL(TEMP),XX,XX,XX,XX,XX,GII,XX,XX,XX,XX,XX,XX)
1202 !         M0(06) = M0(06)*EXP(LN10*(GI0-GII))
1203 !      ENDIF
1204 !
1205       M0(07) = AWSA(IRH)      ! 2H-SO4
1206 !      IF (M0(07) .LT. 100.0) THEN
1207 !         IC = 3.0*M0(07)
1208 !         CALL KMTAB(IC,298.0,     XX,XX,XX,XX,XX,XX,GI0,XX,XX,XX,XX,XX)
1209 !         CALL KMTAB(IC,SNGL(TEMP),XX,XX,XX,XX,XX,XX,GII,XX,XX,XX,XX,XX)
1210 !         M0(07) = M0(07)*EXP(LN10*(GI0-GII))
1211 !      ENDIF
1212 !
1213       M0(08) = AWSA(IRH)      ! H-HSO4
1214 !CC      IF (M0(08) .LT. 100.0) THEN     ! These are redundant, because M0(8) is not used
1215 !CC         IC = M0(08)
1216 !CC         CALL KMTAB(IC,298.0,     XX,XX,XX,XX,XX,XX,XX,GI0,XX,XX,XX,XX)
1217 !CC         CALL KMTAB(IC,SNGL(TEMP),XX,XX,XX,XX,XX,XX,XX,GI0,XX,XX,XX,XX)
1218 !CCCCC         CALL KMTAB(IC,SNGL(TEMP),XX,XX,XX,XX,XX,XX,XX,GII,XX,XX,XX,XX)
1219 !CC         M0(08) = M0(08)*EXP(LN10*(GI0-GII))
1220 !CC      ENDIF
1221 !
1222       M0(09) = AWAB(IRH)      ! NH4HSO4
1223 !      IF (M0(09) .LT. 100.0) THEN
1224 !         IC = M0(09)
1225 !         CALL KMTAB(IC,298.0,     XX,XX,XX,XX,XX,XX,XX,XX,GI0,XX,XX,XX)
1226 !         CALL KMTAB(IC,SNGL(TEMP),XX,XX,XX,XX,XX,XX,XX,XX,GII,XX,XX,XX)
1227 !         M0(09) = M0(09)*EXP(LN10*(GI0-GII))
1228 !      ENDIF
1229 !
1230       M0(12) = AWSB(IRH)      ! NAHSO4
1231 !      IF (M0(12) .LT. 100.0) THEN
1232 !         IC = M0(12)
1233 !         CALL KMTAB(IC,298.0,     XX,XX,XX,XX,XX,XX,XX,XX,XX,XX,XX,GI0)
1234 !         CALL KMTAB(IC,SNGL(TEMP),XX,XX,XX,XX,XX,XX,XX,XX,XX,XX,XX,GII)
1235 !         M0(12) = M0(12)*EXP(LN10*(GI0-GII))
1236 !      ENDIF
1237 !
1238       M0(13) = AWLC(IRH)      ! (NH4)3H(SO4)2
1239 !      IF (M0(13) .LT. 100.0) THEN
1240 !         IC     = 4.0*M0(13)
1241 !         CALL KMTAB(IC,298.0,     XX,XX,XX,GI0,XX,XX,XX,XX,GII,XX,XX,XX)
1242 !         G130   = 0.2*(3.0*GI0+2.0*GII)
1243 !         CALL KMTAB(IC,SNGL(TEMP),XX,XX,XX,GI0,XX,XX,XX,XX,GII,XX,XX,XX)
1244 !         G13I   = 0.2*(3.0*GI0+2.0*GII)
1245 !         M0(13) = M0(13)*EXP(LN10*SNGL(G130-G13I))
1246 !      ENDIF
1247 !
1248 ! *** OTHER INITIALIZATIONS *********************************************
1249 !
1250       ICLACT  = 0
1251       CALAOU  = .TRUE.
1252       CALAIN  = .TRUE.
1253       FRST    = .TRUE.
1254       SCASE   = '??'
1255       SULRATW = 2.D0
1256       NOFER   = 0
1257       STKOFL  =.FALSE.
1258       DO 60 I=1,NERRMX
1259          ERRSTK(I) =-999
1260          ERRMSG(I) = 'MESSAGE N/A'
1261    60 CONTINUE
1262 !
1263 ! *** END OF SUBROUTINE ISOINIT3 *******************************************
1264 !
1265       END
1266
1267 !=======================================================================
1268 !
1269 ! *** ISORROPIA CODE
1270 ! *** SUBROUTINE ADJUST
1271 ! *** ADJUSTS FOR MASS BALANCE BETWEEN VOLATILE SPECIES AND SULFATE
1272 !     FIRST CALCULATE THE EXCESS OF EACH PRECURSOR, AND IF IT EXISTS, THEN
1273 !     ADJUST SEQUENTIALY AEROSOL PHASE SPECIES WHICH CONTAIN THE EXCESS
1274 !     PRECURSOR.
1275 !
1276 ! *** COPYRIGHT 1996-2006, UNIVERSITY OF MIAMI, CARNEGIE MELLON UNIVERSITY,
1277 ! *** GEORGIA INSTITUTE OF TECHNOLOGY
1278 ! *** WRITTEN BY ATHANASIOS NENES
1279 ! *** UPDATED BY CHRISTOS FOUNTOUKIS
1280 !
1281 !=======================================================================
1282 !
1283       SUBROUTINE ADJUST (WI)
1284       USE ISRPIA
1285       IMPLICIT DOUBLE PRECISION (A-H, O-Z)
1286 !      INCLUDE 'ISRPIA.INC'
1287       DOUBLE PRECISION WI(*)
1288 !
1289 ! *** FOR AMMONIUM *****************************************************
1290 !
1291       IF (IPROB.EQ.0) THEN         ! Calculate excess (solution - input)
1292          EXNH4 = GNH3 + MOLAL(3) + CNH4CL + CNH4NO3 + CNH4HS4   &
1293                       + 2D0*CNH42S4       + 3D0*CLC   &
1294                 -WI(3)
1295       ELSE
1296          EXNH4 = MOLAL(3) + CNH4CL + CNH4NO3 + CNH4HS4 + 2D0*CNH42S4   &
1297                           + 3D0*CLC   &
1298                 -WI(3)
1299
1300       ENDIF
1301       EXNH4 = MAX(EXNH4,ZERO)
1302       IF (EXNH4.LT.TINY) GOTO 20    ! No excess NH4, go to next precursor
1303 !
1304       IF (MOLAL(3).GT.EXNH4) THEN   ! Adjust aqueous phase NH4
1305          MOLAL(3) = MOLAL(3) - EXNH4
1306          GOTO 20
1307       ELSE
1308          EXNH4    = EXNH4 - MOLAL(3)
1309          MOLAL(3) = ZERO
1310       ENDIF
1311 !
1312       IF (CNH4CL.GT.EXNH4) THEN     ! Adjust NH4Cl(s)
1313          CNH4CL   = CNH4CL - EXNH4  ! more solid than excess
1314          GHCL     = GHCL   + EXNH4  ! evaporate Cl to gas phase
1315          GOTO 20
1316       ELSE                          ! less solid than excess
1317          GHCL     = GHCL   + CNH4CL ! evaporate into gas phase
1318          EXNH4    = EXNH4  - CNH4CL ! reduce excess
1319          CNH4CL   = ZERO            ! zero salt concentration
1320       ENDIF
1321 !
1322       IF (CNH4NO3.GT.EXNH4) THEN    ! Adjust NH4NO3(s)
1323          CNH4NO3  = CNH4NO3- EXNH4  ! more solid than excess
1324          GHNO3    = GHNO3  + EXNH4  ! evaporate NO3 to gas phase
1325          GOTO 20
1326       ELSE                          ! less solid than excess
1327          GHNO3    = GHNO3  + CNH4NO3! evaporate into gas phase
1328          EXNH4    = EXNH4  - CNH4NO3! reduce excess
1329          CNH4NO3  = ZERO            ! zero salt concentration
1330       ENDIF
1331 !
1332       IF (CLC.GT.3d0*EXNH4) THEN    ! Adjust (NH4)3H(SO4)2(s)
1333          CLC      = CLC - EXNH4/3d0 ! more solid than excess
1334          GOTO 20
1335       ELSE                          ! less solid than excess
1336          EXNH4    = EXNH4 - 3d0*CLC ! reduce excess
1337          CLC      = ZERO            ! zero salt concentration
1338       ENDIF
1339 !
1340       IF (CNH4HS4.GT.EXNH4) THEN    ! Adjust NH4HSO4(s)
1341          CNH4HS4  = CNH4HS4- EXNH4  ! more solid than excess
1342          GOTO 20
1343       ELSE                          ! less solid than excess
1344          EXNH4    = EXNH4  - CNH4HS4! reduce excess
1345          CNH4HS4  = ZERO            ! zero salt concentration
1346       ENDIF
1347 !
1348       IF (CNH42S4.GT.EXNH4) THEN    ! Adjust (NH4)2SO4(s)
1349          CNH42S4  = CNH42S4- EXNH4  ! more solid than excess
1350          GOTO 20
1351       ELSE                          ! less solid than excess
1352          EXNH4    = EXNH4  - CNH42S4! reduce excess
1353          CNH42S4  = ZERO            ! zero salt concentration
1354       ENDIF
1355 !
1356 ! *** FOR NITRATE ******************************************************
1357 !
1358  20   IF (IPROB.EQ.0) THEN         ! Calculate excess (solution - input)
1359          EXNO3 = GHNO3 + MOLAL(7) + CNH4NO3   &
1360                 -WI(4)
1361       ELSE
1362          EXNO3 = MOLAL(7) + CNH4NO3   &
1363                 -WI(4)
1364       ENDIF
1365       EXNO3 = MAX(EXNO3,ZERO)
1366       IF (EXNO3.LT.TINY) GOTO 30    ! No excess NO3, go to next precursor
1367 !
1368       IF (MOLAL(7).GT.EXNO3) THEN   ! Adjust aqueous phase NO3
1369          MOLAL(7) = MOLAL(7) - EXNO3
1370          GOTO 30
1371       ELSE
1372          EXNO3    = EXNO3 - MOLAL(7)
1373          MOLAL(7) = ZERO
1374       ENDIF
1375 !
1376       IF (CNH4NO3.GT.EXNO3) THEN    ! Adjust NH4NO3(s)
1377          CNH4NO3  = CNH4NO3- EXNO3  ! more solid than excess
1378          GNH3     = GNH3   + EXNO3  ! evaporate NO3 to gas phase
1379          GOTO 30
1380       ELSE                          ! less solid than excess
1381          GNH3     = GNH3   + CNH4NO3! evaporate into gas phase
1382          EXNO3    = EXNO3  - CNH4NO3! reduce excess
1383          CNH4NO3  = ZERO            ! zero salt concentration
1384       ENDIF
1385 !
1386 ! *** FOR CHLORIDE *****************************************************
1387 !
1388  30   IF (IPROB.EQ.0) THEN         ! Calculate excess (solution - input)
1389          EXCl = GHCL + MOLAL(4) + CNH4CL   &
1390                -WI(5)
1391       ELSE
1392          EXCl = MOLAL(4) + CNH4CL   &
1393                -WI(5)
1394       ENDIF
1395       EXCl = MAX(EXCl,ZERO)
1396       IF (EXCl.LT.TINY) GOTO 40    ! No excess Cl, go to next precursor
1397 !
1398       IF (MOLAL(4).GT.EXCL) THEN   ! Adjust aqueous phase Cl
1399          MOLAL(4) = MOLAL(4) - EXCL
1400          GOTO 40
1401       ELSE
1402          EXCL     = EXCL - MOLAL(4)
1403          MOLAL(4) = ZERO
1404       ENDIF
1405 !
1406       IF (CNH4CL.GT.EXCL) THEN      ! Adjust NH4Cl(s)
1407          CNH4CL   = CNH4CL - EXCL   ! more solid than excess
1408          GHCL     = GHCL   + EXCL   ! evaporate Cl to gas phase
1409          GOTO 40
1410       ELSE                          ! less solid than excess
1411          GHCL     = GHCL   + CNH4CL ! evaporate into gas phase
1412          EXCL     = EXCL   - CNH4CL ! reduce excess
1413          CNH4CL   = ZERO            ! zero salt concentration
1414       ENDIF
1415 !
1416 ! *** FOR SULFATE ******************************************************
1417 !
1418  40   EXS4 = MOLAL(5) + MOLAL(6) + 2.d0*CLC + CNH42S4 + CNH4HS4 +   &
1419              CNA2SO4  + CNAHSO4 - WI(2)
1420       EXS4 = MAX(EXS4,ZERO)        ! Calculate excess (solution - input)
1421       IF (EXS4.LT.TINY) GOTO 50    ! No excess SO4, return
1422 !
1423       IF (MOLAL(6).GT.EXS4) THEN   ! Adjust aqueous phase HSO4
1424          MOLAL(6) = MOLAL(6) - EXS4
1425          GOTO 50
1426       ELSE
1427          EXS4     = EXS4 - MOLAL(6)
1428          MOLAL(6) = ZERO
1429       ENDIF
1430 !
1431       IF (MOLAL(5).GT.EXS4) THEN   ! Adjust aqueous phase SO4
1432          MOLAL(5) = MOLAL(5) - EXS4
1433          GOTO 50
1434       ELSE
1435          EXS4     = EXS4 - MOLAL(5)
1436          MOLAL(5) = ZERO
1437       ENDIF
1438 !
1439       IF (CLC.GT.2d0*EXS4) THEN     ! Adjust (NH4)3H(SO4)2(s)
1440          CLC      = CLC - EXS4/2d0  ! more solid than excess
1441          GNH3     = GNH3 +1.5d0*EXS4! evaporate NH3 to gas phase
1442          GOTO 50
1443       ELSE                          ! less solid than excess
1444          GNH3     = GNH3 + 1.5d0*CLC! evaporate NH3 to gas phase
1445          EXS4     = EXS4 - 2d0*CLC  ! reduce excess
1446          CLC      = ZERO            ! zero salt concentration
1447       ENDIF
1448 !
1449       IF (CNH4HS4.GT.EXS4) THEN     ! Adjust NH4HSO4(s)
1450          CNH4HS4  = CNH4HS4 - EXS4  ! more solid than excess
1451          GNH3     = GNH3 + EXS4     ! evaporate NH3 to gas phase
1452          GOTO 50
1453       ELSE                          ! less solid than excess
1454          GNH3     = GNH3 + CNH4HS4  ! evaporate NH3 to gas phase
1455          EXS4     = EXS4  - CNH4HS4 ! reduce excess
1456          CNH4HS4  = ZERO            ! zero salt concentration
1457       ENDIF
1458 !
1459       IF (CNH42S4.GT.EXS4) THEN     ! Adjust (NH4)2SO4(s)
1460          CNH42S4  = CNH42S4- EXS4   ! more solid than excess
1461          GNH3     = GNH3 + 2.d0*EXS4! evaporate NH3 to gas phase
1462          GOTO 50
1463       ELSE                          ! less solid than excess
1464          GNH3     = GNH3+2.d0*CNH42S4 ! evaporate NH3 to gas phase
1465          EXS4     = EXS4  - CNH42S4 ! reduce excess
1466          CNH42S4  = ZERO            ! zero salt concentration
1467       ENDIF
1468 !
1469 ! *** RETURN **********************************************************
1470 !
1471  50   RETURN
1472       END
1473
1474 !=======================================================================
1475 !
1476 ! *** ISORROPIA CODE
1477 ! *** FUNCTION GETASR
1478 ! *** CALCULATES THE LIMITING NH4+/SO4 RATIO OF A SULFATE POOR SYSTEM
1479 !     (i.e. SULFATE RATIO = 2.0) FOR GIVEN SO4 LEVEL AND RH
1480 !
1481 ! *** COPYRIGHT 1996-2006, UNIVERSITY OF MIAMI, CARNEGIE MELLON UNIVERSITY,
1482 ! *** GEORGIA INSTITUTE OF TECHNOLOGY
1483 ! *** WRITTEN BY ATHANASIOS NENES
1484 ! *** UPDATED BY CHRISTOS FOUNTOUKIS
1485 !
1486 !=======================================================================
1487 !
1488       DOUBLE PRECISION FUNCTION GETASR (SO4I, RHI)
1489       USE ASRC
1490       PARAMETER (NSO4S=14, NRHS=20, NASRD=NSO4S*NRHS)
1491 !      COMMON /ASRC/ ASRAT(NASRD), ASSO4(NSO4S)
1492       DOUBLE PRECISION SO4I, RHI
1493 !CC
1494 !CC *** SOLVE USING FULL COMPUTATIONS, NOT LOOK-UP TABLES **************
1495 !CC
1496 !CC         W(2) = WAER(2)
1497 !CC         W(3) = WAER(2)*2.0001D0
1498 !CC         CALL CALCA2
1499 !CC         SULRATW = MOLAL(3)/WAER(2)
1500 !CC         CALL INIT1 (WI, RHI, TEMPI)   ! Re-initialize COMMON BLOCK
1501 !
1502 ! *** CALCULATE INDICES ************************************************
1503 !
1504       RAT    = SO4I/1.E-9
1505       A1     = INT(ALOG10(RAT))                   ! Magnitude of RAT
1506       IA1    = INT(RAT/2.5/10.0**A1)
1507 !
1508       INDS   = 4.0*A1 + MIN(IA1,4)
1509       INDS   = MIN(MAX(0, INDS), NSO4S-1) + 1     ! SO4 component of IPOS
1510 !
1511       INDR   = INT(99.0-RHI*100.0) + 1
1512       INDR   = MIN(MAX(1, INDR), NRHS)            ! RH component of IPOS
1513 !
1514 ! *** GET VALUE AND RETURN *********************************************
1515 !
1516       INDSL  = INDS
1517       INDSH  = MIN(INDSL+1, NSO4S)
1518       IPOSL  = (INDSL-1)*NRHS + INDR              ! Low position in array
1519       IPOSH  = (INDSH-1)*NRHS + INDR              ! High position in array
1520 !
1521       WF     = (SO4I-ASSO4(INDSL))/(ASSO4(INDSH)-ASSO4(INDSL) + 1e-7)
1522       WF     = MIN(MAX(WF, 0.0), 1.0)
1523 !
1524       GETASR = WF*ASRAT(IPOSH) + (1.0-WF)*ASRAT(IPOSL)
1525 !
1526 ! *** END OF FUNCTION GETASR *******************************************
1527 !
1528       RETURN
1529       END
1530
1531
1532
1533 !=======================================================================
1534 !
1535 ! *** ISORROPIA CODE
1536 ! *** BLOCK DATA AERSR
1537 ! *** CONTAINS DATA FOR AEROSOL SULFATE RATIO ARRAY NEEDED IN FUNCTION
1538 !     GETASR
1539 !
1540 ! *** COPYRIGHT 1996-2006, UNIVERSITY OF MIAMI, CARNEGIE MELLON UNIVERSITY,
1541 ! *** GEORGIA INSTITUTE OF TECHNOLOGY
1542 ! *** WRITTEN BY ATHANASIOS NENES
1543 ! *** UPDATED BY CHRISTOS FOUNTOUKIS
1544 !
1545 !=======================================================================
1546 !
1547       BLOCK DATA AERSR
1548       USE ASRC
1549       PARAMETER (NSO4S=14, NRHS=20, NASRD=NSO4S*NRHS)
1550 !      COMMON /ASRC/ ASRAT(NASRD), ASSO4(NSO4S)
1551 !
1552 !      DATA ASSO4/1.0E-9, 2.5E-9, 5.0E-9, 7.5E-9, 1.0E-8,   &
1553 !                 2.5E-8, 5.0E-8, 7.5E-8, 1.0E-7, 2.5E-7,    &
1554 !                 5.0E-7, 7.5E-7, 1.0E-6, 5.0E-6/
1555 !
1556 !      DATA (ASRAT(I), I=1,100)/   &
1557 !       1.020464, 0.9998130, 0.9960167, 0.9984423, 1.004004,   &
1558 !       1.010885,  1.018356,  1.026726,  1.034268, 1.043846,   &
1559 !       1.052933,  1.062230,  1.062213,  1.080050, 1.088350,   &
1560 !       1.096603,  1.104289,  1.111745,  1.094662, 1.121594,   &
1561 !       1.268909,  1.242444,  1.233815,  1.232088, 1.234020,   &
1562 !       1.238068,  1.243455,  1.250636,  1.258734, 1.267543,   &
1563 !       1.276948,  1.286642,  1.293337,  1.305592, 1.314726,   &
1564 !       1.323463,  1.333258,  1.343604,  1.344793, 1.355571,   &
1565 !       1.431463,  1.405204,  1.395791,  1.393190, 1.394403,   &
1566 !       1.398107,  1.403811,  1.411744,  1.420560, 1.429990,   &
1567 !       1.439742,  1.449507,  1.458986,  1.468403, 1.477394,   &
1568 !       1.487373,  1.495385,  1.503854,  1.512281, 1.520394,   &
1569 !       1.514464,  1.489699,  1.480686,  1.478187, 1.479446,   &
1570 !       1.483310,  1.489316,  1.497517,  1.506501, 1.515816,   &
1571 !       1.524724,  1.533950,  1.542758,  1.551730, 1.559587,   &
1572 !       1.568343,  1.575610,  1.583140,  1.590440, 1.596481,   &
1573 !       1.567743,  1.544426,  1.535928,  1.533645, 1.535016,   &
1574 !       1.539003,  1.545124,  1.553283,  1.561886, 1.570530,   &
1575 !       1.579234,  1.587813,  1.595956,  1.603901, 1.611349,   &
1576 !       1.618833,  1.625819,  1.632543,  1.639032, 1.645276/
1577
1578 !      DATA (ASRAT(I), I=101,200)/   &
1579 !       1.707390,  1.689553,  1.683198,  1.681810, 1.683490,   &
1580 !       1.687477,  1.693148,  1.700084,  1.706917, 1.713507,   &
1581 !       1.719952,  1.726190,  1.731985,  1.737544, 1.742673,   &
1582 !       1.747756,  1.752431,  1.756890,  1.761141, 1.765190,   &
1583 !       1.785657,  1.771851,  1.767063,  1.766229, 1.767901,   &
1584 !       1.771455,  1.776223,  1.781769,  1.787065, 1.792081,   &
1585 !       1.796922,  1.801561,  1.805832,  1.809896, 1.813622,   &
1586 !       1.817292,  1.820651,  1.823841,  1.826871, 1.829745,   &
1587 !       1.822215,  1.810497,  1.806496,  1.805898, 1.807480,   &
1588 !       1.810684,  1.814860,  1.819613,  1.824093, 1.828306,   &
1589 !       1.832352,  1.836209,  1.839748,  1.843105, 1.846175,   &
1590 !       1.849192,  1.851948,  1.854574,  1.857038, 1.859387,   &
1591 !       1.844588,  1.834208,  1.830701,  1.830233, 1.831727,   &
1592 !       1.834665,  1.838429,  1.842658,  1.846615, 1.850321,   &
1593 !       1.853869,  1.857243,  1.860332,  1.863257, 1.865928,   &
1594 !       1.868550,  1.870942,  1.873208,  1.875355, 1.877389,   &
1595 !       1.899556,  1.892637,  1.890367,  1.890165, 1.891317,   &
1596 !       1.893436,  1.896036,  1.898872,  1.901485, 1.903908,   &
1597 !       1.906212,  1.908391,  1.910375,  1.912248, 1.913952,   &
1598 !       1.915621,  1.917140,  1.918576,  1.919934, 1.921220/
1599
1600 !      DATA (ASRAT(I), I=201,280)/   &
1601 !       1.928264,  1.923245,  1.921625,  1.921523, 1.922421,   &
1602 !       1.924016,  1.925931,  1.927991,  1.929875, 1.931614,   &
1603 !       1.933262,  1.934816,  1.936229,  1.937560, 1.938769,   &
1604 !       1.939951,  1.941026,  1.942042,  1.943003, 1.943911,   &
1605 !       1.941205,  1.937060,  1.935734,  1.935666, 1.936430,   &
1606 !       1.937769,  1.939359,  1.941061,  1.942612, 1.944041,   &
1607 !       1.945393,  1.946666,  1.947823,  1.948911, 1.949900,   &
1608 !       1.950866,  1.951744,  1.952574,  1.953358, 1.954099,   &
1609 !       1.948985,  1.945372,  1.944221,  1.944171, 1.944850,   &
1610 !       1.946027,  1.947419,  1.948902,  1.950251, 1.951494,   &
1611 !       1.952668,  1.953773,  1.954776,  1.955719, 1.956576,   &
1612 !       1.957413,  1.958174,  1.958892,  1.959571, 1.960213,   &
1613 !       1.977193,  1.975540,  1.975023,  1.975015, 1.975346,   &
1614 !       1.975903,  1.976547,  1.977225,  1.977838, 1.978401,   &
1615 !       1.978930,  1.979428,  1.979879,  1.980302, 1.980686,   &
1616 !       1.981060,  1.981401,  1.981722,  1.982025, 1.982312/
1617 !!
1618 ! *** END OF BLOCK DATA AERSR ******************************************
1619 !
1620        END
1621
1622
1623 !=======================================================================
1624 !
1625 ! *** ISORROPIA CODE
1626 ! *** SUBROUTINE CALCHA
1627 ! *** CALCULATES CHLORIDES SPECIATION
1628 !
1629 !     HYDROCHLORIC ACID IN THE LIQUID PHASE IS ASSUMED A MINOR SPECIES,
1630 !     AND DOES NOT SIGNIFICANTLY PERTURB THE HSO4-SO4 EQUILIBRIUM. THE
1631 !     HYDROCHLORIC ACID DISSOLVED IS CALCULATED FROM THE
1632 !     HCL(G) <-> (H+) + (CL-)
1633 !     EQUILIBRIUM, USING THE (H+) FROM THE SULFATES.
1634 !
1635 ! *** COPYRIGHT 1996-2006, UNIVERSITY OF MIAMI, CARNEGIE MELLON UNIVERSITY,
1636 ! *** GEORGIA INSTITUTE OF TECHNOLOGY
1637 ! *** WRITTEN BY ATHANASIOS NENES
1638 ! *** UPDATED BY CHRISTOS FOUNTOUKIS
1639 !
1640 !=======================================================================
1641 !
1642       SUBROUTINE CALCHA
1643       USE ISRPIA
1644       IMPLICIT DOUBLE PRECISION (A-H, O-Z)
1645 !      INCLUDE 'ISRPIA.INC'
1646       DOUBLE PRECISION KAPA
1647 !C      CHARACTER ERRINF*40
1648 !
1649 ! *** CALCULATE HCL DISSOLUTION *****************************************
1650 !
1651       X    = W(5)
1652       DELT = 0.0d0
1653       IF (WATER.GT.TINY) THEN
1654          KAPA = MOLAL(1)
1655          ALFA = XK3*R*TEMP*(WATER/GAMA(11))**2.0
1656          DIAK = SQRT( (KAPA+ALFA)**2.0 + 4.0*ALFA*X)
1657          DELT = 0.5*(-(KAPA+ALFA) + DIAK)
1658 !C         IF (DELT/KAPA.GT.0.1d0) THEN
1659 !C            WRITE (ERRINF,'(1PE10.3)') DELT/KAPA*100.0
1660 !C            CALL PUSHERR (0033, ERRINF)
1661 !C         ENDIF
1662       ENDIF
1663 !
1664 ! *** CALCULATE HCL SPECIATION IN THE GAS PHASE *************************
1665 !
1666       GHCL     = MAX(X-DELT, 0.0d0)  ! GAS HCL
1667 !
1668 ! *** CALCULATE HCL SPECIATION IN THE LIQUID PHASE **********************
1669 !
1670       MOLAL(4) = DELT                ! CL-
1671       MOLAL(1) = MOLAL(1) + DELT     ! H+
1672 !
1673       RETURN
1674 !
1675 ! *** END OF SUBROUTINE CALCHA ******************************************
1676 !
1677       END
1678
1679
1680
1681
1682
1683 !=======================================================================
1684 !
1685 ! *** ISORROPIA CODE
1686 ! *** SUBROUTINE CALCHAP
1687 ! *** CALCULATES CHLORIDES SPECIATION
1688 !
1689 !     HYDROCHLORIC ACID IN THE LIQUID PHASE IS ASSUMED A MINOR SPECIES,
1690 !     THAT DOES NOT SIGNIFICANTLY PERTURB THE HSO4-SO4 EQUILIBRIUM.
1691 !     THE HYDROCHLORIC ACID DISSOLVED IS CALCULATED FROM THE
1692 !     HCL(G) -> HCL(AQ)   AND  HCL(AQ) ->  (H+) + (CL-)
1693 !     EQUILIBRIA, USING (H+) FROM THE SULFATES.
1694 !
1695 !     THIS IS THE VERSION USED BY THE INVERSE PROBLEM SOVER
1696 !
1697 ! *** COPYRIGHT 1996-2006, UNIVERSITY OF MIAMI, CARNEGIE MELLON UNIVERSITY,
1698 ! *** GEORGIA INSTITUTE OF TECHNOLOGY
1699 ! *** WRITTEN BY ATHANASIOS NENES
1700 ! *** UPDATED BY CHRISTOS FOUNTOUKIS
1701 !
1702 !=======================================================================
1703 !
1704       SUBROUTINE CALCHAP
1705       USE ISRPIA
1706       IMPLICIT DOUBLE PRECISION (A-H, O-Z)
1707 !      INCLUDE 'ISRPIA.INC'
1708 !
1709 ! *** IS THERE A LIQUID PHASE? ******************************************
1710 !
1711       IF (WATER.LE.TINY) RETURN
1712 !
1713 ! *** CALCULATE HCL SPECIATION IN THE GAS PHASE *************************
1714 !
1715       CALL CALCCLAQ (MOLAL(4), MOLAL(1), DELT)
1716       ALFA     = XK3*R*TEMP*(WATER/GAMA(11))**2.0
1717       GASAQ(3) = DELT
1718       MOLAL(1) = MOLAL(1) - DELT
1719       MOLAL(4) = MOLAL(4) - DELT
1720       GHCL     = MOLAL(1)*MOLAL(4)/ALFA
1721 !
1722       RETURN
1723 !
1724 ! *** END OF SUBROUTINE CALCHAP *****************************************
1725 !
1726       END
1727
1728
1729 !=======================================================================
1730 !
1731 ! *** ISORROPIA CODE
1732 ! *** SUBROUTINE CALCNA
1733 ! *** CALCULATES NITRATES SPECIATION
1734 !
1735 !     NITRIC ACID IN THE LIQUID PHASE IS ASSUMED A MINOR SPECIES, THAT
1736 !     DOES NOT SIGNIFICANTLY PERTURB THE HSO4-SO4 EQUILIBRIUM. THE NITRIC
1737 !     ACID DISSOLVED IS CALCULATED FROM THE HNO3(G) -> (H+) + (NO3-)
1738 !     EQUILIBRIUM, USING THE (H+) FROM THE SULFATES.
1739 !
1740 ! *** COPYRIGHT 1996-2006, UNIVERSITY OF MIAMI, CARNEGIE MELLON UNIVERSITY,
1741 ! *** GEORGIA INSTITUTE OF TECHNOLOGY
1742 ! *** WRITTEN BY ATHANASIOS NENES
1743 ! *** UPDATED BY CHRISTOS FOUNTOUKIS
1744 !
1745 !=======================================================================
1746 !
1747       SUBROUTINE CALCNA
1748       USE ISRPIA
1749       IMPLICIT DOUBLE PRECISION (A-H, O-Z)
1750 !      INCLUDE 'ISRPIA.INC'
1751       DOUBLE PRECISION KAPA
1752 !C      CHARACTER ERRINF*40
1753 !
1754 ! *** CALCULATE HNO3 DISSOLUTION ****************************************
1755 !
1756       X    = W(4)
1757       DELT = 0.0d0
1758       IF (WATER.GT.TINY) THEN
1759          KAPA = MOLAL(1)
1760          ALFA = XK4*R*TEMP*(WATER/GAMA(10))**2.0
1761          DIAK = SQRT( (KAPA+ALFA)**2.0 + 4.0*ALFA*X)
1762          DELT = 0.5*(-(KAPA+ALFA) + DIAK)
1763 !C         IF (DELT/KAPA.GT.0.1d0) THEN
1764 !C            WRITE (ERRINF,'(1PE10.3)') DELT/KAPA*100.0
1765 !C            CALL PUSHERR (0019, ERRINF)    ! WARNING ERROR: NO SOLUTION
1766 !C         ENDIF
1767       ENDIF
1768 !
1769 ! *** CALCULATE HNO3 SPECIATION IN THE GAS PHASE ************************
1770 !
1771       GHNO3    = MAX(X-DELT, 0.0d0)  ! GAS HNO3
1772 !
1773 ! *** CALCULATE HNO3 SPECIATION IN THE LIQUID PHASE *********************
1774 !
1775       MOLAL(7) = DELT                ! NO3-
1776       MOLAL(1) = MOLAL(1) + DELT     ! H+
1777 !
1778       RETURN
1779 !
1780 ! *** END OF SUBROUTINE CALCNA ******************************************
1781 !
1782       END
1783
1784
1785
1786 !=======================================================================
1787 !
1788 ! *** ISORROPIA CODE
1789 ! *** SUBROUTINE CALCNAP
1790 ! *** CALCULATES NITRATES SPECIATION
1791 !
1792 !     NITRIC ACID IN THE LIQUID PHASE IS ASSUMED A MINOR SPECIES, THAT
1793 !     DOES NOT SIGNIFICANTLY PERTURB THE HSO4-SO4 EQUILIBRIUM. THE NITRIC
1794 !     ACID DISSOLVED IS CALCULATED FROM THE HNO3(G) -> HNO3(AQ) AND
1795 !     HNO3(AQ) -> (H+) + (CL-) EQUILIBRIA, USING (H+) FROM THE SULFATES.
1796 !
1797 !     THIS IS THE VERSION USED BY THE INVERSE PROBLEM SOVER
1798 !
1799 ! *** COPYRIGHT 1996-2006, UNIVERSITY OF MIAMI, CARNEGIE MELLON UNIVERSITY,
1800 ! *** GEORGIA INSTITUTE OF TECHNOLOGY
1801 ! *** WRITTEN BY ATHANASIOS NENES
1802 ! *** UPDATED BY CHRISTOS FOUNTOUKIS
1803 !
1804 !=======================================================================
1805 !
1806       SUBROUTINE CALCNAP
1807       USE ISRPIA
1808       IMPLICIT DOUBLE PRECISION (A-H, O-Z)
1809 !      INCLUDE 'ISRPIA.INC'
1810 !
1811 ! *** IS THERE A LIQUID PHASE? ******************************************
1812 !
1813       IF (WATER.LE.TINY) RETURN
1814 !
1815 ! *** CALCULATE HNO3 SPECIATION IN THE GAS PHASE ************************
1816 !
1817       CALL CALCNIAQ (MOLAL(7), MOLAL(1), DELT)
1818       ALFA     = XK4*R*TEMP*(WATER/GAMA(10))**2.0
1819       GASAQ(3) = DELT
1820       MOLAL(1) = MOLAL(1) - DELT
1821       MOLAL(7) = MOLAL(7) - DELT
1822       GHNO3    = MOLAL(1)*MOLAL(7)/ALFA
1823
1824 !      write (*,*) ALFA, MOLAL(1), MOLAL(7), GHNO3, DELT
1825 !
1826       RETURN
1827 !
1828 ! *** END OF SUBROUTINE CALCNAP *****************************************
1829 !
1830       END
1831
1832 !=======================================================================
1833 !
1834 ! *** ISORROPIA CODE
1835 ! *** SUBROUTINE CALCNH3
1836 ! *** CALCULATES AMMONIA IN GAS PHASE
1837 !
1838 !     AMMONIA IN THE GAS PHASE IS ASSUMED A MINOR SPECIES, THAT
1839 !     DOES NOT SIGNIFICANTLY PERTURB THE AEROSOL EQUILIBRIUM.
1840 !     AMMONIA GAS IS CALCULATED FROM THE NH3(g) + (H+)(l) <==> (NH4+)(l)
1841 !     EQUILIBRIUM, USING (H+), (NH4+) FROM THE AEROSOL SOLUTION.
1842 !
1843 !     THIS IS THE VERSION USED BY THE DIRECT PROBLEM
1844 !
1845 ! *** COPYRIGHT 1996-2006, UNIVERSITY OF MIAMI, CARNEGIE MELLON UNIVERSITY,
1846 ! *** GEORGIA INSTITUTE OF TECHNOLOGY
1847 ! *** WRITTEN BY ATHANASIOS NENES
1848 ! *** UPDATED BY CHRISTOS FOUNTOUKIS
1849 !
1850 !=======================================================================
1851 !
1852       SUBROUTINE CALCNH3
1853       USE ISRPIA
1854       IMPLICIT DOUBLE PRECISION (A-H, O-Z)
1855 !      INCLUDE 'ISRPIA.INC'
1856 !
1857 ! *** IS THERE A LIQUID PHASE? ******************************************
1858 !
1859       IF (WATER.LE.TINY) RETURN
1860 !
1861 ! *** CALCULATE NH3 SUBLIMATION *****************************************
1862 !
1863       A1   = (XK2/XKW)*R*TEMP*(GAMA(10)/GAMA(5))**2.0
1864       CHI1 = MOLAL(3)
1865       CHI2 = MOLAL(1)
1866 !
1867       BB   =(CHI2 + ONE/A1)          ! a=1; b!=1; c!=1
1868       CC   =-CHI1/A1
1869       DIAK = SQRT(BB*BB - 4.D0*CC)   ! Always > 0
1870       PSI  = 0.5*(-BB + DIAK)        ! One positive root
1871       PSI  = MAX(TINY, MIN(PSI,CHI1))! Constrict in acceptible range
1872 !
1873 ! *** CALCULATE NH3 SPECIATION IN THE GAS PHASE *************************
1874 !
1875       GNH3     = PSI                 ! GAS HNO3
1876 !
1877 ! *** CALCULATE NH3 AFFECT IN THE LIQUID PHASE **************************
1878 !
1879       MOLAL(3) = CHI1 - PSI          ! NH4+
1880       MOLAL(1) = CHI2 + PSI          ! H+
1881 !
1882       RETURN
1883 !
1884 ! *** END OF SUBROUTINE CALCNH3 *****************************************
1885 !
1886       END
1887
1888
1889
1890 !=======================================================================
1891 !
1892 ! *** ISORROPIA CODE
1893 ! *** SUBROUTINE CALCNH3P
1894 ! *** CALCULATES AMMONIA IN GAS PHASE
1895 !
1896 !     AMMONIA GAS IS CALCULATED FROM THE NH3(g) + (H+)(l) <==> (NH4+)(l)
1897 !     EQUILIBRIUM, USING (H+), (NH4+) FROM THE AEROSOL SOLUTION.
1898 !
1899 !     THIS IS THE VERSION USED BY THE INVERSE PROBLEM SOLVER
1900 !
1901 ! *** COPYRIGHT 1996-2006, UNIVERSITY OF MIAMI, CARNEGIE MELLON UNIVERSITY,
1902 ! *** GEORGIA INSTITUTE OF TECHNOLOGY
1903 ! *** WRITTEN BY ATHANASIOS NENES
1904 ! *** UPDATED BY CHRISTOS FOUNTOUKIS
1905 !
1906 !=======================================================================
1907 !
1908       SUBROUTINE CALCNH3P
1909       USE ISRPIA
1910       IMPLICIT DOUBLE PRECISION (A-H, O-Z)
1911 !      INCLUDE 'ISRPIA.INC'
1912 !
1913 ! *** IS THERE A LIQUID PHASE? ******************************************
1914 !
1915       IF (WATER.LE.TINY) RETURN
1916 !
1917 ! *** CALCULATE NH3 GAS PHASE CONCENTRATION *****************************
1918 !
1919       A1   = (XK2/XKW)*R*TEMP*(GAMA(10)/GAMA(5))**2.0
1920       GNH3 = MOLAL(3)/MOLAL(1)/A1
1921 !
1922       RETURN
1923 !
1924 ! *** END OF SUBROUTINE CALCNH3P ****************************************
1925 !
1926       END
1927
1928
1929 !=======================================================================
1930 !
1931 ! *** ISORROPIA CODE
1932 ! *** SUBROUTINE CALCNHA
1933 !
1934 !     THIS SUBROUTINE CALCULATES THE DISSOLUTION OF HCL, HNO3 AT
1935 !     THE PRESENCE OF (H,SO4). HCL, HNO3 ARE CONSIDERED MINOR SPECIES,
1936 !     THAT DO NOT SIGNIFICANTLY AFFECT THE EQUILIBRIUM POINT.
1937 !
1938 ! *** COPYRIGHT 1996-2006, UNIVERSITY OF MIAMI, CARNEGIE MELLON UNIVERSITY,
1939 ! *** GEORGIA INSTITUTE OF TECHNOLOGY
1940 ! *** WRITTEN BY ATHANASIOS NENES
1941 ! *** UPDATED BY CHRISTOS FOUNTOUKIS
1942 !
1943 !=======================================================================
1944 !
1945       SUBROUTINE CALCNHA
1946       USE ISRPIA
1947       IMPLICIT DOUBLE PRECISION (A-H, O-Z)
1948 !      INCLUDE 'ISRPIA.INC'
1949       DOUBLE PRECISION M1, M2, M3
1950       CHARACTER ERRINF*40
1951 !
1952 ! *** SPECIAL CASE; WATER=ZERO ******************************************
1953 !
1954       IF (WATER.LE.TINY) THEN
1955          GOTO 55
1956 !
1957 ! *** SPECIAL CASE; HCL=HNO3=ZERO ***************************************
1958 !
1959       ELSEIF (W(5).LE.TINY .AND. W(4).LE.TINY) THEN
1960          GOTO 60
1961 !
1962 ! *** SPECIAL CASE; HCL=ZERO ********************************************
1963 !
1964       ELSE IF (W(5).LE.TINY) THEN
1965          CALL CALCNA              ! CALL HNO3 DISSOLUTION ROUTINE
1966          GOTO 60
1967 !
1968 ! *** SPECIAL CASE; HNO3=ZERO *******************************************
1969 !
1970       ELSE IF (W(4).LE.TINY) THEN
1971          CALL CALCHA              ! CALL HCL DISSOLUTION ROUTINE
1972          GOTO 60
1973       ENDIF
1974 !
1975 ! *** CALCULATE EQUILIBRIUM CONSTANTS ***********************************
1976 !
1977       A3 = XK4*R*TEMP*(WATER/GAMA(10))**2.0   ! HNO3
1978       A4 = XK3*R*TEMP*(WATER/GAMA(11))**2.0   ! HCL
1979 !
1980 ! *** CALCULATE CUBIC EQUATION COEFFICIENTS *****************************
1981 !
1982       DELCL = ZERO
1983       DELNO = ZERO
1984 !
1985       OMEGA = MOLAL(1)       ! H+
1986       CHI3  = W(4)           ! HNO3
1987       CHI4  = W(5)           ! HCL
1988 !
1989       C1    = A3*CHI3
1990       C2    = A4*CHI4
1991       C3    = A3 - A4
1992 !
1993       M1    = (C1 + C2 + (OMEGA+A4)*C3)/C3
1994       M2    = ((OMEGA+A4)*C2 - A4*C3*CHI4)/C3
1995       M3    =-A4*C2*CHI4/C3
1996 !
1997 ! *** CALCULATE ROOTS ***************************************************
1998 !
1999       CALL POLY3 (M1, M2, M3, DELCL, ISLV) ! HCL DISSOLUTION
2000       IF (ISLV.NE.0) THEN
2001          DELCL = TINY       ! TINY AMOUNTS OF HCL ASSUMED WHEN NO ROOT
2002          WRITE (ERRINF,'(1PE7.1)') TINY
2003          CALL PUSHERR (0022, ERRINF)    ! WARNING ERROR: NO SOLUTION
2004       ENDIF
2005       DELCL = MIN(DELCL, CHI4)
2006 !
2007       DELNO = C1*DELCL/(C2 + C3*DELCL)
2008       DELNO = MIN(DELNO, CHI3)
2009 !
2010       IF (DELCL.LT.ZERO .OR. DELNO.LT.ZERO .OR.   &
2011          DELCL.GT.CHI4 .OR. DELNO.GT.CHI3       ) THEN
2012          DELCL = TINY  ! TINY AMOUNTS OF HCL ASSUMED WHEN NO ROOT
2013          DELNO = TINY
2014          WRITE (ERRINF,'(1PE7.1)') TINY
2015          CALL PUSHERR (0022, ERRINF)    ! WARNING ERROR: NO SOLUTION
2016       ENDIF
2017 !CC
2018 !CC *** COMPARE DELTA TO TOTAL H+ ; ESTIMATE EFFECT TO HSO4 ***************
2019 !CC
2020 !C      IF ((DELCL+DELNO)/MOLAL(1).GT.0.1d0) THEN
2021 !C         WRITE (ERRINF,'(1PE10.3)') (DELCL+DELNO)/MOLAL(1)*100.0
2022 !C         CALL PUSHERR (0021, ERRINF)
2023 !C      ENDIF
2024 !
2025 ! *** EFFECT ON LIQUID PHASE ********************************************
2026 !
2027 50    MOLAL(1) = MOLAL(1) + (DELNO+DELCL)  ! H+   CHANGE
2028       MOLAL(4) = MOLAL(4) + DELCL          ! CL-  CHANGE
2029       MOLAL(7) = MOLAL(7) + DELNO          ! NO3- CHANGE
2030 !
2031 ! *** EFFECT ON GAS PHASE ***********************************************
2032 !
2033 55    GHCL     = MAX(W(5) - MOLAL(4), TINY)
2034       GHNO3    = MAX(W(4) - MOLAL(7), TINY)
2035 !
2036 60    RETURN
2037 !
2038 ! *** END OF SUBROUTINE CALCNHA *****************************************
2039 !
2040       END
2041
2042
2043
2044 !=======================================================================
2045 !
2046 ! *** ISORROPIA CODE
2047 ! *** SUBROUTINE CALCNHP
2048 !
2049 !     THIS SUBROUTINE CALCULATES THE GAS PHASE NITRIC AND HYDROCHLORIC
2050 !     ACID. CONCENTRATIONS ARE CALCULATED FROM THE DISSOLUTION
2051 !     EQUILIBRIA, USING (H+), (Cl-), (NO3-) IN THE AEROSOL PHASE.
2052 !
2053 !     THIS IS THE VERSION USED BY THE INVERSE PROBLEM SOLVER
2054 !
2055 ! *** COPYRIGHT 1996-2006, UNIVERSITY OF MIAMI, CARNEGIE MELLON UNIVERSITY,
2056 ! *** GEORGIA INSTITUTE OF TECHNOLOGY
2057 ! *** WRITTEN BY ATHANASIOS NENES
2058 ! *** UPDATED BY CHRISTOS FOUNTOUKIS
2059 !
2060 !=======================================================================
2061 !
2062       SUBROUTINE CALCNHP
2063       USE ISRPIA
2064       IMPLICIT DOUBLE PRECISION (A-H, O-Z)
2065 !      INCLUDE 'ISRPIA.INC'
2066 !
2067 ! *** IS THERE A LIQUID PHASE? ******************************************
2068 !
2069       IF (WATER.LE.TINY) RETURN
2070 !
2071 ! *** CALCULATE EQUILIBRIUM CONSTANTS ***********************************
2072 !
2073       A3       = XK3*R*TEMP*(WATER/GAMA(11))**2.0
2074       A4       = XK4*R*TEMP*(WATER/GAMA(10))**2.0
2075       MOLAL(1) = MOLAL(1) + WAER(4) + WAER(5)  ! H+ increases because NO3, Cl are added.
2076 !
2077 ! *** CALCULATE CONCENTRATIONS ******************************************
2078 ! *** ASSUME THAT 'DELT' FROM HNO3 >> 'DELT' FROM HCL
2079 !
2080       CALL CALCNIAQ (WAER(4), MOLAL(1)+MOLAL(7)+MOLAL(4), DELT)
2081       MOLAL(1) = MOLAL(1) - DELT
2082       MOLAL(7) = WAER(4)  - DELT  ! NO3- = Waer(4) minus any turned into (HNO3aq)
2083       GASAQ(3) = DELT
2084 !
2085       CALL CALCCLAQ (WAER(5), MOLAL(1)+MOLAL(7)+MOLAL(4), DELT)
2086       MOLAL(1) = MOLAL(1) - DELT
2087       MOLAL(4) = WAER(5)  - DELT  ! Cl- = Waer(4) minus any turned into (HNO3aq)
2088       GASAQ(2) = DELT
2089 !
2090       GHNO3    = MOLAL(1)*MOLAL(7)/A4
2091       GHCL     = MOLAL(1)*MOLAL(4)/A3
2092 !
2093       RETURN
2094 !
2095 ! *** END OF SUBROUTINE CALCNHP *****************************************
2096 !
2097       END
2098
2099 !=======================================================================
2100 !
2101 ! *** ISORROPIA CODE
2102 ! *** SUBROUTINE CALCAMAQ
2103 ! *** THIS SUBROUTINE CALCULATES THE NH3(aq) GENERATED FROM (H,NH4+).
2104 !
2105 ! *** COPYRIGHT 1996-2006, UNIVERSITY OF MIAMI, CARNEGIE MELLON UNIVERSITY,
2106 ! *** GEORGIA INSTITUTE OF TECHNOLOGY
2107 ! *** WRITTEN BY ATHANASIOS NENES
2108 ! *** UPDATED BY CHRISTOS FOUNTOUKIS
2109 !
2110 !=======================================================================
2111 !
2112       SUBROUTINE CALCAMAQ (NH4I, OHI, DELT)
2113       USE ISRPIA
2114       IMPLICIT DOUBLE PRECISION (A-H, O-Z)
2115 !      INCLUDE 'ISRPIA.INC'
2116       DOUBLE PRECISION NH4I
2117 !C      CHARACTER ERRINF*40
2118 !
2119 ! *** EQUILIBRIUM CONSTANTS
2120 !
2121       A22  = XK22/XKW/WATER*(GAMA(8)/GAMA(9))**2. ! GAMA(NH3) ASSUMED 1
2122       AKW  = XKW *RH*WATER*WATER
2123 !
2124 ! *** FIND ROOT
2125 !
2126       OM1  = NH4I
2127       OM2  = OHI
2128       BB   =-(OM1+OM2+A22*AKW)
2129       CC   = OM1*OM2
2130       DD   = SQRT(BB*BB-4.D0*CC)
2131
2132       DEL1 = 0.5D0*(-BB - DD)
2133       DEL2 = 0.5D0*(-BB + DD)
2134 !
2135 ! *** GET APPROPRIATE ROOT.
2136 !
2137       IF (DEL1.LT.ZERO) THEN
2138          IF (DEL2.GT.NH4I .OR. DEL2.GT.OHI) THEN
2139             DELT = ZERO
2140          ELSE
2141             DELT = DEL2
2142          ENDIF
2143       ELSE
2144          DELT = DEL1
2145       ENDIF
2146 !C
2147 !C *** COMPARE DELTA TO TOTAL NH4+ ; ESTIMATE EFFECT *********************
2148 !C
2149 !C      IF (DELTA/HYD.GT.0.1d0) THEN
2150 !C         WRITE (ERRINF,'(1PE10.3)') DELTA/HYD*100.0
2151 !C         CALL PUSHERR (0020, ERRINF)
2152 !C      ENDIF
2153 !
2154       RETURN
2155 !
2156 ! *** END OF SUBROUTINE CALCAMAQ ****************************************
2157 !
2158       END
2159
2160
2161
2162 !=======================================================================
2163 !
2164 ! *** ISORROPIA CODE
2165 ! *** SUBROUTINE CALCAMAQ2
2166 !
2167 !     THIS SUBROUTINE CALCULATES THE NH3(aq) GENERATED FROM (H,NH4+).
2168 !
2169 ! *** COPYRIGHT 1996-2006, UNIVERSITY OF MIAMI, CARNEGIE MELLON UNIVERSITY,
2170 ! *** GEORGIA INSTITUTE OF TECHNOLOGY
2171 ! *** WRITTEN BY ATHANASIOS NENES
2172 ! *** UPDATED BY CHRISTOS FOUNTOUKIS
2173 !
2174 !=======================================================================
2175 !
2176       SUBROUTINE CALCAMAQ2 (GGNH3, NH4I, OHI, NH3AQ)
2177       USE ISRPIA
2178       IMPLICIT DOUBLE PRECISION (A-H, O-Z)
2179 !      INCLUDE 'ISRPIA.INC'
2180       DOUBLE PRECISION NH4I, NH3AQ
2181 !
2182 ! *** EQUILIBRIUM CONSTANTS
2183 !
2184       A22  = XK22/XKW/WATER*(GAMA(8)/GAMA(9))**2. ! GAMA(NH3) ASSUMED 1
2185       AKW  = XKW *RH*WATER*WATER
2186 !
2187 ! *** FIND ROOT
2188 !
2189       ALF1 = NH4I - GGNH3
2190       ALF2 = GGNH3
2191       BB   = ALF1 + A22*AKW
2192       CC   =-A22*AKW*ALF2
2193       DEL  = 0.5D0*(-BB + SQRT(BB*BB-4.D0*CC))
2194 !
2195 ! *** ADJUST CONCENTRATIONS
2196 !
2197       NH4I  = ALF1 + DEL
2198       OHI   = DEL
2199       IF (OHI.LE.TINY) OHI = SQRT(AKW)   ! If solution is neutral.
2200       NH3AQ = ALF2 - DEL
2201 !
2202       RETURN
2203 !
2204 ! *** END OF SUBROUTINE CALCAMAQ2 ****************************************
2205 !
2206       END
2207
2208
2209
2210 !=======================================================================
2211 !
2212 ! *** ISORROPIA CODE
2213 ! *** SUBROUTINE CALCCLAQ
2214 !
2215 !     THIS SUBROUTINE CALCULATES THE HCL(aq) GENERATED FROM (H+,CL-).
2216 !
2217 ! *** COPYRIGHT 1996-2006, UNIVERSITY OF MIAMI, CARNEGIE MELLON UNIVERSITY,
2218 ! *** GEORGIA INSTITUTE OF TECHNOLOGY
2219 ! *** WRITTEN BY ATHANASIOS NENES
2220 ! *** UPDATED BY CHRISTOS FOUNTOUKIS
2221 !
2222 !=======================================================================
2223 !
2224       SUBROUTINE CALCCLAQ (CLI, HI, DELT)
2225       USE ISRPIA
2226       IMPLICIT DOUBLE PRECISION (A-H, O-Z)
2227 !      INCLUDE 'ISRPIA.INC'
2228       DOUBLE PRECISION CLI
2229 !
2230 ! *** EQUILIBRIUM CONSTANTS
2231 !
2232       A32  = XK32*WATER/(GAMA(11))**2. ! GAMA(HCL) ASSUMED 1
2233 !
2234 ! *** FIND ROOT
2235 !
2236       OM1  = CLI
2237       OM2  = HI
2238       BB   =-(OM1+OM2+A32)
2239       CC   = OM1*OM2
2240       DD   = SQRT(BB*BB-4.D0*CC)
2241
2242       DEL1 = 0.5D0*(-BB - DD)
2243       DEL2 = 0.5D0*(-BB + DD)
2244 !
2245 ! *** GET APPROPRIATE ROOT.
2246 !
2247       IF (DEL1.LT.ZERO) THEN
2248          IF (DEL2.LT.ZERO .OR. DEL2.GT.CLI .OR. DEL2.GT.HI) THEN
2249             DELT = ZERO
2250          ELSE
2251             DELT = DEL2
2252          ENDIF
2253       ELSE
2254          DELT = DEL1
2255       ENDIF
2256 !
2257       RETURN
2258 !
2259 ! *** END OF SUBROUTINE CALCCLAQ ****************************************
2260 !
2261       END
2262
2263
2264
2265 !=======================================================================
2266 !
2267 ! *** ISORROPIA CODE
2268 ! *** SUBROUTINE CALCCLAQ2
2269 !
2270 !     THIS SUBROUTINE CALCULATES THE HCL(aq) GENERATED FROM (H+,CL-).
2271 !
2272 ! *** COPYRIGHT 1996-2006, UNIVERSITY OF MIAMI, CARNEGIE MELLON UNIVERSITY,
2273 ! *** GEORGIA INSTITUTE OF TECHNOLOGY
2274 ! *** WRITTEN BY ATHANASIOS NENES
2275 ! *** UPDATED BY CHRISTOS FOUNTOUKIS
2276 !
2277 !=======================================================================
2278 !
2279       SUBROUTINE CALCCLAQ2 (GGCL, CLI, HI, CLAQ)
2280       USE ISRPIA
2281       IMPLICIT DOUBLE PRECISION (A-H, O-Z)
2282 !      INCLUDE 'ISRPIA.INC'
2283       DOUBLE PRECISION CLI
2284 !
2285 ! *** EQUILIBRIUM CONSTANTS
2286 !
2287       A32  = XK32*WATER/(GAMA(11))**2. ! GAMA(HCL) ASSUMED 1
2288       AKW  = XKW *RH*WATER*WATER
2289 !
2290 ! *** FIND ROOT
2291 !
2292       ALF1  = CLI - GGCL
2293       ALF2  = GGCL
2294       COEF  = (ALF1+A32)
2295       DEL1  = 0.5*(-COEF + SQRT(COEF*COEF+4.D0*A32*ALF2))
2296 !
2297 ! *** CORRECT CONCENTRATIONS
2298 !
2299       CLI  = ALF1 + DEL1
2300       HI   = DEL1
2301       IF (HI.LE.TINY) HI = SQRT(AKW)   ! If solution is neutral.
2302       CLAQ = ALF2 - DEL1
2303 !
2304       RETURN
2305 !
2306 ! *** END OF SUBROUTINE CALCCLAQ2 ****************************************
2307 !
2308       END
2309
2310
2311
2312 !=======================================================================
2313 !
2314 ! *** ISORROPIA CODE
2315 ! *** SUBROUTINE CALCNIAQ
2316 !
2317 !     THIS SUBROUTINE CALCULATES THE HNO3(aq) GENERATED FROM (H,NO3-).
2318 !
2319 ! *** COPYRIGHT 1996-2006, UNIVERSITY OF MIAMI, CARNEGIE MELLON UNIVERSITY,
2320 ! *** GEORGIA INSTITUTE OF TECHNOLOGY
2321 ! *** WRITTEN BY ATHANASIOS NENES
2322 ! *** UPDATED BY CHRISTOS FOUNTOUKIS
2323 !
2324 !=======================================================================
2325 !
2326       SUBROUTINE CALCNIAQ (NO3I, HI, DELT)
2327       USE ISRPIA
2328       IMPLICIT DOUBLE PRECISION (A-H, O-Z)
2329 !      INCLUDE 'ISRPIA.INC'
2330       DOUBLE PRECISION NO3I, HI, DELT
2331 !
2332 ! *** EQUILIBRIUM CONSTANTS
2333 !
2334       A42  = XK42*WATER/(GAMA(10))**2. ! GAMA(HNO3) ASSUMED 1
2335 !
2336 ! *** FIND ROOT
2337 !
2338       OM1  = NO3I
2339       OM2  = HI
2340       BB   =-(OM1+OM2+A42)
2341       CC   = OM1*OM2
2342       DD   = SQRT(BB*BB-4.D0*CC)
2343
2344       DEL1 = 0.5D0*(-BB - DD)
2345       DEL2 = 0.5D0*(-BB + DD)
2346 !
2347 ! *** GET APPROPRIATE ROOT.
2348 !
2349       IF (DEL1.LT.ZERO .OR. DEL1.GT.HI .OR. DEL1.GT.NO3I) THEN
2350          print *, DELT
2351          DELT = ZERO
2352       ELSE
2353          DELT = DEL1
2354          RETURN
2355       ENDIF
2356 !
2357       IF (DEL2.LT.ZERO .OR. DEL2.GT.NO3I .OR. DEL2.GT.HI) THEN
2358          DELT = ZERO
2359       ELSE
2360          DELT = DEL2
2361       ENDIF
2362 !
2363       RETURN
2364 !
2365 ! *** END OF SUBROUTINE CALCNIAQ ****************************************
2366 !
2367       END
2368
2369
2370
2371 !=======================================================================
2372 !
2373 ! *** ISORROPIA CODE
2374 ! *** SUBROUTINE CALCNIAQ2
2375 !
2376 !     THIS SUBROUTINE CALCULATES THE UNDISSOCIATED HNO3(aq)
2377 !
2378 ! *** COPYRIGHT 1996-2006, UNIVERSITY OF MIAMI, CARNEGIE MELLON UNIVERSITY,
2379 ! *** GEORGIA INSTITUTE OF TECHNOLOGY
2380 ! *** WRITTEN BY ATHANASIOS NENES
2381 ! *** UPDATED BY CHRISTOS FOUNTOUKIS
2382 !
2383 !=======================================================================
2384 !
2385       SUBROUTINE CALCNIAQ2 (GGNO3, NO3I, HI, NO3AQ)
2386       USE ISRPIA
2387       IMPLICIT DOUBLE PRECISION (A-H, O-Z)
2388 !      INCLUDE 'ISRPIA.INC'
2389       DOUBLE PRECISION NO3I, NO3AQ
2390 !
2391 ! *** EQUILIBRIUM CONSTANTS
2392 !
2393       A42  = XK42*WATER/(GAMA(10))**2. ! GAMA(HNO3) ASSUMED 1
2394       AKW  = XKW *RH*WATER*WATER
2395 !
2396 ! *** FIND ROOT
2397 !
2398       ALF1  = NO3I - GGNO3
2399       ALF2  = GGNO3
2400       ALF3  = HI
2401 !
2402       BB    = ALF3 + ALF1 + A42
2403       CC    = ALF3*ALF1 - A42*ALF2
2404       DEL1  = 0.5*(-BB + SQRT(BB*BB-4.D0*CC))
2405 !
2406 ! *** CORRECT CONCENTRATIONS
2407 !
2408       NO3I  = ALF1 + DEL1
2409       HI    = ALF3 + DEL1
2410       IF (HI.LE.TINY) HI = SQRT(AKW)   ! If solution is neutral.
2411       NO3AQ = ALF2 - DEL1
2412 !
2413       RETURN
2414 !
2415 ! *** END OF SUBROUTINE CALCNIAQ2 ****************************************
2416 !
2417       END
2418
2419
2420 !=======================================================================
2421 !
2422 ! *** ISORROPIA CODE
2423 ! *** SUBROUTINE CALCMR
2424 ! *** THIS SUBROUTINE CALCULATES:
2425 !     1. ION PAIR CONCENTRATIONS (FROM [MOLAR] ARRAY)
2426 !     2. WATER CONTENT OF LIQUID AEROSOL PHASE (FROM ZSR CORRELATION)
2427 !
2428 ! *** COPYRIGHT 1996-2006, UNIVERSITY OF MIAMI, CARNEGIE MELLON UNIVERSITY,
2429 ! *** GEORGIA INSTITUTE OF TECHNOLOGY
2430 ! *** WRITTEN BY ATHANASIOS NENES
2431 ! *** UPDATED BY CHRISTOS FOUNTOUKIS
2432 !
2433 !=======================================================================
2434 !
2435       SUBROUTINE CALCMR
2436       USE SOLUT
2437       USE ISRPIA
2438       IMPLICIT DOUBLE PRECISION (A-H, O-Z)
2439 !      INCLUDE 'ISRPIA.INC'
2440 !      COMMON /SOLUT/ CHI1, CHI2, CHI3, CHI4, CHI5, CHI6, CHI7, CHI8,   &
2441 !                     PSI1, PSI2, PSI3, PSI4, PSI5, PSI6, PSI7, PSI8,   &
2442 !                     A1,   A2,   A3,   A4,   A5,   A6,   A7,   A8
2443       CHARACTER SC*1
2444 !
2445 ! *** CALCULATE ION PAIR CONCENTRATIONS ACCORDING TO SPECIFIC CASE ****
2446 !
2447       SC =SCASE(1:1)                   ! SULRAT & SODRAT case
2448 !
2449 ! *** NH4-SO4 SYSTEM ; SULFATE POOR CASE
2450 !
2451       IF (SC.EQ.'A') THEN
2452          MOLALR(4) = MOLAL(5)+MOLAL(6) ! (NH4)2SO4 - CORRECT FOR SO4 TO HSO4
2453 !
2454 ! *** NH4-SO4 SYSTEM ; SULFATE RICH CASE ; NO FREE ACID
2455 !
2456       ELSE IF (SC.EQ.'B') THEN
2457          SO4I  = MOLAL(5)-MOLAL(1)     ! CORRECT FOR HSO4 DISSOCIATION
2458          HSO4I = MOLAL(6)+MOLAL(1)
2459          IF (SO4I.LT.HSO4I) THEN
2460             MOLALR(13) = SO4I                   ! [LC] = [SO4]
2461             MOLALR(9)  = MAX(HSO4I-SO4I, ZERO)  ! NH4HSO4
2462          ELSE
2463             MOLALR(13) = HSO4I                  ! [LC] = [HSO4]
2464             MOLALR(4)  = MAX(SO4I-HSO4I, ZERO)  ! (NH4)2SO4
2465          ENDIF
2466 !
2467 ! *** NH4-SO4 SYSTEM ; SULFATE RICH CASE ; FREE ACID
2468 !
2469       ELSE IF (SC.EQ.'C') THEN
2470          MOLALR(9) = MOLAL(3)                     ! NH4HSO4
2471          MOLALR(7) = MAX(W(2)-W(3), ZERO)         ! H2SO4
2472 !
2473 ! *** NH4-SO4-NO3 SYSTEM ; SULFATE POOR CASE
2474 !
2475       ELSE IF (SC.EQ.'D') THEN
2476          MOLALR(4) = MOLAL(5) + MOLAL(6)          ! (NH4)2SO4
2477          AML5      = MOLAL(3)-2.D0*MOLALR(4)      ! "free" NH4
2478          MOLALR(5) = MAX(MIN(AML5,MOLAL(7)), ZERO)! NH4NO3 = MIN("free", NO3)
2479 !
2480 ! *** NH4-SO4-NO3 SYSTEM ; SULFATE RICH CASE ; NO FREE ACID
2481 !
2482       ELSE IF (SC.EQ.'E') THEN
2483          SO4I  = MAX(MOLAL(5)-MOLAL(1),ZERO)      ! FROM HSO4 DISSOCIATION
2484          HSO4I = MOLAL(6)+MOLAL(1)
2485          IF (SO4I.LT.HSO4I) THEN
2486             MOLALR(13) = SO4I                     ! [LC] = [SO4]
2487             MOLALR(9)  = MAX(HSO4I-SO4I, ZERO)    ! NH4HSO4
2488          ELSE
2489             MOLALR(13) = HSO4I                    ! [LC] = [HSO4]
2490             MOLALR(4)  = MAX(SO4I-HSO4I, ZERO)    ! (NH4)2SO4
2491          ENDIF
2492 !
2493 ! *** NH4-SO4-NO3 SYSTEM ; SULFATE RICH CASE ; FREE ACID
2494 !
2495       ELSE IF (SC.EQ.'F') THEN
2496          MOLALR(9) = MOLAL(3)                              ! NH4HSO4
2497          MOLALR(7) = MAX(MOLAL(5)+MOLAL(6)-MOLAL(3),ZERO)  ! H2SO4
2498 !
2499 ! *** NA-NH4-SO4-NO3-CL SYSTEM ; SULFATE POOR ; SODIUM POOR CASE
2500 !
2501       ELSE IF (SC.EQ.'G') THEN
2502          MOLALR(2) = 0.5*MOLAL(2)                          ! NA2SO4
2503          TOTS4     = MOLAL(5)+MOLAL(6)                     ! Total SO4
2504          MOLALR(4) = MAX(TOTS4 - MOLALR(2), ZERO)          ! (NH4)2SO4
2505          FRNH4     = MAX(MOLAL(3) - 2.D0*MOLALR(4), ZERO)
2506          MOLALR(5) = MIN(MOLAL(7),FRNH4)                   ! NH4NO3
2507          FRNH4     = MAX(FRNH4 - MOLALR(5), ZERO)
2508          MOLALR(6) = MIN(MOLAL(4), FRNH4)                  ! NH4CL
2509 !
2510 ! *** NA-NH4-SO4-NO3-CL SYSTEM ; SULFATE POOR ; SODIUM RICH CASE
2511 ! *** RETREIVE DISSOLVED SALTS DIRECTLY FROM COMMON BLOCK /SOLUT/
2512 !
2513       ELSE IF (SC.EQ.'H') THEN
2514          MOLALR(1) = PSI7                                  ! NACL
2515          MOLALR(2) = PSI1                                  ! NA2SO4
2516          MOLALR(3) = PSI8                                  ! NANO3
2517          MOLALR(4) = ZERO                                  ! (NH4)2SO4
2518          FRNO3     = MAX(MOLAL(7) - MOLALR(3), ZERO)       ! "FREE" NO3
2519          FRCL      = MAX(MOLAL(4) - MOLALR(1), ZERO)       ! "FREE" CL
2520          MOLALR(5) = MIN(MOLAL(3),FRNO3)                   ! NH4NO3
2521          FRNH4     = MAX(MOLAL(3) - MOLALR(5), ZERO)       ! "FREE" NH3
2522          MOLALR(6) = MIN(FRCL, FRNH4)                      ! NH4CL
2523 !
2524 ! *** NA-NH4-SO4-NO3-CL SYSTEM ; SULFATE RICH CASE ; NO FREE ACID
2525 ! *** RETREIVE DISSOLVED SALTS DIRECTLY FROM COMMON BLOCK /SOLUT/
2526 !
2527       ELSE IF (SC.EQ.'I') THEN
2528          MOLALR(04) = PSI5                                 ! (NH4)2SO4
2529          MOLALR(02) = PSI4                                 ! NA2SO4
2530          MOLALR(09) = PSI1                                 ! NH4HSO4
2531          MOLALR(12) = PSI3                                 ! NAHSO4
2532          MOLALR(13) = PSI2                                 ! LC
2533 !
2534 ! *** NA-NH4-SO4-NO3-CL SYSTEM ; SULFATE RICH CASE ; FREE ACID
2535 !
2536       ELSE IF (SC.EQ.'J') THEN
2537          MOLALR(09) = MOLAL(3)                             ! NH4HSO4
2538          MOLALR(12) = MOLAL(2)                             ! NAHSO4
2539          MOLALR(07) = MOLAL(5)+MOLAL(6)-MOLAL(3)-MOLAL(2)  ! H2SO4
2540          MOLALR(07) = MAX(MOLALR(07),ZERO)
2541 !
2542 ! ======= REVERSE PROBLEMS ===========================================
2543 !
2544 ! *** NH4-SO4-NO3 SYSTEM ; SULFATE POOR CASE
2545 !
2546       ELSE IF (SC.EQ.'N') THEN
2547          MOLALR(4) = MOLAL(5) + MOLAL(6)          ! (NH4)2SO4
2548          AML5      = WAER(3)-2.D0*MOLALR(4)       ! "free" NH4
2549          MOLALR(5) = MAX(MIN(AML5,WAER(4)), ZERO) ! NH4NO3 = MIN("free", NO3)
2550 !
2551 ! *** NH4-SO4-NO3-NA-CL SYSTEM ; SULFATE POOR, SODIUM POOR CASE
2552 !
2553       ELSE IF (SC.EQ.'Q') THEN
2554          MOLALR(2) = PSI1                                  ! NA2SO4
2555          MOLALR(4) = PSI6                                  ! (NH4)2SO4
2556          MOLALR(5) = PSI5                                  ! NH4NO3
2557          MOLALR(6) = PSI4                                  ! NH4CL
2558 !
2559 ! *** NH4-SO4-NO3-NA-CL SYSTEM ; SULFATE POOR, SODIUM RICH CASE
2560 !
2561       ELSE IF (SC.EQ.'R') THEN
2562          MOLALR(1) = PSI3                                  ! NACL
2563          MOLALR(2) = PSI1                                  ! NA2SO4
2564          MOLALR(3) = PSI2                                  ! NANO3
2565          MOLALR(4) = ZERO                                  ! (NH4)2SO4
2566          MOLALR(5) = PSI5                                  ! NH4NO3
2567          MOLALR(6) = PSI4                                  ! NH4CL
2568 !
2569 ! *** UNKNOWN CASE
2570 !
2571       ELSE
2572          CALL PUSHERR (1001, ' ') ! FATAL ERROR: CASE NOT SUPPORTED
2573       ENDIF
2574 !
2575 ! *** CALCULATE WATER CONTENT ; ZSR CORRELATION ***********************
2576 !
2577       WATER = ZERO
2578       DO 10 I=1,NPAIR
2579          WATER = WATER + MOLALR(I)/M0(I)
2580 10    CONTINUE
2581       WATER = MAX(WATER, TINY)
2582 !
2583       RETURN
2584 !
2585 ! *** END OF SUBROUTINE CALCMR ******************************************
2586 !
2587       END
2588 !=======================================================================
2589 !
2590 ! *** ISORROPIA CODE
2591 ! *** SUBROUTINE CALCMDRH
2592 !
2593 !     THIS IS THE CASE WHERE THE RELATIVE HUMIDITY IS IN THE MUTUAL
2594 !     DRH REGION. THE SOLUTION IS ASSUMED TO BE THE SUM OF TWO WEIGHTED
2595 !     SOLUTIONS ; THE 'DRY' SOLUTION (SUBROUTINE DRYCASE) AND THE
2596 !     'SATURATED LIQUID' SOLUTION (SUBROUTINE LIQCASE).
2597 !
2598 ! *** COPYRIGHT 1996-2006, UNIVERSITY OF MIAMI, CARNEGIE MELLON UNIVERSITY,
2599 ! *** GEORGIA INSTITUTE OF TECHNOLOGY
2600 ! *** WRITTEN BY ATHANASIOS NENES
2601 ! *** UPDATED BY CHRISTOS FOUNTOUKIS
2602 !
2603 !=======================================================================
2604 !
2605       SUBROUTINE CALCMDRH (RHI, RHDRY, RHLIQ, DRYCASE, LIQCASE)
2606       USE ISRPIA
2607       IMPLICIT DOUBLE PRECISION (A-H, O-Z)
2608 !      INCLUDE 'ISRPIA.INC'
2609       EXTERNAL DRYCASE, LIQCASE
2610 !
2611 ! *** FIND WEIGHT FACTOR **********************************************
2612 !
2613       IF (WFTYP.EQ.0) THEN
2614          WF = ONE
2615       ELSEIF (WFTYP.EQ.1) THEN
2616          WF = 0.5D0
2617       ELSE
2618          WF = (RHLIQ-RHI)/(RHLIQ-RHDRY)
2619       ENDIF
2620       ONEMWF  = ONE - WF
2621 !
2622 ! *** FIND FIRST SECTION ; DRY ONE ************************************
2623 !
2624       CALL DRYCASE
2625       IF (ABS(ONEMWF).LE.1D-5) GOTO 200  ! DRY AEROSOL
2626 !
2627       CNH42SO = CNH42S4                  ! FIRST (DRY) SOLUTION
2628       CNH4HSO = CNH4HS4
2629       CLCO    = CLC
2630       CNH4N3O = CNH4NO3
2631       CNH4CLO = CNH4CL
2632       CNA2SO  = CNA2SO4
2633       CNAHSO  = CNAHSO4
2634       CNANO   = CNANO3
2635       CNACLO  = CNACL
2636       GNH3O   = GNH3
2637       GHNO3O  = GHNO3
2638       GHCLO   = GHCL
2639 !
2640 ! *** FIND SECOND SECTION ; DRY & LIQUID ******************************
2641 !
2642       CNH42S4 = ZERO
2643       CNH4HS4 = ZERO
2644       CLC     = ZERO
2645       CNH4NO3 = ZERO
2646       CNH4CL  = ZERO
2647       CNA2SO4 = ZERO
2648       CNAHSO4 = ZERO
2649       CNANO3  = ZERO
2650       CNACL   = ZERO
2651       GNH3    = ZERO
2652       GHNO3   = ZERO
2653       GHCL    = ZERO
2654       CALL LIQCASE                   ! SECOND (LIQUID) SOLUTION
2655 !
2656 ! *** ADJUST THINGS FOR THE CASE THAT THE LIQUID SUB PREDICTS DRY AEROSOL
2657 !
2658       IF (WATER.LE.TINY) THEN
2659          DO 100 I=1,NIONS
2660             MOLAL(I)= ZERO           ! Aqueous phase
2661   100    CONTINUE
2662          WATER   = ZERO
2663 !
2664          CNH42S4 = CNH42SO           ! Solid phase
2665          CNA2SO4 = CNA2SO
2666          CNAHSO4 = CNAHSO
2667          CNH4HS4 = CNH4HSO
2668          CLC     = CLCO
2669          CNH4NO3 = CNH4N3O
2670          CNANO3  = CNANO
2671          CNACL   = CNACLO
2672          CNH4CL  = CNH4CLO
2673 !
2674          GNH3    = GNH3O             ! Gas phase
2675          GHNO3   = GHNO3O
2676          GHCL    = GHCLO
2677 !
2678          GOTO 200
2679       ENDIF
2680 !
2681 ! *** FIND SALT DISSOLUTIONS BETWEEN DRY & LIQUID SOLUTIONS.
2682 !
2683       DAMSUL  = CNH42SO - CNH42S4
2684       DSOSUL  = CNA2SO  - CNA2SO4
2685       DAMBIS  = CNH4HSO - CNH4HS4
2686       DSOBIS  = CNAHSO  - CNAHSO4
2687       DLC     = CLCO    - CLC
2688       DAMNIT  = CNH4N3O - CNH4NO3
2689       DAMCHL  = CNH4CLO - CNH4CL
2690       DSONIT  = CNANO   - CNANO3
2691       DSOCHL  = CNACLO  - CNACL
2692 !
2693 ! *** FIND GAS DISSOLUTIONS BETWEEN DRY & LIQUID SOLUTIONS.
2694 !
2695       DAMG    = GNH3O   - GNH3
2696       DHAG    = GHCLO   - GHCL
2697       DNAG    = GHNO3O  - GHNO3
2698 !
2699 ! *** FIND SOLUTION AT MDRH BY WEIGHTING DRY & LIQUID SOLUTIONS.
2700 !
2701 !     LIQUID
2702 !
2703       MOLAL(1)= ONEMWF*MOLAL(1)                                 ! H+
2704       MOLAL(2)= ONEMWF*(2.D0*DSOSUL + DSOBIS + DSONIT + DSOCHL) ! NA+
2705       MOLAL(3)= ONEMWF*(2.D0*DAMSUL + DAMG   + DAMBIS + DAMCHL +   &
2706                         3.D0*DLC    + DAMNIT )                  ! NH4+
2707       MOLAL(4)= ONEMWF*(     DAMCHL + DSOCHL + DHAG)            ! CL-
2708       MOLAL(5)= ONEMWF*(     DAMSUL + DSOSUL + DLC - MOLAL(6))  ! SO4-- !VB 17 Sept 2001
2709       MOLAL(6)= ONEMWF*(   MOLAL(6) + DSOBIS + DAMBIS + DLC)    ! HSO4-
2710       MOLAL(7)= ONEMWF*(     DAMNIT + DSONIT + DNAG)            ! NO3-
2711       WATER   = ONEMWF*WATER
2712 !
2713 !     SOLID
2714 !
2715       CNH42S4 = WF*CNH42SO + ONEMWF*CNH42S4
2716       CNA2SO4 = WF*CNA2SO  + ONEMWF*CNA2SO4
2717       CNAHSO4 = WF*CNAHSO  + ONEMWF*CNAHSO4
2718       CNH4HS4 = WF*CNH4HSO + ONEMWF*CNH4HS4
2719       CLC     = WF*CLCO    + ONEMWF*CLC
2720       CNH4NO3 = WF*CNH4N3O + ONEMWF*CNH4NO3
2721       CNANO3  = WF*CNANO   + ONEMWF*CNANO3
2722       CNACL   = WF*CNACLO  + ONEMWF*CNACL
2723       CNH4CL  = WF*CNH4CLO + ONEMWF*CNH4CL
2724 !
2725 !     GAS
2726 !
2727       GNH3    = WF*GNH3O   + ONEMWF*GNH3
2728       GHNO3   = WF*GHNO3O  + ONEMWF*GHNO3
2729       GHCL    = WF*GHCLO   + ONEMWF*GHCL
2730 !
2731 ! *** RETURN POINT
2732 !
2733 200   RETURN
2734 !
2735 ! *** END OF SUBROUTINE CALCMDRH ****************************************
2736 !
2737       END
2738
2739
2740
2741
2742
2743
2744 !=======================================================================
2745 !
2746 ! *** ISORROPIA CODE
2747 ! *** SUBROUTINE CALCMDRP
2748 !
2749 !     THIS IS THE CASE WHERE THE RELATIVE HUMIDITY IS IN THE MUTUAL
2750 !     DRH REGION. THE SOLUTION IS ASSUMED TO BE THE SUM OF TWO WEIGHTED
2751 !     SOLUTIONS ; THE 'DRY' SOLUTION (SUBROUTINE DRYCASE) AND THE
2752 !     'SATURATED LIQUID' SOLUTION (SUBROUTINE LIQCASE).
2753 !
2754 ! *** COPYRIGHT 1996-2006, UNIVERSITY OF MIAMI, CARNEGIE MELLON UNIVERSITY,
2755 ! *** GEORGIA INSTITUTE OF TECHNOLOGY
2756 ! *** WRITTEN BY ATHANASIOS NENES
2757 ! *** UPDATED BY CHRISTOS FOUNTOUKIS
2758 !
2759 !=======================================================================
2760 !
2761       SUBROUTINE CALCMDRP (RHI, RHDRY, RHLIQ, DRYCASE, LIQCASE)
2762       USE ISRPIA
2763       IMPLICIT DOUBLE PRECISION (A-H, O-Z)
2764 !      INCLUDE 'ISRPIA.INC'
2765       EXTERNAL DRYCASE, LIQCASE
2766 !
2767 ! *** FIND WEIGHT FACTOR **********************************************
2768 !
2769       IF (WFTYP.EQ.0) THEN
2770          WF = ONE
2771       ELSEIF (WFTYP.EQ.1) THEN
2772          WF = 0.5D0
2773       ELSE
2774          WF = (RHLIQ-RHI)/(RHLIQ-RHDRY)
2775       ENDIF
2776       ONEMWF  = ONE - WF
2777 !
2778 ! *** FIND FIRST SECTION ; DRY ONE ************************************
2779 !
2780       CALL DRYCASE
2781       IF (ABS(ONEMWF).LE.1D-5) GOTO 200  ! DRY AEROSOL
2782 !
2783       CNH42SO = CNH42S4              ! FIRST (DRY) SOLUTION
2784       CNH4HSO = CNH4HS4
2785       CLCO    = CLC
2786       CNH4N3O = CNH4NO3
2787       CNH4CLO = CNH4CL
2788       CNA2SO  = CNA2SO4
2789       CNAHSO  = CNAHSO4
2790       CNANO   = CNANO3
2791       CNACLO  = CNACL
2792 !
2793 ! *** FIND SECOND SECTION ; DRY & LIQUID ******************************
2794 !
2795       CNH42S4 = ZERO
2796       CNH4HS4 = ZERO
2797       CLC     = ZERO
2798       CNH4NO3 = ZERO
2799       CNH4CL  = ZERO
2800       CNA2SO4 = ZERO
2801       CNAHSO4 = ZERO
2802       CNANO3  = ZERO
2803       CNACL   = ZERO
2804       GNH3    = ZERO
2805       GHNO3   = ZERO
2806       GHCL    = ZERO
2807       CALL LIQCASE                   ! SECOND (LIQUID) SOLUTION
2808 !
2809 ! *** ADJUST THINGS FOR THE CASE THAT THE LIQUID SUB PREDICTS DRY AEROSOL
2810 !
2811       IF (WATER.LE.TINY) THEN
2812          WATER = ZERO
2813          DO 100 I=1,NIONS
2814             MOLAL(I)= ZERO
2815  100     CONTINUE
2816          CALL DRYCASE
2817          GOTO 200
2818       ENDIF
2819 !
2820 ! *** FIND SALT DISSOLUTIONS BETWEEN DRY & LIQUID SOLUTIONS.
2821 !
2822       DAMBIS  = CNH4HSO - CNH4HS4
2823       DSOBIS  = CNAHSO  - CNAHSO4
2824       DLC     = CLCO    - CLC
2825 !
2826 ! *** FIND SOLUTION AT MDRH BY WEIGHTING DRY & LIQUID SOLUTIONS.
2827 !
2828 ! *** SOLID
2829 !
2830       CNH42S4 = WF*CNH42SO + ONEMWF*CNH42S4
2831       CNA2SO4 = WF*CNA2SO  + ONEMWF*CNA2SO4
2832       CNAHSO4 = WF*CNAHSO  + ONEMWF*CNAHSO4
2833       CNH4HS4 = WF*CNH4HSO + ONEMWF*CNH4HS4
2834       CLC     = WF*CLCO    + ONEMWF*CLC
2835       CNH4NO3 = WF*CNH4N3O + ONEMWF*CNH4NO3
2836       CNANO3  = WF*CNANO   + ONEMWF*CNANO3
2837       CNACL   = WF*CNACLO  + ONEMWF*CNACL
2838       CNH4CL  = WF*CNH4CLO + ONEMWF*CNH4CL
2839 !
2840 ! *** LIQUID
2841 !
2842       WATER   = ONEMWF*WATER
2843 !
2844       MOLAL(2)= WAER(1) - 2.D0*CNA2SO4 - CNAHSO4 - CNANO3 -        &
2845                                CNACL                            ! NA+
2846       MOLAL(3)= WAER(3) - 2.D0*CNH42S4 - CNH4HS4 - CNH4CL -    &
2847                           3.D0*CLC     - CNH4NO3                ! NH4+
2848       MOLAL(4)= WAER(5) - CNACL - CNH4CL                        ! CL-
2849       MOLAL(7)= WAER(4) - CNANO3 - CNH4NO3                      ! NO3-
2850       MOLAL(6)= ONEMWF*(MOLAL(6) + DSOBIS + DAMBIS + DLC)       ! HSO4-
2851       MOLAL(5)= WAER(2) - MOLAL(6) - CLC - CNH42S4 - CNA2SO4    ! SO4--
2852 !
2853       A8      = XK1*WATER/GAMA(7)*(GAMA(8)/GAMA(7))**2.
2854       IF (MOLAL(5).LE.TINY) THEN
2855          HIEQ = SQRT(XKW *RH*WATER*WATER)  ! Neutral solution
2856       ELSE
2857          HIEQ = A8*MOLAL(6)/MOLAL(5)
2858       ENDIF
2859       HIEN    = MOLAL(4) + MOLAL(7) + MOLAL(6) + 2.D0*MOLAL(5) -   &
2860                 MOLAL(2) - MOLAL(3)
2861       MOLAL(1)= MAX (HIEQ, HIEN)                                ! H+
2862 !
2863 ! *** GAS (ACTIVITY COEFS FROM LIQUID SOLUTION)
2864 !
2865       A2      = (XK2/XKW)*R*TEMP*(GAMA(10)/GAMA(5))**2. ! NH3  <==> NH4+
2866       A3      = XK4 *R*TEMP*(WATER/GAMA(10))**2.        ! HNO3 <==> NO3-
2867       A4      = XK3 *R*TEMP*(WATER/GAMA(11))**2.        ! HCL  <==> CL-
2868 !
2869       GNH3    = MOLAL(3)/MAX(MOLAL(1),TINY)/A2
2870       GHNO3   = MOLAL(1)*MOLAL(7)/A3
2871       GHCL    = MOLAL(1)*MOLAL(4)/A4
2872 !
2873 200   RETURN
2874 !
2875 ! *** END OF SUBROUTINE CALCMDRP ****************************************
2876 !
2877       END
2878 !=======================================================================
2879 !
2880 ! *** ISORROPIA CODE
2881 ! *** SUBROUTINE CALCHS4
2882 ! *** THIS SUBROUTINE CALCULATES THE HSO4 GENERATED FROM (H,SO4).
2883 !
2884 ! *** COPYRIGHT 1996-2006, UNIVERSITY OF MIAMI, CARNEGIE MELLON UNIVERSITY,
2885 ! *** GEORGIA INSTITUTE OF TECHNOLOGY
2886 ! *** WRITTEN BY ATHANASIOS NENES
2887 ! *** UPDATED BY CHRISTOS FOUNTOUKIS
2888 !
2889 !=======================================================================
2890 !
2891       SUBROUTINE CALCHS4 (HI, SO4I, HSO4I, DELTA)
2892       USE ISRPIA
2893       IMPLICIT DOUBLE PRECISION (A-H, O-Z)
2894 !      INCLUDE 'ISRPIA.INC'
2895 !C      CHARACTER ERRINF*40
2896 !
2897 ! *** IF TOO LITTLE WATER, DONT SOLVE
2898 !
2899       IF (WATER.LE.1d1*TINY) THEN
2900          DELTA = ZERO
2901          RETURN
2902       ENDIF
2903 !
2904 ! *** CALCULATE HSO4 SPECIATION *****************************************
2905 !
2906       A8 = XK1*WATER/GAMA(7)*(GAMA(8)/GAMA(7))**2.
2907 !
2908       BB =-(HI + SO4I + A8)
2909       CC = HI*SO4I - HSO4I*A8
2910       DD = BB*BB - 4.D0*CC
2911 !
2912       IF (DD.GE.ZERO) THEN
2913          SQDD   = SQRT(DD)
2914          DELTA1 = 0.5*(-BB + SQDD)
2915          DELTA2 = 0.5*(-BB - SQDD)
2916          IF (HSO4I.LE.TINY) THEN
2917             DELTA = DELTA2
2918          ELSEIF( HI*SO4I .GE. A8*HSO4I ) THEN
2919             DELTA = DELTA2
2920          ELSEIF( HI*SO4I .LT. A8*HSO4I ) THEN
2921             DELTA = DELTA1
2922          ELSE
2923             DELTA = ZERO
2924          ENDIF
2925       ELSE
2926          DELTA  = ZERO
2927       ENDIF
2928 !CC
2929 !CC *** COMPARE DELTA TO TOTAL H+ ; ESTIMATE EFFECT OF HSO4 ***************
2930 !CC
2931 !C      HYD = MAX(HI, MOLAL(1))
2932 !C      IF (HYD.GT.TINY) THEN
2933 !C         IF (DELTA/HYD.GT.0.1d0) THEN
2934 !C            WRITE (ERRINF,'(1PE10.3)') DELTA/HYD*100.0
2935 !C            CALL PUSHERR (0020, ERRINF)
2936 !C         ENDIF
2937 !C      ENDIF
2938 !
2939       RETURN
2940 !
2941 ! *** END OF SUBROUTINE CALCHS4 *****************************************
2942 !
2943       END
2944
2945
2946 !=======================================================================
2947 !
2948 ! *** ISORROPIA CODE
2949 ! *** SUBROUTINE CALCPH
2950 !
2951 ! *** COPYRIGHT 1996-2006, UNIVERSITY OF MIAMI, CARNEGIE MELLON UNIVERSITY,
2952 ! *** GEORGIA INSTITUTE OF TECHNOLOGY
2953 ! *** WRITTEN BY ATHANASIOS NENES
2954 ! *** UPDATED BY CHRISTOS FOUNTOUKIS
2955 !
2956 !=======================================================================
2957 !
2958       SUBROUTINE CALCPH (GG, HI, OHI)
2959       USE ISRPIA
2960       IMPLICIT DOUBLE PRECISION (A-H, O-Z)
2961 !      INCLUDE 'ISRPIA.INC'
2962 !
2963       AKW  = XKW *RH*WATER*WATER
2964       CN   = SQRT(AKW)
2965 !
2966 ! *** GG = (negative charge) - (positive charge)
2967 !
2968       IF (GG.GT.TINY) THEN                        ! H+ in excess
2969          BB =-GG
2970          CC =-AKW
2971          DD = BB*BB - 4.D0*CC
2972          HI = MAX(0.5D0*(-BB + SQRT(DD)),CN)
2973          OHI= AKW/HI
2974       ELSE                                        ! OH- in excess
2975          BB = GG
2976          CC =-AKW
2977          DD = BB*BB - 4.D0*CC
2978          OHI= MAX(0.5D0*(-BB + SQRT(DD)),CN)
2979          HI = AKW/OHI
2980       ENDIF
2981 !
2982       RETURN
2983 !
2984 ! *** END OF SUBROUTINE CALCPH ******************************************
2985 !
2986       END
2987
2988 !=======================================================================
2989 !
2990 ! *** ISORROPIA CODE
2991 ! *** SUBROUTINE CALCACT
2992 ! *** CALCULATES MULTI-COMPONENET ACTIVITY COEFFICIENTS FROM BROMLEYS
2993 !     METHOD. THE BINARY ACTIVITY COEFFICIENTS ARE CALCULATED BY
2994 !     KUSIK-MEISNER RELATION (SUBROUTINE KMTAB or SUBROUTINE KMFUL).
2995 !
2996 ! *** COPYRIGHT 1996-2006, UNIVERSITY OF MIAMI, CARNEGIE MELLON UNIVERSITY,
2997 ! *** GEORGIA INSTITUTE OF TECHNOLOGY
2998 ! *** WRITTEN BY ATHANASIOS NENES
2999 ! *** UPDATED BY CHRISTOS FOUNTOUKIS
3000 !
3001 !=======================================================================
3002 !
3003       SUBROUTINE CALCACT
3004       USE ISRPIA
3005       IMPLICIT DOUBLE PRECISION (A-H, O-Z)
3006 !      INCLUDE 'ISRPIA.INC'
3007 !
3008       REAL EX10, URF
3009       REAL G0(3,4),ZPL,ZMI,AGAMA,SION,H,CH,F1(3),F2(4)
3010       DOUBLE PRECISION MPL, XIJ, YJI
3011 !      PARAMETER (URF=0.5)
3012       PARAMETER (LN10=2.30258509299404568402D0)
3013 !
3014       G(I,J)= (F1(I)/Z(I) + F2(J)/Z(J+3)) / (Z(I)+Z(J+3)) - H
3015 !
3016 ! *** SAVE ACTIVITIES IN OLD ARRAY *************************************
3017 !
3018       IF (FRST) THEN               ! Outer loop
3019          DO 10 I=1,NPAIR
3020             GAMOU(I) = GAMA(I)
3021 10       CONTINUE
3022       ENDIF
3023 !
3024       DO 20 I=1,NPAIR              ! Inner loop
3025          GAMIN(I) = GAMA(I)
3026 20    CONTINUE
3027 !
3028 ! *** CALCULATE IONIC ACTIVITY OF SOLUTION *****************************
3029 !
3030       IONIC=0.0
3031       DO 30 I=1,NIONS
3032          IONIC=IONIC + MOLAL(I)*Z(I)*Z(I)
3033 30    CONTINUE
3034       IONIC = MAX(MIN(0.5*IONIC/WATER,500.d0), TINY)
3035 !
3036 ! *** CALCULATE BINARY ACTIVITY COEFFICIENTS ***************************
3037 !
3038 !  G0(1,1)=G11;G0(1,2)=G07;G0(1,3)=G08;G0(1,4)=G10;G0(2,1)=G01;G0(2,2)=G02
3039 !  G0(2,3)=G12;G0(2,4)=G03;G0(3,1)=G06;G0(3,2)=G04;G0(3,3)=G09;G0(3,4)=G05
3040 !
3041       IF (IACALC.EQ.0) THEN              ! K.M.; FULL
3042          CALL KMFUL (IONIC, SNGL(TEMP),G0(2,1),G0(2,2),G0(2,4),   &
3043                      G0(3,2),G0(3,4),G0(3,1),G0(1,2),G0(1,3),G0(3,3),   &
3044                      G0(1,4),G0(1,1),G0(2,3))
3045       ELSE                               ! K.M.; TABULATED
3046          CALL KMTAB (IONIC, SNGL(TEMP),G0(2,1),G0(2,2),G0(2,4),   &
3047                      G0(3,2),G0(3,4),G0(3,1),G0(1,2),G0(1,3),G0(3,3),   &
3048                      G0(1,4),G0(1,1),G0(2,3))
3049       ENDIF
3050 !
3051 ! *** CALCULATE MULTICOMPONENT ACTIVITY COEFFICIENTS *******************
3052 !
3053       AGAMA = 0.511*(298.0/TEMP)**1.5    ! Debye Huckel const. at T
3054       SION  = SQRT(IONIC)
3055       H     = AGAMA*SION/(1+SION)
3056 !
3057       DO 100 I=1,3
3058          F1(I)=0.0
3059          F2(I)=0.0
3060 100   CONTINUE
3061       F2(4)=0.0
3062 !
3063       DO 110 I=1,3
3064          ZPL = Z(I)
3065          MPL = MOLAL(I)/WATER
3066          DO 110 J=1,4
3067             ZMI   = Z(J+3)
3068             CH    = 0.25*(ZPL+ZMI)*(ZPL+ZMI)/IONIC
3069             XIJ   = CH*MPL
3070             YJI   = CH*MOLAL(J+3)/WATER
3071             F1(I) = F1(I) + SNGL(YJI*(G0(I,J) + ZPL*ZMI*H))
3072             F2(J) = F2(J) + SNGL(XIJ*(G0(I,J) + ZPL*ZMI*H))
3073 110   CONTINUE
3074 !
3075 ! *** LOG10 OF ACTIVITY COEFFICIENTS ***********************************
3076 !
3077       GAMA(01) = G(2,1)*ZZ(01)                     ! NACL
3078       GAMA(02) = G(2,2)*ZZ(02)                     ! NA2SO4
3079       GAMA(03) = G(2,4)*ZZ(03)                     ! NANO3
3080       GAMA(04) = G(3,2)*ZZ(04)                     ! (NH4)2SO4
3081       GAMA(05) = G(3,4)*ZZ(05)                     ! NH4NO3
3082       GAMA(06) = G(3,1)*ZZ(06)                     ! NH4CL
3083       GAMA(07) = G(1,2)*ZZ(07)                     ! 2H-SO4
3084       GAMA(08) = G(1,3)*ZZ(08)                     ! H-HSO4
3085       GAMA(09) = G(3,3)*ZZ(09)                     ! NH4HSO4
3086       GAMA(10) = G(1,4)*ZZ(10)                     ! HNO3
3087       GAMA(11) = G(1,1)*ZZ(11)                     ! HCL
3088       GAMA(12) = G(2,3)*ZZ(12)                     ! NAHSO4
3089       GAMA(13) = 0.20*(3.0*GAMA(04)+2.0*GAMA(09))  ! LC ; SCAPE
3090 !C      GAMA(13) = 0.50*(GAMA(04)+GAMA(09))          ! LC ; SEQUILIB
3091 !C      GAMA(13) = 0.25*(3.0*GAMA(04)+GAMA(07))      ! LC ; AIM
3092 !
3093 ! *** CONVERT LOG (GAMA) COEFFICIENTS TO GAMA **************************
3094 !
3095       DO 200 I=1,NPAIR
3096          GAMA(I)=MAX(-11.0d0, MIN(GAMA(I),11.0d0) ) ! F77 LIBRARY ROUTINE
3097 !         GAMA(I)=10.0**GAMA(I)
3098          GAMA(I)=EXP(LN10*GAMA(I))
3099 !C         GAMA(I)=EX10(SNGL(GAMA(I)), 5.0)    ! CUTOFF SET TO [-5,5]
3100 !         GAMA(I) = GAMIN(I)*(1.0-URF) + URF*GAMA(I)  ! Under-relax GAMA's
3101   200 CONTINUE
3102 !
3103 ! *** SETUP ACTIVITY CALCULATION FLAGS *********************************
3104 !
3105 ! OUTER CALCULATION LOOP ; ONLY IF FRST=.TRUE.
3106 !
3107       IF (FRST) THEN
3108          ERROU = ZERO                    ! CONVERGENCE CRITERION
3109          DO 210 I=1,NPAIR
3110             ERROU=MAX(ERROU, ABS((GAMOU(I)-GAMA(I))/GAMOU(I)))
3111 210      CONTINUE
3112          CALAOU = ERROU .GE. EPSACT      ! SETUP FLAGS
3113          FRST   =.FALSE.
3114       ENDIF
3115 !
3116 ! INNER CALCULATION LOOP ; ALWAYS
3117 !
3118       ERRIN = ZERO                       ! CONVERGENCE CRITERION
3119       DO 220 I=1,NPAIR
3120          ERRIN = MAX (ERRIN, ABS((GAMIN(I)-GAMA(I))/GAMIN(I)))
3121 220   CONTINUE
3122       CALAIN = ERRIN .GE. EPSACT
3123 !
3124       ICLACT = ICLACT + 1                ! Increment ACTIVITY call counter
3125 !
3126 ! *** END OF SUBROUTINE ACTIVITY ****************************************
3127 !
3128       RETURN
3129       END
3130
3131
3132 !=======================================================================
3133 !
3134 ! *** ISORROPIA CODE
3135 ! *** SUBROUTINE RSTGAM
3136 ! *** RESETS ACTIVITY COEFFICIENT ARRAYS TO DEFAULT VALUE OF 0.1
3137 !
3138 ! *** COPYRIGHT 1996-2006, UNIVERSITY OF MIAMI, CARNEGIE MELLON UNIVERSITY,
3139 ! *** GEORGIA INSTITUTE OF TECHNOLOGY
3140 ! *** WRITTEN BY ATHANASIOS NENES
3141 ! *** UPDATED BY CHRISTOS FOUNTOUKIS
3142 !
3143 !=======================================================================
3144 !
3145       SUBROUTINE RSTGAM
3146       USE ISRPIA
3147       IMPLICIT DOUBLE PRECISION (A-H, O-Z)
3148 !      INCLUDE 'ISRPIA.INC'
3149 !
3150       DO 10 I=1, NPAIR
3151          GAMA(I) = 0.1
3152 10    CONTINUE
3153 !
3154 ! *** END OF SUBROUTINE RSTGAM ******************************************
3155 !
3156       RETURN
3157       END
3158 !=======================================================================
3159 !
3160 ! *** ISORROPIA CODE
3161 ! *** SUBROUTINE KMFUL
3162 ! *** CALCULATES BINARY ACTIVITY COEFFICIENTS BY KUSIK-MEISSNER METHOD.
3163 !
3164 ! *** COPYRIGHT 1996-2006, UNIVERSITY OF MIAMI, CARNEGIE MELLON UNIVERSITY,
3165 ! *** GEORGIA INSTITUTE OF TECHNOLOGY
3166 ! *** WRITTEN BY ATHANASIOS NENES
3167 ! *** UPDATED BY CHRISTOS FOUNTOUKIS
3168 !
3169 !=======================================================================
3170 !
3171       SUBROUTINE KMFUL (IONIC,TEMP,G01,G02,G03,G04,G05,G06,G07,G08,G09,   &
3172                         G10,G11,G12)
3173       REAL Ionic, TEMP
3174       DATA Z01,Z02,Z03,Z04,Z05,Z06,Z07,Z08,Z10,Z11   &
3175           /1,  2,  1,  2,  1,  1,  2,  1,  1,  1/
3176 !
3177       SION = SQRT(IONIC)
3178 !
3179 ! *** Coefficients at 25 oC
3180 !
3181       CALL MKBI(2.230, IONIC, SION, Z01, G01)
3182       CALL MKBI(-0.19, IONIC, SION, Z02, G02)
3183       CALL MKBI(-0.39, IONIC, SION, Z03, G03)
3184       CALL MKBI(-0.25, IONIC, SION, Z04, G04)
3185       CALL MKBI(-1.15, IONIC, SION, Z05, G05)
3186       CALL MKBI(0.820, IONIC, SION, Z06, G06)
3187       CALL MKBI(-.100, IONIC, SION, Z07, G07)
3188       CALL MKBI(8.000, IONIC, SION, Z08, G08)
3189       CALL MKBI(2.600, IONIC, SION, Z10, G10)
3190       CALL MKBI(6.000, IONIC, SION, Z11, G11)
3191 !
3192 ! *** Correct for T other than 298 K
3193 !
3194       TI  = TEMP-273.0
3195       TC  = TI-25.0
3196       IF (ABS(TC) .GT. 1.0) THEN
3197          CF1 = 1.125-0.005*TI
3198          CF2 = (0.125-0.005*TI)*(0.039*IONIC**0.92-0.41*SION/(1.+SION))
3199          G01 = CF1*G01 - CF2*Z01
3200          G02 = CF1*G02 - CF2*Z02
3201          G03 = CF1*G03 - CF2*Z03
3202          G04 = CF1*G04 - CF2*Z04
3203          G05 = CF1*G05 - CF2*Z05
3204          G06 = CF1*G06 - CF2*Z06
3205          G07 = CF1*G07 - CF2*Z07
3206          G08 = CF1*G08 - CF2*Z08
3207          G10 = CF1*G10 - CF2*Z10
3208          G11 = CF1*G11 - CF2*Z11
3209       ENDIF
3210 !
3211       G09 = G06 + G08 - G11
3212       G12 = G01 + G08 - G11
3213 !
3214 ! *** Return point ; End of subroutine
3215 !
3216       RETURN
3217       END
3218
3219
3220 !=======================================================================
3221 !
3222 ! *** ISORROPIA CODE
3223 ! *** SUBROUTINE MKBI
3224 ! *** CALCULATES BINARY ACTIVITY COEFFICIENTS BY KUSIK-MEISSNER METHOD.
3225 !
3226 ! *** COPYRIGHT 1996-2006, UNIVERSITY OF MIAMI, CARNEGIE MELLON UNIVERSITY,
3227 ! *** GEORGIA INSTITUTE OF TECHNOLOGY
3228 ! *** WRITTEN BY ATHANASIOS NENES
3229 ! *** UPDATED BY CHRISTOS FOUNTOUKIS
3230 !
3231 !=======================================================================
3232 !
3233       SUBROUTINE MKBI(Q,IONIC,SION,ZIP,BI)
3234 !
3235       REAL IONIC
3236 !
3237       B=.75-.065*Q
3238       C= 1.0
3239       IF (IONIC.LT.6.0) C=1.+.055*Q*EXP(-.023*IONIC*IONIC*IONIC)
3240       XX=-0.5107*SION/(1.+C*SION)
3241       BI=(1.+B*(1.+.1*IONIC)**Q-B)
3242       BI=ZIP*ALOG10(BI) + ZIP*XX
3243 !
3244       RETURN
3245       END
3246 !=======================================================================
3247 !
3248 ! *** ISORROPIA CODE
3249 ! *** SUBROUTINE KMTAB
3250 ! *** CALCULATES BINARY ACTIVITY COEFFICIENTS BY KUSIK-MEISSNER METHOD.
3251 !     THE COMPUTATIONS HAVE BEEN PERFORMED AND THE RESULTS ARE STORED IN
3252 !     LOOKUP TABLES. THE IONIC ACTIVITY 'IONIC' IS INPUT, AND THE ARRAY
3253 !     'BINARR' IS RETURNED WITH THE BINARY COEFFICIENTS.
3254 !
3255 ! *** COPYRIGHT 1996-2006, UNIVERSITY OF MIAMI, CARNEGIE MELLON UNIVERSITY,
3256 ! *** GEORGIA INSTITUTE OF TECHNOLOGY
3257 ! *** WRITTEN BY ATHANASIOS NENES
3258 ! *** UPDATED BY CHRISTOS FOUNTOUKIS
3259 !
3260 !=======================================================================
3261 !
3262       SUBROUTINE KMTAB (IN,TEMP,G01,G02,G03,G04,G05,G06,G07,G08,G09,G10,   &
3263                                 G11,G12)
3264       REAL IN, Temp
3265 !
3266 ! *** Find temperature range
3267 !
3268       IND = NINT((TEMP-198.0)/25.0) + 1
3269       IND = MIN(MAX(IND,1),6)
3270 !
3271 ! *** Call appropriate routine
3272 !
3273       IF (IND.EQ.1) THEN
3274          CALL KM198 (IN,G01,G02,G03,G04,G05,G06,G07,G08,G09,G10,G11,G12)
3275       ELSEIF (IND.EQ.2) THEN
3276          CALL KM223 (IN,G01,G02,G03,G04,G05,G06,G07,G08,G09,G10,G11,G12)
3277       ELSEIF (IND.EQ.3) THEN
3278          CALL KM248 (IN,G01,G02,G03,G04,G05,G06,G07,G08,G09,G10,G11,G12)
3279       ELSEIF (IND.EQ.4) THEN
3280          CALL KM273 (IN,G01,G02,G03,G04,G05,G06,G07,G08,G09,G10,G11,G12)
3281       ELSEIF (IND.EQ.5) THEN
3282          CALL KM298 (IN,G01,G02,G03,G04,G05,G06,G07,G08,G09,G10,G11,G12)
3283       ELSE
3284          CALL KM323 (IN,G01,G02,G03,G04,G05,G06,G07,G08,G09,G10,G11,G12)
3285       ENDIF
3286 !
3287 ! *** Return point; End of subroutine
3288 !
3289       RETURN
3290       END
3291
3292
3293       INTEGER FUNCTION IBACPOS(IN)
3294 !
3295 !     Compute the index in the binary activity coefficient array
3296 !     based on the input ionic strength.
3297 !
3298 !     Chris Nolte, 6/16/05
3299 !
3300       implicit none
3301       real IN
3302       IF (IN .LE. 0.300000E+02) THEN
3303          ibacpos = MIN(NINT( 0.200000E+02*IN) + 1, 600)
3304       ELSE
3305          ibacpos =   600+NINT( 0.200000E+01*IN- 0.600000E+02)
3306       ENDIF
3307       ibacpos = min(ibacpos, 741)
3308       return
3309       end
3310
3311 !=======================================================================
3312 !
3313 ! *** ISORROPIA CODE
3314 ! *** SUBROUTINE KM198
3315 ! *** CALCULATES BINARY ACTIVITY COEFFICIENTS BY KUSIK-MEISSNER METHOD.
3316 !     THE COMPUTATIONS HAVE BEEN PERFORMED AND THE RESULTS ARE STORED IN
3317 !     LOOKUP TABLES. THE IONIC ACTIVITY 'IN' IS INPUT, AND THE ARRAY
3318 !     'BINARR' IS RETURNED WITH THE BINARY COEFFICIENTS.
3319 !
3320 !     TEMPERATURE IS 198K
3321 !
3322 ! *** COPYRIGHT 1996-2006, UNIVERSITY OF MIAMI, CARNEGIE MELLON UNIVERSITY,
3323 ! *** GEORGIA INSTITUTE OF TECHNOLOGY
3324 ! *** WRITTEN BY ATHANASIOS NENES
3325 ! *** UPDATED BY CHRISTOS FOUNTOUKIS
3326 !
3327 !=======================================================================
3328 !
3329       SUBROUTINE KM198 (IN,G01,G02,G03,G04,G05,G06,G07,G08,G09,G10,   &
3330                            G11,G12)
3331        USE KMC198
3332 !
3333 ! *** Common block definition
3334 !
3335 !      COMMON /KMC198/   &
3336 !      BNC01M(  741),BNC02M(  741),BNC03M(  741),BNC04M(  741),   &
3337 !      BNC05M(  741),BNC06M(  741),BNC07M(  741),BNC08M(  741),   &
3338 !      BNC09M(  741),BNC10M(  741),BNC11M(  741),BNC12M(  741),   &
3339 !      BNC13M(  741)
3340       REAL IN
3341 !
3342 ! *** Find position in arrays for binary activity coefficients
3343 !
3344       ipos = ibacpos(IN)
3345 !
3346 ! *** Assign values to return array
3347 !
3348       G01 = BNC01M(ipos)
3349       G02 = BNC02M(ipos)
3350       G03 = BNC03M(ipos)
3351       G04 = BNC04M(ipos)
3352       G05 = BNC05M(ipos)
3353       G06 = BNC06M(ipos)
3354       G07 = BNC07M(ipos)
3355       G08 = BNC08M(ipos)
3356       G09 = BNC09M(ipos)
3357       G10 = BNC10M(ipos)
3358       G11 = BNC11M(ipos)
3359       G12 = BNC12M(ipos)
3360 !
3361 ! *** Return point ; End of subroutine
3362 !
3363       RETURN
3364       END
3365
3366
3367 !=======================================================================
3368 !
3369 ! *** ISORROPIA CODE
3370 ! *** SUBROUTINE KM223
3371 ! *** CALCULATES BINARY ACTIVITY COEFFICIENTS BY KUSIK-MEISSNER METHOD.
3372 !     THE COMPUTATIONS HAVE BEEN PERFORMED AND THE RESULTS ARE STORED IN
3373 !     LOOKUP TABLES. THE IONIC ACTIVITY 'IN' IS INPUT, AND THE ARRAY
3374 !     'BINARR' IS RETURNED WITH THE BINARY COEFFICIENTS.
3375 !
3376 !     TEMPERATURE IS 223K
3377 !
3378 ! *** COPYRIGHT 1996-2006, UNIVERSITY OF MIAMI, CARNEGIE MELLON UNIVERSITY,
3379 ! *** GEORGIA INSTITUTE OF TECHNOLOGY
3380 ! *** WRITTEN BY ATHANASIOS NENES
3381 ! *** UPDATED BY CHRISTOS FOUNTOUKIS
3382 !
3383 !=======================================================================
3384 !
3385       SUBROUTINE KM223 (IN,G01,G02,G03,G04,G05,G06,G07,G08,G09,G10,   &
3386                            G11,G12)
3387       USE KMC223
3388 !
3389 ! *** Common block definition
3390 !
3391 !      COMMON /KMC223/   &
3392 !      BNC01M(  741),BNC02M(  741),BNC03M(  741),BNC04M(  741),   &
3393 !      BNC05M(  741),BNC06M(  741),BNC07M(  741),BNC08M(  741),   &
3394 !      BNC09M(  741),BNC10M(  741),BNC11M(  741),BNC12M(  741),   &
3395 !      BNC13M(  741)
3396       REAL IN
3397 !
3398 ! *** Find position in arrays for binary activity coefficients
3399 !
3400       ipos = ibacpos(IN)
3401 !
3402 ! *** Assign values to return array
3403 !
3404       G01 = BNC01M(ipos)
3405       G02 = BNC02M(ipos)
3406       G03 = BNC03M(ipos)
3407       G04 = BNC04M(ipos)
3408       G05 = BNC05M(ipos)
3409       G06 = BNC06M(ipos)
3410       G07 = BNC07M(ipos)
3411       G08 = BNC08M(ipos)
3412       G09 = BNC09M(ipos)
3413       G10 = BNC10M(ipos)
3414       G11 = BNC11M(ipos)
3415       G12 = BNC12M(ipos)
3416 !
3417 ! *** Return point ; End of subroutine
3418 !
3419       RETURN
3420       END
3421
3422
3423
3424 !=======================================================================
3425 !
3426 ! *** ISORROPIA CODE
3427 ! *** SUBROUTINE KM248
3428 ! *** CALCULATES BINARY ACTIVITY COEFFICIENTS BY KUSIK-MEISSNER METHOD.
3429 !     THE COMPUTATIONS HAVE BEEN PERFORMED AND THE RESULTS ARE STORED IN
3430 !     LOOKUP TABLES. THE IONIC ACTIVITY 'IN' IS INPUT, AND THE ARRAY
3431 !     'BINARR' IS RETURNED WITH THE BINARY COEFFICIENTS.
3432 !
3433 !     TEMPERATURE IS 248K
3434 !
3435 ! *** COPYRIGHT 1996-2006, UNIVERSITY OF MIAMI, CARNEGIE MELLON UNIVERSITY,
3436 ! *** GEORGIA INSTITUTE OF TECHNOLOGY
3437 ! *** WRITTEN BY ATHANASIOS NENES
3438 ! *** UPDATED BY CHRISTOS FOUNTOUKIS
3439 !
3440 !=======================================================================
3441 !
3442       SUBROUTINE KM248 (IN,G01,G02,G03,G04,G05,G06,G07,G08,G09,G10,   &
3443                            G11,G12)
3444        USE KMC248
3445 !
3446 ! *** Common block definition
3447 !
3448 !      COMMON /KMC248/   &
3449 !      BNC01M(  741),BNC02M(  741),BNC03M(  741),BNC04M(  741),   &
3450 !      BNC05M(  741),BNC06M(  741),BNC07M(  741),BNC08M(  741),   &
3451 !      BNC09M(  741),BNC10M(  741),BNC11M(  741),BNC12M(  741),   &
3452 !      BNC13M(  741)
3453       REAL IN
3454 !
3455 ! *** Find position in arrays for binary activity coefficients
3456 !
3457       ipos = ibacpos(IN)
3458 !
3459 ! *** Assign values to return array
3460 !
3461       G01 = BNC01M(ipos)
3462       G02 = BNC02M(ipos)
3463       G03 = BNC03M(ipos)
3464       G04 = BNC04M(ipos)
3465       G05 = BNC05M(ipos)
3466       G06 = BNC06M(ipos)
3467       G07 = BNC07M(ipos)
3468       G08 = BNC08M(ipos)
3469       G09 = BNC09M(ipos)
3470       G10 = BNC10M(ipos)
3471       G11 = BNC11M(ipos)
3472       G12 = BNC12M(ipos)
3473 !
3474 ! *** Return point ; End of subroutine
3475 !
3476       RETURN
3477       END
3478
3479
3480 !=======================================================================
3481 !
3482 ! *** ISORROPIA CODE
3483 ! *** SUBROUTINE KM273
3484 ! *** CALCULATES BINARY ACTIVITY COEFFICIENTS BY KUSIK-MEISSNER METHOD.
3485 !     THE COMPUTATIONS HAVE BEEN PERFORMED AND THE RESULTS ARE STORED IN
3486 !     LOOKUP TABLES. THE IONIC ACTIVITY 'IN' IS INPUT, AND THE ARRAY
3487 !     'BINARR' IS RETURNED WITH THE BINARY COEFFICIENTS.
3488 !
3489 !     TEMPERATURE IS 273K
3490 !
3491 ! *** COPYRIGHT 1996-2006, UNIVERSITY OF MIAMI, CARNEGIE MELLON UNIVERSITY,
3492 ! *** GEORGIA INSTITUTE OF TECHNOLOGY
3493 ! *** WRITTEN BY ATHANASIOS NENES
3494 ! *** UPDATED BY CHRISTOS FOUNTOUKIS
3495 !
3496 !=======================================================================
3497 !
3498       SUBROUTINE KM273 (IN,G01,G02,G03,G04,G05,G06,G07,G08,G09,G10,   &
3499                            G11,G12)
3500       USE KMC273
3501 !
3502 ! *** Common block definition
3503 !
3504 !      COMMON /KMC273/   &
3505 !      BNC01M(  741),BNC02M(  741),BNC03M(  741),BNC04M(  741),   &
3506 !      BNC05M(  741),BNC06M(  741),BNC07M(  741),BNC08M(  741),   &
3507 !      BNC09M(  741),BNC10M(  741),BNC11M(  741),BNC12M(  741),   &
3508 !      BNC13M(  741)
3509       REAL IN
3510 !
3511 ! *** Find position in arrays for binary activity coefficients
3512 !
3513       ipos = ibacpos(IN)
3514 !
3515 ! *** Assign values to return array
3516 !
3517       G01 = BNC01M(ipos)
3518       G02 = BNC02M(ipos)
3519       G03 = BNC03M(ipos)
3520       G04 = BNC04M(ipos)
3521       G05 = BNC05M(ipos)
3522       G06 = BNC06M(ipos)
3523       G07 = BNC07M(ipos)
3524       G08 = BNC08M(ipos)
3525       G09 = BNC09M(ipos)
3526       G10 = BNC10M(ipos)
3527       G11 = BNC11M(ipos)
3528       G12 = BNC12M(ipos)
3529 !
3530 ! *** Return point ; End of subroutine
3531 !
3532       RETURN
3533       END
3534
3535
3536 !=======================================================================
3537 !
3538 ! *** ISORROPIA CODE
3539 ! *** SUBROUTINE KM298
3540 ! *** CALCULATES BINARY ACTIVITY COEFFICIENTS BY KUSIK-MEISSNER METHOD.
3541 !     THE COMPUTATIONS HAVE BEEN PERFORMED AND THE RESULTS ARE STORED IN
3542 !     LOOKUP TABLES. THE IONIC ACTIVITY 'IN' IS INPUT, AND THE ARRAY
3543 !     'BINARR' IS RETURNED WITH THE BINARY COEFFICIENTS.
3544 !
3545 !     TEMPERATURE IS 298K
3546 !
3547 ! *** COPYRIGHT 1996-2006, UNIVERSITY OF MIAMI, CARNEGIE MELLON UNIVERSITY,
3548 ! *** GEORGIA INSTITUTE OF TECHNOLOGY
3549 ! *** WRITTEN BY ATHANASIOS NENES
3550 ! *** UPDATED BY CHRISTOS FOUNTOUKIS
3551 !
3552 !=======================================================================
3553 !
3554       SUBROUTINE KM298 (IN,G01,G02,G03,G04,G05,G06,G07,G08,G09,G10,   &
3555                            G11,G12)
3556       USE KMC298
3557 !
3558 ! *** Common block definition
3559 !
3560 !      COMMON /KMC298/   &
3561 !      BNC01M(  741),BNC02M(  741),BNC03M(  741),BNC04M(  741),   &
3562 !      BNC05M(  741),BNC06M(  741),BNC07M(  741),BNC08M(  741),   &
3563 !      BNC09M(  741),BNC10M(  741),BNC11M(  741),BNC12M(  741),   &
3564 !      BNC13M(  741)
3565       REAL IN
3566 !
3567 ! *** Find position in arrays for binary activity coefficients
3568 !
3569       ipos = ibacpos(IN)
3570 !
3571 ! *** Assign values to return array
3572 !
3573       G01 = BNC01M(ipos)
3574       G02 = BNC02M(ipos)
3575       G03 = BNC03M(ipos)
3576       G04 = BNC04M(ipos)
3577       G05 = BNC05M(ipos)
3578       G06 = BNC06M(ipos)
3579       G07 = BNC07M(ipos)
3580       G08 = BNC08M(ipos)
3581       G09 = BNC09M(ipos)
3582       G10 = BNC10M(ipos)
3583       G11 = BNC11M(ipos)
3584       G12 = BNC12M(ipos)
3585 !
3586 ! *** Return point ; End of subroutine
3587 !
3588       RETURN
3589       END
3590
3591
3592 !=======================================================================
3593 !
3594 ! *** ISORROPIA CODE
3595 ! *** SUBROUTINE KM323
3596 ! *** CALCULATES BINARY ACTIVITY COEFFICIENTS BY KUSIK-MEISSNER METHOD.
3597 !     THE COMPUTATIONS HAVE BEEN PERFORMED AND THE RESULTS ARE STORED IN
3598 !     LOOKUP TABLES. THE IONIC ACTIVITY 'IN' IS INPUT, AND THE ARRAY
3599 !     'BINARR' IS RETURNED WITH THE BINARY COEFFICIENTS.
3600 !
3601 !     TEMPERATURE IS 323K
3602 !
3603 ! *** COPYRIGHT 1996-2006, UNIVERSITY OF MIAMI, CARNEGIE MELLON UNIVERSITY,
3604 ! *** GEORGIA INSTITUTE OF TECHNOLOGY
3605 ! *** WRITTEN BY ATHANASIOS NENES
3606 ! *** UPDATED BY CHRISTOS FOUNTOUKIS
3607 !
3608 !=======================================================================
3609 !
3610       SUBROUTINE KM323 (IN,G01,G02,G03,G04,G05,G06,G07,G08,G09,G10,   &
3611                            G11,G12)
3612       USE KMC323
3613 !
3614 ! *** Common block definition
3615 !
3616 !      COMMON /KMC323/   &
3617 !      BNC01M(  741),BNC02M(  741),BNC03M(  741),BNC04M(  741),   &
3618 !      BNC05M(  741),BNC06M(  741),BNC07M(  741),BNC08M(  741),   &
3619 !      BNC09M(  741),BNC10M(  741),BNC11M(  741),BNC12M(  741),   &
3620 !      BNC13M(  741)
3621       REAL IN
3622 !
3623 ! *** Find position in arrays for binary activity coefficients
3624 !
3625       ipos = ibacpos(IN)
3626 !
3627 ! *** Assign values to return array
3628 !
3629       G01 = BNC01M(ipos)
3630       G02 = BNC02M(ipos)
3631       G03 = BNC03M(ipos)
3632       G04 = BNC04M(ipos)
3633       G05 = BNC05M(ipos)
3634       G06 = BNC06M(ipos)
3635       G07 = BNC07M(ipos)
3636       G08 = BNC08M(ipos)
3637       G09 = BNC09M(ipos)
3638       G10 = BNC10M(ipos)
3639       G11 = BNC11M(ipos)
3640       G12 = BNC12M(ipos)
3641 !
3642 ! *** Return point ; End of subroutine
3643 !
3644       RETURN
3645       END
3646
3647
3648
3649 !  *** TEMP = 198.0
3650
3651       BLOCK DATA KMCF198
3652        USE KMC198
3653 !
3654 !  *** Common block definition
3655 !
3656 !      COMMON /KMC198/   &
3657 !      BNC01M(  741),BNC02M(  741),BNC03M(  741),BNC04M(  741),   &
3658 !      BNC05M(  741),BNC06M(  741),BNC07M(  741),BNC08M(  741),   &
3659 !      BNC09M(  741),BNC10M(  741),BNC11M(  741),BNC12M(  741),   &
3660 !      BNC13M(  741)
3661 !
3662       END
3663
3664 !  *** TEMP = 223.0
3665
3666       BLOCK DATA KMCF223
3667       USE KMC223
3668 !
3669 !  *** Common block definition
3670 !
3671 !      COMMON /KMC223/   &
3672 !      BNC01M(  741),BNC02M(  741),BNC03M(  741),BNC04M(  741),   &
3673 !      BNC05M(  741),BNC06M(  741),BNC07M(  741),BNC08M(  741),   &
3674 !      BNC09M(  741),BNC10M(  741),BNC11M(  741),BNC12M(  741),   &
3675 !      BNC13M(  741)
3676 !
3677       END
3678
3679 !  *** TEMP = 248.0
3680
3681       BLOCK DATA KMCF248
3682        USE KMC248
3683 !
3684 !  *** Common block definition
3685 !
3686 !      COMMON /KMC248/   &
3687 !      BNC01M(  741),BNC02M(  741),BNC03M(  741),BNC04M(  741),   &
3688 !      BNC05M(  741),BNC06M(  741),BNC07M(  741),BNC08M(  741),   &
3689 !      BNC09M(  741),BNC10M(  741),BNC11M(  741),BNC12M(  741),   &
3690 !      BNC13M(  741)
3691       END
3692
3693 !  *** TEMP = 273.0
3694
3695       BLOCK DATA KMCF273
3696       USE KMC273
3697 !
3698 !  *** Common block definition
3699 !
3700 !      COMMON /KMC273/   &
3701 !      BNC01M(  741),BNC02M(  741),BNC03M(  741),BNC04M(  741),   &
3702 !      BNC05M(  741),BNC06M(  741),BNC07M(  741),BNC08M(  741),   &
3703 !      BNC09M(  741),BNC10M(  741),BNC11M(  741),BNC12M(  741),   &
3704 !      BNC13M(  741)
3705       END
3706
3707 !  *** TEMP = 298.0
3708
3709       BLOCK DATA KMCF298
3710       USE KMC298
3711 !
3712 !  *** Common block definition
3713 !
3714 !      COMMON /KMC298/   &
3715 !      BNC01M(  741),BNC02M(  741),BNC03M(  741),BNC04M(  741),   &
3716 !      BNC05M(  741),BNC06M(  741),BNC07M(  741),BNC08M(  741),   &
3717 !      BNC09M(  741),BNC10M(  741),BNC11M(  741),BNC12M(  741),   &
3718 !      BNC13M(  741)
3719       END
3720
3721 !  *** TEMP = 323.0
3722
3723       BLOCK DATA KMCF323
3724       USE KMC323
3725 !
3726 !  *** Common block definition
3727 !
3728 !      COMMON /KMC323/   &
3729 !      BNC01M(  741),BNC02M(  741),BNC03M(  741),BNC04M(  741),   &
3730 !      BNC05M(  741),BNC06M(  741),BNC07M(  741),BNC08M(  741),   &
3731 !      BNC09M(  741),BNC10M(  741),BNC11M(  741),BNC12M(  741),   &
3732 !      BNC13M(  741)
3733       END
3734
3735 !C*************************************************************************
3736 !C
3737 !C  TOOLBOX LIBRARY v.1.0 (May 1995)
3738 !C
3739 !C  Program unit   : SUBROUTINE CHRBLN
3740 !C  Purpose        : Position of last non-blank character in a string
3741 !C  Author         : Athanasios Nenes
3742 !C
3743 !C  ======================= ARGUMENTS / USAGE =============================
3744 !C
3745 !C  STR        is the CHARACTER variable containing the string examined
3746 !C  IBLK       is a INTEGER variable containing the position of last non
3747 !C             blank character. If string is all spaces (ie '   '), then
3748 !C             the value returned is 1.
3749 !C
3750 !C  EXAMPLE:
3751 !C             STR = 'TEST1.DAT     '
3752 !C             CALL CHRBLN (STR, IBLK)
3753 !C
3754 !C  after execution of this code segment, "IBLK" has the value "9", which
3755 !C  is the position of the last non-blank character of "STR".
3756 !C
3757 !C***********************************************************************
3758 !C
3759       SUBROUTINE CHRBLN (STR, IBLK)
3760 !C
3761 !C***********************************************************************
3762       CHARACTER*(*) STR
3763 !
3764       IBLK = 1                       ! Substring pointer (default=1)
3765       ILEN = LEN(STR)                ! Length of string
3766       DO 10 i=ILEN,1,-1
3767          IF (STR(i:i).NE.' ' .AND. STR(i:i).NE.CHAR(0)) THEN
3768             IBLK = i
3769             RETURN
3770          ENDIF
3771 10    CONTINUE
3772       RETURN
3773 !
3774       END
3775
3776
3777 !C*************************************************************************
3778 !C
3779 !C  TOOLBOX LIBRARY v.1.0 (May 1995)
3780 !C
3781 !C  Program unit   : SUBROUTINE SHFTRGHT
3782 !C  Purpose        : RIGHT-JUSTIFICATION FUNCTION ON A STRING
3783 !C  Author         : Athanasios Nenes
3784 !C
3785 !C  ======================= ARGUMENTS / USAGE =============================
3786 !C
3787 !C  STRING     is the CHARACTER variable with the string to be justified
3788 !C
3789 !C  EXAMPLE:
3790 !C             STRING    = 'AAAA    '
3791 !C             CALL SHFTRGHT (STRING)
3792 !C
3793 !C  after execution of this code segment, STRING contains the value
3794 !C  '    AAAA'.
3795 !C
3796 !C*************************************************************************
3797 !C
3798       SUBROUTINE SHFTRGHT (CHR)
3799 !C
3800 !C***********************************************************************
3801       CHARACTER CHR*(*)
3802 !
3803       I1  = LEN(CHR)             ! Total length of string
3804       CALL CHRBLN(CHR,I2)        ! Position of last non-blank character
3805       IF (I2.EQ.I1) RETURN
3806 !
3807       DO 10 I=I2,1,-1            ! Shift characters
3808          CHR(I1+I-I2:I1+I-I2) = CHR(I:I)
3809          CHR(I:I) = ' '
3810 10    CONTINUE
3811       RETURN
3812 !
3813       END
3814
3815
3816
3817
3818 !C*************************************************************************
3819 !C
3820 !C  TOOLBOX LIBRARY v.1.0 (May 1995)
3821 !C
3822 !C  Program unit   : SUBROUTINE RPLSTR
3823 !C  Purpose        : REPLACE CHARACTERS OCCURING IN A STRING
3824 !C  Author         : Athanasios Nenes
3825 !C
3826 !C  ======================= ARGUMENTS / USAGE =============================
3827 !C
3828 !C  STRING     is the CHARACTER variable with the string to be edited
3829 !C  OLD        is the old character which is to be replaced
3830 !C  NEW        is the new character which OLD is to be replaced with
3831 !C  IERR       is 0 if everything went well, is 1 if 'NEW' contains 'OLD'.
3832 !C             In this case, this is invalid, and no change is done.
3833 !C
3834 !C  EXAMPLE:
3835 !C             STRING    = 'AAAA'
3836 !C             OLD       = 'A'
3837 !C             NEW       = 'B'
3838 !C             CALL RPLSTR (STRING, OLD, NEW)
3839 !C
3840 !C  after execution of this code segment, STRING contains the value
3841 !C  'BBBB'.
3842 !C
3843 !C*************************************************************************
3844 !C
3845       SUBROUTINE RPLSTR (STRING, OLD, NEW, IERR)
3846 !C
3847 !C***********************************************************************
3848       CHARACTER STRING*(*), OLD*(*), NEW*(*)
3849 !
3850 ! *** INITIALIZE ********************************************************
3851 !
3852       ILO = LEN(OLD)
3853 !
3854 ! *** CHECK AND SEE IF 'NEW' CONTAINS 'OLD', WHICH CANNOT ***************
3855 !
3856       IP = INDEX(NEW,OLD)
3857       IF (IP.NE.0) THEN
3858          IERR = 1
3859          RETURN
3860       ELSE
3861          IERR = 0
3862       ENDIF
3863 !
3864 ! *** PROCEED WITH REPLACING *******************************************
3865 !
3866 10    IP = INDEX(STRING,OLD)      ! SEE IF 'OLD' EXISTS IN 'STRING'
3867       IF (IP.EQ.0) RETURN         ! 'OLD' DOES NOT EXIST ; RETURN
3868       STRING(IP:IP+ILO-1) = NEW   ! REPLACE SUBSTRING 'OLD' WITH 'NEW'
3869       GOTO 10                     ! GO FOR NEW OCCURANCE OF 'OLD'
3870 !
3871       END
3872
3873
3874 !C*************************************************************************
3875 !C
3876 !C  TOOLBOX LIBRARY v.1.0 (May 1995)
3877 !C
3878 !C  Program unit   : SUBROUTINE INPTD
3879 !C  Purpose        : Prompts user for a value (DOUBLE). A default value
3880 !C                   is provided, so if user presses <Enter>, the default
3881 !C                   is used.
3882 !C  Author         : Athanasios Nenes
3883 !C
3884 !C  ======================= ARGUMENTS / USAGE =============================
3885 !C
3886 !C  VAR        is the DOUBLE PRECISION variable which value is to be saved
3887 !C  DEF        is a DOUBLE PRECISION variable, with the default value of VAR.
3888 !C  PROMPT     is a CHARACTER varible containing the prompt string.
3889 !C  PRFMT      is a CHARACTER variable containing the FORMAT specifier
3890 !C             for the default value DEF.
3891 !C  IERR       is an INTEGER error flag, and has the values:
3892 !C             0 - No error detected.
3893 !C             1 - Invalid FORMAT and/or Invalid default value.
3894 !C             2 - Bad value specified by user
3895 !C
3896 !C  EXAMPLE:
3897 !C             CALL INPTD (VAR, 1.0D0, 'Give value for A ', '*', Ierr)
3898 !C
3899 !C  after execution of this code segment, the user is prompted for the
3900 !C  value of variable VAR. If <Enter> is pressed (ie no value is specified)
3901 !C  then 1.0 is assigned to VAR. The default value is displayed in free-
3902 !C  format. The error status is specified by variable Ierr
3903 !C
3904 !C***********************************************************************
3905 !C
3906       SUBROUTINE INPTD (VAR, DEF, PROMPT, PRFMT, IERR)
3907 !C
3908 !C***********************************************************************
3909       CHARACTER PROMPT*(*), PRFMT*(*), BUFFER*128
3910       DOUBLE PRECISION DEF, VAR
3911       INTEGER IERR
3912 !
3913       IERR = 0
3914 !
3915 ! *** WRITE DEFAULT VALUE TO WORK BUFFER *******************************
3916 !
3917       WRITE (BUFFER, FMT=PRFMT, ERR=10) DEF
3918       CALL CHRBLN (BUFFER, IEND)
3919 !
3920 ! *** PROMPT USER FOR INPUT AND READ IT ********************************
3921 !
3922       WRITE (*,*) PROMPT,' [',BUFFER(1:IEND),']: '
3923       READ  (*, '(A)', ERR=20, END=20) BUFFER
3924       CALL CHRBLN (BUFFER,IEND)
3925 !
3926 ! *** READ DATA OR SET DEFAULT ? ****************************************
3927 !
3928       IF (IEND.EQ.1 .AND. BUFFER(1:1).EQ.' ') THEN
3929          VAR = DEF
3930       ELSE
3931          READ (BUFFER, *, ERR=20, END=20) VAR
3932       ENDIF
3933 !
3934 ! *** RETURN POINT ******************************************************
3935 !
3936 30    RETURN
3937 !
3938 ! *** ERROR HANDLER *****************************************************
3939 !
3940 10    IERR = 1       ! Bad FORMAT and/or bad default value
3941       GOTO 30
3942 !
3943 20    IERR = 2       ! Bad number given by user
3944       GOTO 30
3945 !
3946       END
3947
3948
3949 !C*************************************************************************
3950 !C
3951 !C  TOOLBOX LIBRARY v.1.0 (May 1995)
3952 !C
3953 !C  Program unit   : SUBROUTINE Pushend
3954 !C  Purpose        : Positions the pointer of a sequential file at its end
3955 !C                   Simulates the ACCESS='APPEND' clause of a F77L OPEN
3956 !C                   statement with Standard Fortran commands.
3957 !C
3958 !C  ======================= ARGUMENTS / USAGE =============================
3959 !C
3960 !C  Iunit      is a INTEGER variable, the file unit which the file is
3961 !C             connected to.
3962 !C
3963 !C  EXAMPLE:
3964 !C             CALL PUSHEND (10)
3965 !C
3966 !C  after execution of this code segment, the pointer of unit 10 is
3967 !C  pushed to its end.
3968 !C
3969 !C***********************************************************************
3970 !C
3971       SUBROUTINE Pushend (Iunit)
3972 !C
3973 !C***********************************************************************
3974 !
3975       LOGICAL OPNED
3976 !
3977 ! *** INQUIRE IF Iunit CONNECTED TO FILE ********************************
3978 !
3979       INQUIRE (UNIT=Iunit, OPENED=OPNED)
3980       IF (.NOT.OPNED) GOTO 25
3981 !
3982 ! *** Iunit CONNECTED, PUSH POINTER TO END ******************************
3983 !
3984 10    READ (Iunit,'()', ERR=20, END=20)
3985       GOTO 10
3986 !
3987 ! *** RETURN POINT ******************************************************
3988 !
3989 20    BACKSPACE (Iunit)
3990 25    RETURN
3991       END
3992
3993
3994
3995 !C*************************************************************************
3996 !C
3997 !C  TOOLBOX LIBRARY v.1.0 (May 1995)
3998 !C
3999 !C  Program unit   : SUBROUTINE APPENDEXT
4000 !C  Purpose        : Fix extension in file name string
4001 !C
4002 !C  ======================= ARGUMENTS / USAGE =============================
4003 !C
4004 !C  Filename   is the CHARACTER variable with the file name
4005 !C  Defext     is the CHARACTER variable with extension (including '.',
4006 !C             ex. '.DAT')
4007 !C  Overwrite  is a LOGICAL value, .TRUE. overwrites any existing extension
4008 !C             in "Filename" with "Defext", .FALSE. puts "Defext" only if
4009 !C             there is no extension in "Filename".
4010 !C
4011 !C  EXAMPLE:
4012 !C             FILENAME1 = 'TEST.DAT'
4013 !C             FILENAME2 = 'TEST.DAT'
4014 !C             CALL APPENDEXT (FILENAME1, '.TXT', .FALSE.)
4015 !C             CALL APPENDEXT (FILENAME2, '.TXT', .TRUE. )
4016 !C
4017 !C  after execution of this code segment, "FILENAME1" has the value
4018 !C  'TEST.DAT', while "FILENAME2" has the value 'TEST.TXT'
4019 !C
4020 !C***********************************************************************
4021 !C
4022       SUBROUTINE Appendext (Filename, Defext, Overwrite)
4023 !C
4024 !C***********************************************************************
4025       CHARACTER*(*) Filename, Defext
4026       LOGICAL       Overwrite
4027 !
4028       CALL CHRBLN (Filename, Iend)
4029       IF (Filename(1:1).EQ.' ' .AND. Iend.EQ.1) RETURN  ! Filename empty
4030       Idot = INDEX (Filename, '.')                      ! Append extension ?
4031       IF (Idot.EQ.0) Filename = Filename(1:Iend)//Defext
4032       IF (Overwrite .AND. Idot.NE.0)   &
4033                     Filename = Filename(:Idot-1)//Defext
4034       RETURN
4035       END
4036
4037
4038
4039
4040
4041 !=======================================================================
4042 !
4043 ! *** ISORROPIA CODE
4044 ! *** SUBROUTINE POLY3
4045 ! *** FINDS THE REAL ROOTS OF THE THIRD ORDER ALGEBRAIC EQUATION:
4046 !     X**3 + A1*X**2 + A2*X + A3 = 0.0
4047 !     THE EQUATION IS SOLVED ANALYTICALLY.
4048 !
4049 !     PARAMETERS A1, A2, A3 ARE SPECIFIED BY THE USER. THE MINIMUM
4050 !     NONEGATIVE ROOT IS RETURNED IN VARIABLE 'ROOT'. IF NO ROOT IS
4051 !     FOUND (WHICH IS GREATER THAN ZERO), ROOT HAS THE VALUE 1D30.
4052 !     AND THE FLAG ISLV HAS A VALUE GREATER THAN ZERO.
4053 !
4054 !     SOLUTION FORMULA IS FOUND IN PAGE 32 OF:
4055 !     MATHEMATICAL HANDBOOK OF FORMULAS AND TABLES
4056 !     SCHAUM'S OUTLINE SERIES
4057 !     MURRAY SPIEGER, McGRAW-HILL, NEW YORK, 1968
4058 !     (GREEK TRANSLATION: BY SOTIRIOS PERSIDES, ESPI, ATHENS, 1976)
4059 !
4060 !     A SPECIAL CASE IS CONSIDERED SEPERATELY ; WHEN A3 = 0, THEN
4061 !     ONE ROOT IS X=0.0, AND THE OTHER TWO FROM THE SOLUTION OF THE
4062 !     QUADRATIC EQUATION X**2 + A1*X + A2 = 0.0
4063 !     THIS SPECIAL CASE IS CONSIDERED BECAUSE THE ANALYTICAL FORMULA
4064 !     DOES NOT YIELD ACCURATE RESULTS (DUE TO NUMERICAL ROUNDOFF ERRORS)
4065 !
4066 ! *** COPYRIGHT 1996-2006, UNIVERSITY OF MIAMI, CARNEGIE MELLON UNIVERSITY,
4067 ! *** GEORGIA INSTITUTE OF TECHNOLOGY
4068 ! *** WRITTEN BY ATHANASIOS NENES
4069 ! *** UPDATED BY CHRISTOS FOUNTOUKIS
4070 !
4071 !=======================================================================
4072 !
4073       SUBROUTINE POLY3 (A1, A2, A3, ROOT, ISLV)
4074 !
4075       IMPLICIT DOUBLE PRECISION (A-H, O-Z)
4076       PARAMETER (EXPON=1.D0/3.D0,     ZERO=0.D0, THET1=120.D0/180.D0,    &
4077                  THET2=240.D0/180.D0, PI=3.14159265358932, EPS=1D-50)
4078       DOUBLE PRECISION  X(3)
4079 !
4080 ! *** SPECIAL CASE : QUADRATIC*X EQUATION *****************************
4081 !
4082       IF (ABS(A3).LE.EPS) THEN
4083          ISLV = 1
4084          IX   = 1
4085          X(1) = ZERO
4086          D    = A1*A1-4.D0*A2
4087          IF (D.GE.ZERO) THEN
4088             IX   = 3
4089             SQD  = SQRT(D)
4090             X(2) = 0.5*(-A1+SQD)
4091             X(3) = 0.5*(-A1-SQD)
4092          ENDIF
4093       ELSE
4094 !
4095 ! *** NORMAL CASE : CUBIC EQUATION ************************************
4096 !
4097 ! DEFINE PARAMETERS Q, R, S, T, D
4098 !
4099          ISLV= 1
4100          Q   = (3.D0*A2 - A1*A1)/9.D0
4101          R   = (9.D0*A1*A2 - 27.D0*A3 - 2.D0*A1*A1*A1)/54.D0
4102          D   = Q*Q*Q + R*R
4103 !
4104 ! *** CALCULATE ROOTS *************************************************
4105 !
4106 !  D < 0, THREE REAL ROOTS
4107 !
4108          IF (D.LT.-EPS) THEN        ! D < -EPS  : D < ZERO
4109             IX   = 3
4110             THET = EXPON*ACOS(R/SQRT(-Q*Q*Q))
4111             COEF = 2.D0*SQRT(-Q)
4112             X(1) = COEF*COS(THET)            - EXPON*A1
4113             X(2) = COEF*COS(THET + THET1*PI) - EXPON*A1
4114             X(3) = COEF*COS(THET + THET2*PI) - EXPON*A1
4115 !
4116 !  D = 0, THREE REAL (ONE DOUBLE) ROOTS
4117 !
4118          ELSE IF (D.LE.EPS) THEN    ! -EPS <= D <= EPS  : D = ZERO
4119             IX   = 2
4120             SSIG = SIGN (1.D0, R)
4121             S    = SSIG*(ABS(R))**EXPON
4122             X(1) = 2.D0*S  - EXPON*A1
4123             X(2) =     -S  - EXPON*A1
4124 !
4125 !  D > 0, ONE REAL ROOT
4126 !
4127          ELSE                       ! D > EPS  : D > ZERO
4128             IX   = 1
4129             SQD  = SQRT(D)
4130             SSIG = SIGN (1.D0, R+SQD)       ! TRANSFER SIGN TO SSIG
4131             TSIG = SIGN (1.D0, R-SQD)
4132             S    = SSIG*(ABS(R+SQD))**EXPON ! EXPONENTIATE ABS()
4133             T    = TSIG*(ABS(R-SQD))**EXPON
4134             X(1) = S + T - EXPON*A1
4135          ENDIF
4136       ENDIF
4137 !
4138 ! *** SELECT APPROPRIATE ROOT *****************************************
4139 !
4140       ROOT = 1.D30
4141       DO 10 I=1,IX
4142          IF (X(I).GT.ZERO) THEN
4143             ROOT = MIN (ROOT, X(I))
4144             ISLV = 0
4145          ENDIF
4146 10    CONTINUE
4147 !
4148 ! *** END OF SUBROUTINE POLY3 *****************************************
4149 !
4150       RETURN
4151       END
4152
4153
4154
4155
4156 !=======================================================================
4157 !
4158 ! *** ISORROPIA CODE
4159 ! *** SUBROUTINE POLY3B
4160 ! *** FINDS A REAL ROOT OF THE THIRD ORDER ALGEBRAIC EQUATION:
4161 !     X**3 + A1*X**2 + A2*X + A3 = 0.0
4162 !     THE EQUATION IS SOLVED NUMERICALLY (BISECTION).
4163 !
4164 !     PARAMETERS A1, A2, A3 ARE SPECIFIED BY THE USER. THE MINIMUM
4165 !     NONEGATIVE ROOT IS RETURNED IN VARIABLE 'ROOT'. IF NO ROOT IS
4166 !     FOUND (WHICH IS GREATER THAN ZERO), ROOT HAS THE VALUE 1D30.
4167 !     AND THE FLAG ISLV HAS A VALUE GREATER THAN ZERO.
4168 !
4169 !     RTLW, RTHI DEFINE THE INTERVAL WHICH THE ROOT IS LOOKED FOR.
4170 !
4171 !=======================================================================
4172 !
4173       SUBROUTINE POLY3B (A1, A2, A3, RTLW, RTHI, ROOT, ISLV)
4174 !
4175       IMPLICIT DOUBLE PRECISION (A-H, O-Z)
4176       PARAMETER (ZERO=0.D0, EPS=1D-15, MAXIT=100, NDIV=5)
4177 !
4178       FUNC(X) = X**3.d0 + A1*X**2.0 + A2*X + A3
4179 !
4180 ! *** INITIAL VALUES FOR BISECTION *************************************
4181 !
4182       X1   = RTLW
4183       Y1   = FUNC(X1)
4184       IF (ABS(Y1).LE.EPS) THEN     ! Is low a root?
4185          ROOT = RTLW
4186          GOTO 50
4187       ENDIF
4188 !
4189 ! *** ROOT TRACKING ; FOR THE RANGE OF HI AND LO ***********************
4190 !
4191       DX = (RTHI-RTLW)/REAL(NDIV)
4192       DO 10 I=1,NDIV
4193          X2 = X1+DX
4194          Y2 = FUNC (X2)
4195          IF (SIGN(1.d0,Y1)*SIGN(1.d0,Y2) .LT. ZERO) GOTO 20 ! (Y1*Y2.LT.ZERO)
4196          X1 = X2
4197          Y1 = Y2
4198 10    CONTINUE
4199 !
4200 ! *** NO SUBDIVISION WITH SOLUTION FOUND
4201 !
4202       IF (ABS(Y2) .LT. EPS) THEN   ! X2 is a root
4203          ROOT = X2
4204       ELSE
4205          ROOT = 1.d30
4206          ISLV = 1
4207       ENDIF
4208       GOTO 50
4209 !
4210 ! *** BISECTION *******************************************************
4211 !
4212 20    DO 30 I=1,MAXIT
4213          X3 = 0.5*(X1+X2)
4214          Y3 = FUNC (X3)
4215          IF (SIGN(1.d0,Y1)*SIGN(1.d0,Y3) .LE. ZERO) THEN  ! (Y1*Y3 .LE. ZERO)
4216             Y2    = Y3
4217             X2    = X3
4218          ELSE
4219             Y1    = Y3
4220             X1    = X3
4221          ENDIF
4222          IF (ABS(X2-X1) .LE. EPS*X1) GOTO 40
4223 30    CONTINUE
4224 !
4225 ! *** CONVERGED ; RETURN ***********************************************
4226 !
4227 40    X3   = 0.5*(X1+X2)
4228       Y3   = FUNC (X3)
4229       ROOT = X3
4230       ISLV = 0
4231 !
4232 50    RETURN
4233 !
4234 ! *** END OF SUBROUTINE POLY3B *****************************************
4235 !
4236       END
4237
4238
4239
4240 !cc      PROGRAM DRIVER
4241 !cc      DOUBLE PRECISION ROOT
4242 !ccC
4243 !cc      CALL POLY3 (-1.d0, 1.d0, -1.d0, ROOT, ISLV)
4244 !cc      IF (ISLV.NE.0) STOP 'Error in POLY3'
4245 !cc      WRITE (*,*) 'Root=', ROOT
4246 !ccC
4247 !cc      CALL POLY3B (-1.d0, 1.d0, -1.d0, -10.d0, 10.d0, ROOT, ISLV)
4248 !cc      IF (ISLV.NE.0) STOP 'Error in POLY3B'
4249 !cc      WRITE (*,*) 'Root=', ROOT
4250 !ccC
4251 !cc      END
4252 !=======================================================================
4253 !
4254 ! *** ISORROPIA CODE
4255 ! *** FUNCTION EX10
4256 ! *** 10^X FUNCTION ; ALTERNATE OF LIBRARY ROUTINE ; USED BECAUSE IT IS
4257 !     MUCH FASTER BUT WITHOUT GREAT LOSS IN ACCURACY. ,
4258 !     MAXIMUM ERROR IS 2%, EXECUTION TIME IS 42% OF THE LIBRARY ROUTINE
4259 !     (ON A 80286/80287 MACHINE, using Lahey FORTRAN 77 v.3.0).
4260 !
4261 !     EXPONENT RANGE IS BETWEEN -K AND K (K IS THE REAL ARGUMENT 'K')
4262 !     MAX VALUE FOR K: 9.999
4263 !     IF X < -K, X IS SET TO -K, IF X > K, X IS SET TO K
4264 !
4265 !     THE EXPONENT IS CALCULATED BY THE PRODUCT ADEC*AINT, WHERE ADEC
4266 !     IS THE MANTISSA AND AINT IS THE MAGNITUDE (EXPONENT). BOTH
4267 !     MANTISSA AND MAGNITUDE ARE PRE-CALCULATED AND STORED IN LOOKUP
4268 !     TABLES ; THIS LEADS TO THE INCREASED SPEED.
4269 !
4270 ! *** COPYRIGHT 1996-2006, UNIVERSITY OF MIAMI, CARNEGIE MELLON UNIVERSITY,
4271 ! *** GEORGIA INSTITUTE OF TECHNOLOGY
4272 ! *** WRITTEN BY ATHANASIOS NENES
4273 ! *** UPDATED BY CHRISTOS FOUNTOUKIS
4274 !
4275 !=======================================================================
4276 !
4277       FUNCTION EX10(X,K)
4278       USE EXPNC
4279       REAL    X, EX10, Y , K
4280       INTEGER K1, K2
4281 !      COMMON /EXPNC/ AINT10(20), ADEC10(200)
4282 !
4283 ! *** LIMIT X TO [-K, K] RANGE *****************************************
4284 !
4285       Y    = MAX(-K, MIN(X,K))   ! MIN: -9.999, MAX: 9.999
4286 !
4287 ! *** GET INTEGER AND DECIMAL PART *************************************
4288 !
4289       K1   = INT(Y)
4290       K2   = INT(100*(Y-K1))
4291 !
4292 ! *** CALCULATE EXP FUNCTION *******************************************
4293 !
4294       EX10 = AINT10(K1+10)*ADEC10(K2+100)
4295 !
4296 ! *** END OF EXP FUNCTION **********************************************
4297 !
4298       RETURN
4299       END
4300
4301
4302 !=======================================================================
4303 !
4304 ! *** ISORROPIA CODE
4305 ! *** BLOCK DATA EXPON
4306 ! *** CONTAINS DATA FOR EXPONENT ARRAYS NEEDED IN FUNCTION EXP10
4307 !
4308 ! *** COPYRIGHT 1996-2006, UNIVERSITY OF MIAMI, CARNEGIE MELLON UNIVERSITY,
4309 ! *** GEORGIA INSTITUTE OF TECHNOLOGY
4310 ! *** WRITTEN BY ATHANASIOS NENES
4311 ! *** UPDATED BY CHRISTOS FOUNTOUKIS
4312 !
4313 !=======================================================================
4314 !
4315       BLOCK DATA EXPON
4316 !
4317 ! *** Common block definition
4318 !
4319       USE EXPNC
4320 !      REAL AINT10, ADEC10
4321 !      COMMON /EXPNC/ AINT10(20), ADEC10(200)
4322 !
4323 ! *** Integer part
4324 !
4325 !      DATA AINT10/   &
4326 !       0.1000E-08, 0.1000E-07, 0.1000E-06, 0.1000E-05, 0.1000E-04,   &
4327 !       0.1000E-03, 0.1000E-02, 0.1000E-01, 0.1000E+00, 0.1000E+01,   &
4328 !       0.1000E+02, 0.1000E+03, 0.1000E+04, 0.1000E+05, 0.1000E+06,   &
4329 !       0.1000E+07, 0.1000E+08, 0.1000E+09, 0.1000E+10, 0.1000E+11   &
4330 !       /
4331 !
4332 ! *** decimal part
4333 !
4334 !      DATA (ADEC10(I),I=1,100)/   &
4335 !       0.1023E+00, 0.1047E+00, 0.1072E+00, 0.1096E+00, 0.1122E+00,   &
4336 !       0.1148E+00, 0.1175E+00, 0.1202E+00, 0.1230E+00, 0.1259E+00,   &
4337 !       0.1288E+00, 0.1318E+00, 0.1349E+00, 0.1380E+00, 0.1413E+00,   &
4338 !       0.1445E+00, 0.1479E+00, 0.1514E+00, 0.1549E+00, 0.1585E+00,   &
4339 !       0.1622E+00, 0.1660E+00, 0.1698E+00, 0.1738E+00, 0.1778E+00,   &
4340 !       0.1820E+00, 0.1862E+00, 0.1905E+00, 0.1950E+00, 0.1995E+00,   &
4341 !       0.2042E+00, 0.2089E+00, 0.2138E+00, 0.2188E+00, 0.2239E+00,   &
4342 !       0.2291E+00, 0.2344E+00, 0.2399E+00, 0.2455E+00, 0.2512E+00,   &
4343 !       0.2570E+00, 0.2630E+00, 0.2692E+00, 0.2754E+00, 0.2818E+00,   &
4344 !       0.2884E+00, 0.2951E+00, 0.3020E+00, 0.3090E+00, 0.3162E+00,   &
4345 !       0.3236E+00, 0.3311E+00, 0.3388E+00, 0.3467E+00, 0.3548E+00,   &
4346 !       0.3631E+00, 0.3715E+00, 0.3802E+00, 0.3890E+00, 0.3981E+00,   &
4347 !       0.4074E+00, 0.4169E+00, 0.4266E+00, 0.4365E+00, 0.4467E+00,   &
4348 !      0.4571E+00, 0.4677E+00, 0.4786E+00, 0.4898E+00, 0.5012E+00,   &
4349 !       0.5129E+00, 0.5248E+00, 0.5370E+00, 0.5495E+00, 0.5623E+00,   &
4350 !       0.5754E+00, 0.5888E+00, 0.6026E+00, 0.6166E+00, 0.6310E+00,   &
4351 !       0.6457E+00, 0.6607E+00, 0.6761E+00, 0.6918E+00, 0.7079E+00,   &
4352 !       0.7244E+00, 0.7413E+00, 0.7586E+00, 0.7762E+00, 0.7943E+00,   &
4353 !       0.8128E+00, 0.8318E+00, 0.8511E+00, 0.8710E+00, 0.8913E+00,   &
4354 !       0.9120E+00, 0.9333E+00, 0.9550E+00, 0.9772E+00, 0.1000E+01/
4355 !
4356 !      DATA (ADEC10(I),I=101,200)/   &
4357 !       0.1023E+01, 0.1047E+01, 0.1072E+01, 0.1096E+01, 0.1122E+01,   &
4358 !       0.1148E+01, 0.1175E+01, 0.1202E+01, 0.1230E+01, 0.1259E+01,   &
4359 !       0.1288E+01, 0.1318E+01, 0.1349E+01, 0.1380E+01, 0.1413E+01,   &
4360 !       0.1445E+01, 0.1479E+01, 0.1514E+01, 0.1549E+01, 0.1585E+01,   &
4361 !       0.1622E+01, 0.1660E+01, 0.1698E+01, 0.1738E+01, 0.1778E+01,   &
4362 !       0.1820E+01, 0.1862E+01, 0.1905E+01, 0.1950E+01, 0.1995E+01,   &
4363 !       0.2042E+01, 0.2089E+01, 0.2138E+01, 0.2188E+01, 0.2239E+01,   &
4364 !       0.2291E+01, 0.2344E+01, 0.2399E+01, 0.2455E+01, 0.2512E+01,   &
4365 !       0.2570E+01, 0.2630E+01, 0.2692E+01, 0.2754E+01, 0.2818E+01,   &
4366 !       0.2884E+01, 0.2951E+01, 0.3020E+01, 0.3090E+01, 0.3162E+01,   &
4367 !       0.3236E+01, 0.3311E+01, 0.3388E+01, 0.3467E+01, 0.3548E+01,   &
4368 !       0.3631E+01, 0.3715E+01, 0.3802E+01, 0.3890E+01, 0.3981E+01,   &
4369 !       0.4074E+01, 0.4169E+01, 0.4266E+01, 0.4365E+01, 0.4467E+01,   &
4370 !       0.4571E+01, 0.4677E+01, 0.4786E+01, 0.4898E+01, 0.5012E+01,   &
4371 !       0.5129E+01, 0.5248E+01, 0.5370E+01, 0.5495E+01, 0.5623E+01,   &
4372 !       0.5754E+01, 0.5888E+01, 0.6026E+01, 0.6166E+01, 0.6310E+01,   &
4373 !       0.6457E+01, 0.6607E+01, 0.6761E+01, 0.6918E+01, 0.7079E+01,   &
4374 !       0.7244E+01, 0.7413E+01, 0.7586E+01, 0.7762E+01, 0.7943E+01,   &
4375 !       0.8128E+01, 0.8318E+01, 0.8511E+01, 0.8710E+01, 0.8913E+01,   &
4376 !       0.9120E+01, 0.9333E+01, 0.9550E+01, 0.9772E+01, 0.1000E+02   &
4377 !       /
4378 !
4379 ! *** END OF BLOCK DATA EXPON ******************************************
4380 !
4381       END
4382 !=======================================================================
4383 !
4384 ! *** ISORROPIA CODE
4385 ! *** SUBROUTINE ISOPLUS
4386 ! *** THIS SUBROUTINE IS THE MASTER ROUTINE FOR THE ISORROPIA-PLUS
4387 !     THERMODYNAMIC EQUILIBRIUM AEROSOL MODEL (VERSION 1.0)
4388 !
4389 ! *** NOTE: THIS SUBROUTINE IS INCLUDED FOR BACKWARD COMPATABILITY ONLY.
4390 !     A PROGRAMMER SHOULD USE THE MORE COMPLETE SUBROUTINE ISOROPIA INSTEAD
4391 !
4392 ! ======================== ARGUMENTS / USAGE ===========================
4393 !
4394 !  INPUT:
4395 !  1. [WI]
4396 !     DOUBLE PRECISION array of length [5].
4397 !     Concentrations, expressed in moles/m3. Depending on the type of
4398 !     problem solved, WI contains either GAS+AEROSOL or AEROSOL only
4399 !     concentratios.
4400 !     WI(1) - sodium
4401 !     WI(2) - sulfate
4402 !     WI(3) - ammonium
4403 !     WI(4) - nitrate
4404 !     WI(5) - chloride
4405 !
4406 !  2. [RHI]
4407 !     DOUBLE PRECISION variable.
4408 !     Ambient relative humidity expressed in a (0,1) scale.
4409 !
4410 !  3. [TEMPI]
4411 !     DOUBLE PRECISION variable.
4412 !     Ambient temperature expressed in Kelvins.
4413 !
4414 !  4. [IPROB]
4415 !     INTEGER variable.
4416 !     The type of problem solved.
4417 !     IPROB = 0  - Forward problem is solved. In this case, array WI
4418 !                  contains GAS and AEROSOL concentrations together.
4419 !     IPROB = 1  - Reverse problem is solved. In this case, array WI
4420 !                  contains AEROSOL concentrations only.
4421 !
4422 !  OUTPUT:
4423 !  1. [GAS]
4424 !     DOUBLE PRECISION array of length [03].
4425 !     Gaseous species concentrations, expressed in moles/m3.
4426 !     GAS(1) - NH3
4427 !     GAS(2) - HNO3
4428 !     GAS(3) - HCl
4429 !
4430 !  2. [AERLIQ]
4431 !     DOUBLE PRECISION array of length [11].
4432 !     Liquid aerosol species concentrations, expressed in moles/m3.
4433 !     AERLIQ(01) - H+(aq)
4434 !     AERLIQ(02) - Na+(aq)
4435 !     AERLIQ(03) - NH4+(aq)
4436 !     AERLIQ(04) - Cl-(aq)
4437 !     AERLIQ(05) - SO4--(aq)
4438 !     AERLIQ(06) - HSO4-(aq)
4439 !     AERLIQ(07) - NO3-(aq)
4440 !     AERLIQ(08) - H2O
4441 !     AERLIQ(09) - NH3(aq) (undissociated)
4442 !     AERLIQ(10) - HNCl(aq) (undissociated)
4443 !     AERLIQ(11) - HNO3(aq) (undissociated)
4444 !
4445 !  3. [AERSLD]
4446 !     DOUBLE PRECISION array of length [09].
4447 !     Solid aerosol species concentrations, expressed in moles/m3.
4448 !     AERSLD(01) - NaNO3(s)
4449 !     AERSLD(02) - NH4NO3(s)
4450 !     AERSLD(03) - NaCl(s)
4451 !     AERSLD(04) - NH4Cl(s)
4452 !     AERSLD(05) - Na2SO4(s)
4453 !     AERSLD(06) - (NH4)2SO4(s)
4454 !     AERSLD(07) - NaHSO4(s)
4455 !     AERSLD(08) - NH4HSO4(s)
4456 !     AERSLD(09) - (NH4)4H(SO4)2(s)
4457 !
4458 !  4. [DRY]
4459 !     LOGICAL variable.
4460 !     Contains information about the physical state of the system.
4461 !     .TRUE. - There is no aqueous phase present
4462 !     .FALSE.- There is an aqueous phase present
4463 !
4464 ! *** COPYRIGHT 1996-2006, UNIVERSITY OF MIAMI, CARNEGIE MELLON UNIVERSITY,
4465 ! *** GEORGIA INSTITUTE OF TECHNOLOGY
4466 ! *** WRITTEN BY ATHANASIOS NENES
4467 ! *** UPDATED BY CHRISTOS FOUNTOUKIS
4468 !
4469 !=======================================================================
4470 !
4471       SUBROUTINE ISOPLUS (WI,  RHI,    TEMPI,  IPROBI,    &
4472                           GAS, AERLIQ, AERSLD, DRYI   )
4473       USE ISRPIA
4474       IMPLICIT DOUBLE PRECISION (A-H, O-Z)
4475 !      INCLUDE 'ISRPIA.INC'
4476       DIMENSION WI(NCOMP), GAS(NGASAQ), AERLIQ(NIONS+NGASAQ+1),   &
4477                 AERSLD(NSLDS)
4478       LOGICAL   DRYI
4479 !
4480 ! *** PROBLEM TYPE (0=FOREWARD, 1=REVERSE) ******************************
4481 !
4482       IPROB = IPROBI
4483 !
4484 ! *** SOLVE FOREWARD PROBLEM ********************************************
4485 !
4486       IF (IPROB.EQ.0) THEN
4487          IF (WI(1)+WI(2)+WI(3)+WI(4)+WI(5) .LE. TINY) THEN ! Everything=0
4488             CALL INIT1 (WI, RHI, TEMPI)
4489          ELSE IF (WI(1)+WI(4)+WI(5) .LE. TINY) THEN        ! Na,Cl,NO3=0
4490             CALL ISRP1F (WI, RHI, TEMPI)
4491          ELSE IF (WI(1)+WI(5) .LE. TINY) THEN              ! Na,Cl=0
4492             CALL ISRP2F (WI, RHI, TEMPI)
4493          ELSE
4494             CALL ISRP3F (WI, RHI, TEMPI)
4495          ENDIF
4496 !
4497 ! *** SOLVE REVERSE PROBLEM *********************************************
4498 !
4499       ELSE
4500          IF (WI(1)+WI(2)+WI(3)+WI(4)+WI(5) .LE. TINY) THEN ! Everything=0
4501             CALL INIT1 (WI, RHI, TEMPI)
4502          ELSE IF (WI(1)+WI(4)+WI(5) .LE. TINY) THEN        ! Na,Cl,NO3=0
4503             CALL ISRP1R (WI, RHI, TEMPI)
4504          ELSE IF (WI(1)+WI(5) .LE. TINY) THEN              ! Na,Cl=0
4505             CALL ISRP2R (WI, RHI, TEMPI)
4506          ELSE
4507             CALL ISRP3R (WI, RHI, TEMPI)
4508          ENDIF
4509       ENDIF
4510 !
4511 ! *** SAVE RESULTS TO ARRAYS (units = mole/m3, kg/m3 for water) *********
4512 !
4513       GAS(1) = GNH3
4514       GAS(2) = GHNO3
4515       GAS(3) = GHCL
4516 !
4517       DO 10 I=1,NIONS
4518          AERLIQ(I) = MOLAL(I)
4519   10  CONTINUE
4520       AERLIQ(NIONS+1) = WATER*1.0D3/18.0D0
4521       DO 20 I=1,NGASAQ
4522          AERLIQ(NIONS+1+I) = GASAQ(I)
4523   20  CONTINUE
4524 !
4525       AERSLD(1) = CNANO3
4526       AERSLD(2) = CNH4NO3
4527       AERSLD(3) = CNACL
4528       AERSLD(4) = CNH4CL
4529       AERSLD(5) = CNA2SO4
4530       AERSLD(6) = CNH42S4
4531       AERSLD(7) = CNAHSO4
4532       AERSLD(8) = CNH4HS4
4533       AERSLD(9) = CLC
4534 !
4535       DRYI = WATER.LE.TINY
4536 !
4537       RETURN
4538 !
4539 ! *** END OF SUBROUTINE ISOPLUS ******************************************
4540 !
4541       END
4542
4543
4544
4545
4546 !=======================================================================
4547 !
4548 ! *** ISORROPIA CODE
4549 ! *** SUBROUTINE ISRPIA
4550 ! *** THIS SUBROUTINE IS THE MASTER ROUTINE FOR THE ISORROPIA-PLUS
4551 !     THERMODYNAMIC EQUILIBRIUM AEROSOL MODEL (VERSIONS 0.x)
4552 !
4553 ! *** NOTE: THIS SUBROUTINE IS INCLUDED FOR BACKWARD COMPATABILITY ONLY.
4554 !     A PROGRAMMER SHOULD USE THE MORE COMPLETE SUBROUTINE ISOROPIA INSTEAD
4555 !
4556 !
4557 !     DEPENDING ON THE INPUT VALUES PROVIDED, THE FOLLOWING MODEL
4558 !     SUBVERSIONS ARE CALLED:
4559 !
4560 !     FOREWARD PROBLEM (IPROB=0):
4561 !     Na      SO4      NH4       NO3      CL       SUBROUTINE CALLED
4562 !     ----    ----     ----      ----     ----     -----------------
4563 !     0.0     >0.0     >0.0       0.0      0.0     SUBROUTINE ISRP1F
4564 !     0.0     >0.0     >0.0      >0.0      0.0     SUBROUTINE ISRP2F
4565 !     >0.0    >0.0     >0.0      >0.0     >0.0     SUBROUTINE ISRP3F
4566 !
4567 !     REVERSE PROBLEM (IPROB=1):
4568 !     Na      SO4      NH4       NO3      CL       SUBROUTINE CALLED
4569 !     ----    ----     ----      ----     ----     -----------------
4570 !     0.0     >0.0     >0.0       0.0      0.0     SUBROUTINE ISRP1R
4571 !     0.0     >0.0     >0.0      >0.0      0.0     SUBROUTINE ISRP2R
4572 !     >0.0    >0.0     >0.0      >0.0     >0.0     SUBROUTINE ISRP3R
4573 !
4574 ! *** COPYRIGHT 1996-2006, UNIVERSITY OF MIAMI, CARNEGIE MELLON UNIVERSITY,
4575 ! *** GEORGIA INSTITUTE OF TECHNOLOGY
4576 ! *** WRITTEN BY ATHANASIOS NENES
4577 ! *** UPDATED BY CHRISTOS FOUNTOUKIS
4578 !
4579 !=======================================================================
4580 !
4581 !      SUBROUTINE ISRPIA (WI, RHI, TEMPI, IPROBI)
4582       SUBROUTINE ISRPIAA (WI, RHI, TEMPI, IPROBI)
4583       USE ISRPIA
4584       IMPLICIT DOUBLE PRECISION (A-H, O-Z)
4585 !      INCLUDE 'ISRPIA.INC'
4586       DIMENSION WI(NCOMP)
4587 !
4588 ! *** PROBLEM TYPE (0=FOREWARD, 1=REVERSE) ******************************
4589 !
4590       IPROB = IPROBI
4591 !
4592 ! *** SOLVE FOREWARD PROBLEM ********************************************
4593 !
4594       IF (IPROB.EQ.0) THEN
4595          IF (WI(1)+WI(2)+WI(3)+WI(4)+WI(5) .LE. TINY) THEN ! Everything=0
4596             CALL INIT1 (WI, RHI, TEMPI)
4597          ELSE IF (WI(1)+WI(4)+WI(5) .LE. TINY) THEN        ! Na,Cl,NO3=0
4598             CALL ISRP1F (WI, RHI, TEMPI)
4599          ELSE IF (WI(1)+WI(5) .LE. TINY) THEN              ! Na,Cl=0
4600             CALL ISRP2F (WI, RHI, TEMPI)
4601          ELSE
4602             CALL ISRP3F (WI, RHI, TEMPI)
4603          ENDIF
4604 !
4605 ! *** SOLVE REVERSE PROBLEM *********************************************
4606 !
4607       ELSE
4608          IF (WI(1)+WI(2)+WI(3)+WI(4)+WI(5) .LE. TINY) THEN ! Everything=0
4609             CALL INIT1 (WI, RHI, TEMPI)
4610          ELSE IF (WI(1)+WI(4)+WI(5) .LE. TINY) THEN        ! Na,Cl,NO3=0
4611             CALL ISRP1R (WI, RHI, TEMPI)
4612          ELSE IF (WI(1)+WI(5) .LE. TINY) THEN              ! Na,Cl=0
4613             CALL ISRP2R (WI, RHI, TEMPI)
4614          ELSE
4615             CALL ISRP3R (WI, RHI, TEMPI)
4616          ENDIF
4617       ENDIF
4618 !
4619 ! *** SETUP 'DRY' FLAG ***************************************************
4620 !
4621       DRYF = WATER.LE.TINY
4622 !
4623 ! *** FIND TOTALS *******************************************************
4624 !
4625       IF (IPROB.EQ.0) THEN
4626          CONTINUE
4627       ELSE
4628          W(1) = WAER(1)
4629          W(2) = WAER(2)
4630          W(3) = WAER(3)
4631          W(4) = WAER(4)
4632          W(5) = WAER(5)
4633 !
4634          IF (.NOT.DRYF) THEN
4635             W(3) = W(3) + GNH3
4636             W(4) = W(4) + GHNO3
4637             W(5) = W(5) + GHCL
4638          ENDIF
4639       ENDIF
4640 !
4641       RETURN
4642 !
4643 ! *** END OF SUBROUTINE ISRPIA *******************************************
4644 !
4645       END
4646 !=======================================================================
4647 !
4648 ! *** ISORROPIA CODE
4649 ! *** SUBROUTINE PUSHERR
4650 ! *** THIS SUBROUTINE SAVES AN ERROR MESSAGE IN THE ERROR STACK
4651 !
4652 ! *** COPYRIGHT 1996-2006, UNIVERSITY OF MIAMI, CARNEGIE MELLON UNIVERSITY,
4653 ! *** GEORGIA INSTITUTE OF TECHNOLOGY
4654 ! *** WRITTEN BY ATHANASIOS NENES
4655 ! *** UPDATED BY CHRISTOS FOUNTOUKIS
4656 !
4657 !=======================================================================
4658 !
4659       SUBROUTINE PUSHERR (IERR,ERRINF)
4660       USE ISRPIA
4661       IMPLICIT DOUBLE PRECISION (A-H, O-Z)
4662 !      INCLUDE 'ISRPIA.INC'
4663       CHARACTER ERRINF*(*)
4664 !
4665 ! *** SAVE ERROR CODE IF THERE IS ANY SPACE ***************************
4666 !
4667       IF (NOFER.LT.NERRMX) THEN
4668          NOFER         = NOFER + 1
4669          ERRSTK(NOFER) = IERR
4670          ERRMSG(NOFER) = ERRINF
4671          STKOFL        =.FALSE.
4672       ELSE
4673          STKOFL        =.TRUE.      ! STACK OVERFLOW
4674       ENDIF
4675 !
4676 ! *** END OF SUBROUTINE PUSHERR ****************************************
4677 !
4678       END
4679
4680
4681
4682 !=======================================================================
4683 !
4684 ! *** ISORROPIA CODE
4685 ! *** SUBROUTINE ISERRINF
4686 ! *** THIS SUBROUTINE OBTAINS A COPY OF THE ERROR STACK (& MESSAGES)
4687 !
4688 ! *** COPYRIGHT 1996-2006, UNIVERSITY OF MIAMI, CARNEGIE MELLON UNIVERSITY,
4689 ! *** GEORGIA INSTITUTE OF TECHNOLOGY
4690 ! *** WRITTEN BY ATHANASIOS NENES
4691 ! *** UPDATED BY CHRISTOS FOUNTOUKIS
4692 !
4693 !=======================================================================
4694 !
4695       SUBROUTINE ISERRINF (ERRSTKI, ERRMSGI, NOFERI, STKOFLI)
4696       USE ISRPIA
4697       IMPLICIT DOUBLE PRECISION (A-H, O-Z)
4698 !      INCLUDE 'ISRPIA.INC'
4699       CHARACTER ERRMSGI*40
4700       INTEGER   ERRSTKI
4701       LOGICAL   STKOFLI
4702       DIMENSION ERRMSGI(NERRMX), ERRSTKI(NERRMX)
4703 !
4704 ! *** OBTAIN WHOLE ERROR STACK ****************************************
4705 !
4706       DO 10 I=1,NOFER              ! Error messages & codes
4707         ERRSTKI(I) = ERRSTK(I)
4708         ERRMSGI(I) = ERRMSG(I)
4709   10  CONTINUE
4710 !
4711       STKOFLI = STKOFL
4712       NOFERI  = NOFER
4713 !
4714       RETURN
4715 !
4716 ! *** END OF SUBROUTINE ISERRINF ***************************************
4717 !
4718       END
4719
4720
4721
4722 !=======================================================================
4723 !
4724 ! *** ISORROPIA CODE
4725 ! *** SUBROUTINE ERRSTAT
4726 ! *** THIS SUBROUTINE REPORTS ERROR MESSAGES TO UNIT 'IO'
4727 !
4728 ! *** COPYRIGHT 1996-2006, UNIVERSITY OF MIAMI, CARNEGIE MELLON UNIVERSITY,
4729 ! *** GEORGIA INSTITUTE OF TECHNOLOGY
4730 ! *** WRITTEN BY ATHANASIOS NENES
4731 ! *** UPDATED BY CHRISTOS FOUNTOUKIS
4732 !
4733 !=======================================================================
4734 !
4735       SUBROUTINE ERRSTAT (IO,IERR,ERRINF)
4736       USE ISRPIA
4737       IMPLICIT DOUBLE PRECISION (A-H, O-Z)
4738 !      INCLUDE 'ISRPIA.INC'
4739       CHARACTER CER*4, NCIS*29, NCIF*27, NSIS*26, NSIF*24, ERRINF*(*)
4740       DATA NCIS /'NO CONVERGENCE IN SUBROUTINE '/,   &
4741            NCIF /'NO CONVERGENCE IN FUNCTION '  /,   &
4742            NSIS /'NO SOLUTION IN SUBROUTINE '   /,   &
4743            NSIF /'NO SOLUTION IN FUNCTION '     /
4744 !
4745 ! *** WRITE ERROR IN CHARACTER *****************************************
4746 !
4747       WRITE (CER,'(I4)') IERR
4748       CALL RPLSTR (CER, ' ', '0',IOK)   ! REPLACE BLANKS WITH ZEROS
4749       CALL CHRBLN (ERRINF, IEND)        ! LAST POSITION OF ERRINF CHAR
4750 !
4751 ! *** WRITE ERROR TYPE (FATAL, WARNING ) *******************************
4752 !
4753       IF (IERR.EQ.0) THEN
4754          WRITE (IO,1000) 'NO ERRORS DETECTED '
4755          GOTO 10
4756 !
4757       ELSE IF (IERR.LT.0) THEN
4758          WRITE (IO,1000) 'ERROR STACK EXHAUSTED '
4759          GOTO 10
4760 !
4761       ELSE IF (IERR.GT.1000) THEN
4762          WRITE (IO,1100) 'FATAL',CER
4763 !
4764       ELSE
4765          WRITE (IO,1100) 'WARNING',CER
4766       ENDIF
4767 !
4768 ! *** WRITE ERROR MESSAGE **********************************************
4769 !
4770 ! FATAL MESSAGES
4771 !
4772       IF (IERR.EQ.1001) THEN
4773          CALL CHRBLN (SCASE, IEND)
4774          WRITE (IO,1000) 'CASE NOT SUPPORTED IN CALCMR ['//SCASE(1:IEND)   &
4775                          //']'
4776 !
4777       ELSEIF (IERR.EQ.1002) THEN
4778          CALL CHRBLN (SCASE, IEND)
4779          WRITE (IO,1000) 'CASE NOT SUPPORTED ['//SCASE(1:IEND)//']'
4780 !
4781 ! WARNING MESSAGES
4782 !
4783       ELSEIF (IERR.EQ.0001) THEN
4784          WRITE (IO,1000) NSIS,ERRINF
4785 !
4786       ELSEIF (IERR.EQ.0002) THEN
4787          WRITE (IO,1000) NCIS,ERRINF
4788 !
4789       ELSEIF (IERR.EQ.0003) THEN
4790          WRITE (IO,1000) NSIF,ERRINF
4791 !
4792       ELSEIF (IERR.EQ.0004) THEN
4793          WRITE (IO,1000) NCIF,ERRINF
4794 !
4795       ELSE IF (IERR.EQ.0019) THEN
4796          WRITE (IO,1000) 'HNO3(aq) AFFECTS H+, WHICH '//   &
4797                          'MIGHT AFFECT SO4/HSO4 RATIO'
4798          WRITE (IO,1000) 'DIRECT INCREASE IN H+ [',ERRINF(1:IEND),'] %'
4799 !
4800       ELSE IF (IERR.EQ.0020) THEN
4801          IF (W(4).GT.TINY .AND. W(5).GT.TINY) THEN
4802             WRITE (IO,1000) 'HSO4-SO4 EQUILIBRIUM MIGHT AFFECT HNO3,'   &
4803                           //'HCL DISSOLUTION'
4804          ELSE
4805             WRITE (IO,1000) 'HSO4-SO4 EQUILIBRIUM MIGHT AFFECT NH3 '   &
4806                           //'DISSOLUTION'
4807          ENDIF
4808          WRITE (IO,1000) 'DIRECT DECREASE IN H+ [',ERRINF(1:IEND),'] %'
4809 !
4810       ELSE IF (IERR.EQ.0021) THEN
4811          WRITE (IO,1000) 'HNO3(aq),HCL(aq) AFFECT H+, WHICH '//   &
4812                          'MIGHT AFFECT SO4/HSO4 RATIO'
4813          WRITE (IO,1000) 'DIRECT INCREASE IN H+ [',ERRINF(1:IEND),'] %'
4814 !
4815       ELSE IF (IERR.EQ.0022) THEN
4816          WRITE (IO,1000) 'HCL(g) EQUILIBRIUM YIELDS NONPHYSICAL '//   &
4817                          'DISSOLUTION'
4818          WRITE (IO,1000) 'A TINY AMOUNT [',ERRINF(1:IEND),'] IS '//   &
4819                          'ASSUMED TO BE DISSOLVED'
4820 !
4821       ELSEIF (IERR.EQ.0033) THEN
4822          WRITE (IO,1000) 'HCL(aq) AFFECTS H+, WHICH '//   &
4823                          'MIGHT AFFECT SO4/HSO4 RATIO'
4824          WRITE (IO,1000) 'DIRECT INCREASE IN H+ [',ERRINF(1:IEND),'] %'
4825 !
4826       ELSEIF (IERR.EQ.0050) THEN
4827          WRITE (IO,1000) 'TOO MUCH SODIUM GIVEN AS INPUT.'
4828          WRITE (IO,1000) 'REDUCED TO COMPLETELY NEUTRALIZE SO4,Cl,NO3.'
4829          WRITE (IO,1000) 'EXCESS SODIUM IS IGNORED.'
4830 !
4831       ELSE
4832          WRITE (IO,1000) 'NO DIAGNOSTIC MESSAGE AVAILABLE'
4833       ENDIF
4834 !
4835 10    RETURN
4836 !
4837 ! *** FORMAT STATEMENTS *************************************
4838 !
4839 1000  FORMAT (1X,A:A:A:A:A)
4840 1100  FORMAT (1X,A,' ERROR [',A4,']:')
4841 !
4842 ! *** END OF SUBROUTINE ERRSTAT *****************************
4843 !
4844       END
4845 !=======================================================================
4846 !
4847 ! *** ISORROPIA CODE
4848 ! *** SUBROUTINE ISORINF
4849 ! *** THIS SUBROUTINE PROVIDES INFORMATION ABOUT ISORROPIA
4850 !
4851 ! ======================== ARGUMENTS / USAGE ===========================
4852 !
4853 !  OUTPUT:
4854 !  1. [VERSI]
4855 !     CHARACTER*15 variable.
4856 !     Contains version-date information of ISORROPIA
4857 !
4858 !  2. [NCMP]
4859 !     INTEGER variable.
4860 !     The number of components needed in input array WI
4861 !     (or, the number of major species accounted for by ISORROPIA)
4862 !
4863 !  3. [NION]
4864 !     INTEGER variable
4865 !     The number of ions considered in the aqueous phase
4866 !
4867 !  4. [NAQGAS]
4868 !     INTEGER variable
4869 !     The number of undissociated species found in aqueous aerosol
4870 !     phase
4871 !
4872 !  5. [NSOL]
4873 !     INTEGER variable
4874 !     The number of solids considered in the solid aerosol phase
4875 !
4876 !  6. [NERR]
4877 !     INTEGER variable
4878 !     The size of the error stack (maximum number of errors that can
4879 !     be stored before the stack exhausts).
4880 !
4881 !  7. [TIN]
4882 !     DOUBLE PRECISION variable
4883 !     The value used for a very small number.
4884 !
4885 !  8. [GRT]
4886 !     DOUBLE PRECISION variable
4887 !     The value used for a very large number.
4888 !
4889 ! *** COPYRIGHT 1996-2006, UNIVERSITY OF MIAMI, CARNEGIE MELLON UNIVERSITY,
4890 ! *** GEORGIA INSTITUTE OF TECHNOLOGY
4891 ! *** WRITTEN BY ATHANASIOS NENES
4892 ! *** UPDATED BY CHRISTOS FOUNTOUKIS
4893 !
4894 !=======================================================================
4895 !
4896       SUBROUTINE ISORINF (VERSI, NCMP, NION, NAQGAS, NSOL, NERR, TIN,   &
4897                           GRT)
4898       USE ISRPIA
4899       IMPLICIT DOUBLE PRECISION (A-H, O-Z)
4900 !      INCLUDE 'ISRPIA.INC'
4901       CHARACTER VERSI*(*)
4902 !
4903 ! *** ASSIGN INFO *******************************************************
4904 !
4905       VERSI  = VERSION
4906       NCMP   = NCOMP
4907       NION   = NIONS
4908       NAQGAS = NGASAQ
4909       NSOL   = NSLDS
4910       NERR   = NERRMX
4911       TIN    = TINY
4912       GRT    = GREAT
4913 !
4914       RETURN
4915 !
4916 ! *** END OF SUBROUTINE ISORINF *******************************************
4917 !
4918       END