Merge remote-tracking branch 'origin/release-v4.6.1'
[WRF.git] / chem / isofwd.F
blobae6b2db9be24ab054a8ef0193d19c629b84b2a8f
1 !=======================================================================
3 ! *** ISORROPIA CODE
4 ! *** SUBROUTINE ISRP1F
5 ! *** THIS SUBROUTINE IS THE DRIVER ROUTINE FOR THE FOREWARD PROBLEM OF
6 !     AN AMMONIUM-SULFATE AEROSOL SYSTEM.
7 !     THE COMPOSITION REGIME IS DETERMINED BY THE SULFATE RATIO AND BY
8 !     THE AMBIENT RELATIVE HUMIDITY.
10 ! *** COPYRIGHT 1996-2006, UNIVERSITY OF MIAMI, CARNEGIE MELLON UNIVERSITY,
11 ! *** GEORGIA INSTITUTE OF TECHNOLOGY
12 ! *** WRITTEN BY ATHANASIOS NENES
13 ! *** UPDATED BY CHRISTOS FOUNTOUKIS
15 !  REVISION HISTORY:                                                   *
16 !    Original code was provided by Dr. ATHANASIOS NENES, Georgia Tech, in 2000
17 !    Revised by Y. Zhang, AER, Inc. to incorporate v1.5 into MADRID, 2000
18 !    Revised by Y. Zhang and Xiao-Ming Hu to incorporate it along with MADRID into WRF/Chem, 2005
19 !    Updated by Xiao-Ming Hu and Y. Zhang, NCSU to v. 1.7, Oct., 2007
20 !=======================================================================
22       SUBROUTINE ISRP1F (WI, RHI, TEMPI)
23       USE ISRPIA
24       IMPLICIT DOUBLE PRECISION (A-H, O-Z)
25 !      INCLUDE 'ISRPIA.INC'
26       DIMENSION WI(NCOMP)
28 ! *** INITIALIZE ALL VARIABLES IN COMMON BLOCK **************************
30       CALL INIT1 (WI, RHI, TEMPI)
32 ! *** CALCULATE SULFATE RATIO *******************************************
34       SULRAT = W(3)/W(2)
36 ! *** FIND CALCULATION REGIME FROM (SULRAT,RH) **************************
38 ! *** SULFATE POOR
40       IF (2.0.LE.SULRAT) THEN
41       DC   = W(3) - 2.001D0*W(2)  ! For numerical stability
42       W(3) = W(3) + MAX(-DC, ZERO)
44       IF(METSTBL.EQ.1) THEN
45          SCASE = 'A2'
46          CALL CALCA2                 ! Only liquid (metastable)
47       ELSE
49          IF (RH.LT.DRNH42S4) THEN
50             SCASE = 'A1'
51             CALL CALCA1              ! NH42SO4              ; case A1
53          ELSEIF (DRNH42S4.LE.RH) THEN
54             SCASE = 'A2'
55             CALL CALCA2              ! Only liquid          ; case A2
56          ENDIF
57       ENDIF
59 ! *** SULFATE RICH (NO ACID)
61       ELSEIF (1.0.LE.SULRAT .AND. SULRAT.LT.2.0) THEN
63       IF(METSTBL.EQ.1) THEN
64          SCASE = 'B4'
65          CALL CALCB4                 ! Only liquid (metastable)
66       ELSE
68          IF (RH.LT.DRNH4HS4) THEN
69             SCASE = 'B1'
70             CALL CALCB1              ! NH4HSO4,LC,NH42SO4   ; case B1
72          ELSEIF (DRNH4HS4.LE.RH .AND. RH.LT.DRLC) THEN
73             SCASE = 'B2'
74             CALL CALCB2              ! LC,NH42S4            ; case B2
76          ELSEIF (DRLC.LE.RH .AND. RH.LT.DRNH42S4) THEN
77             SCASE = 'B3'
78             CALL CALCB3              ! NH42S4               ; case B3
80          ELSEIF (DRNH42S4.LE.RH) THEN
81             SCASE = 'B4'
82             CALL CALCB4              ! Only liquid          ; case B4
83          ENDIF
84       ENDIF
85       CALL CALCNH3
87 ! *** SULFATE RICH (FREE ACID)
89       ELSEIF (SULRAT.LT.1.0) THEN
91       IF(METSTBL.EQ.1) THEN
92          SCASE = 'C2'
93          CALL CALCC2                 ! Only liquid (metastable)
94       ELSE
96          IF (RH.LT.DRNH4HS4) THEN
97             SCASE = 'C1'
98             CALL CALCC1              ! NH4HSO4              ; case C1
100          ELSEIF (DRNH4HS4.LE.RH) THEN
101             SCASE = 'C2'
102             CALL CALCC2              ! Only liquid          ; case C2
104          ENDIF
105       ENDIF
106       CALL CALCNH3
107       ENDIF
109 ! *** RETURN POINT
111       RETURN
113 ! *** END OF SUBROUTINE ISRP1F *****************************************
115       END SUBROUTINE ISRP1F                 
116 !=======================================================================
118 ! *** ISORROPIA CODE
119 ! *** SUBROUTINE ISRP2F
120 ! *** THIS SUBROUTINE IS THE DRIVER ROUTINE FOR THE FOREWARD PROBLEM OF
121 !     AN AMMONIUM-SULFATE-NITRATE AEROSOL SYSTEM.
122 !     THE COMPOSITION REGIME IS DETERMINED BY THE SULFATE RATIO AND BY
123 !     THE AMBIENT RELATIVE HUMIDITY.
125 ! *** COPYRIGHT 1996-2006, UNIVERSITY OF MIAMI, CARNEGIE MELLON UNIVERSITY,
126 ! *** GEORGIA INSTITUTE OF TECHNOLOGY
127 ! *** WRITTEN BY ATHANASIOS NENES
128 ! *** UPDATED BY CHRISTOS FOUNTOUKIS
130 !=======================================================================
132       SUBROUTINE ISRP2F (WI, RHI, TEMPI)
133       USE ISRPIA
134       IMPLICIT DOUBLE PRECISION (A-H, O-Z)
135 !      INCLUDE 'ISRPIA.INC'
136       DIMENSION WI(NCOMP)
138 ! *** INITIALIZE ALL VARIABLES IN COMMON BLOCK **************************
140       CALL INIT2 (WI, RHI, TEMPI)
142 ! *** CALCULATE SULFATE RATIO *******************************************
144       SULRAT = W(3)/W(2)
146 ! *** FIND CALCULATION REGIME FROM (SULRAT,RH) **************************
148 ! *** SULFATE POOR
150       IF (2.0.LE.SULRAT) THEN
152       IF(METSTBL.EQ.1) THEN
153          SCASE = 'D3'
154          CALL CALCD3                 ! Only liquid (metastable)
155       ELSE
157          IF (RH.LT.DRNH4NO3) THEN
158             SCASE = 'D1'
159             CALL CALCD1              ! NH42SO4,NH4NO3       ; case D1
161          ELSEIF (DRNH4NO3.LE.RH .AND. RH.LT.DRNH42S4) THEN
162             SCASE = 'D2'
163             CALL CALCD2              ! NH42S4               ; case D2
165          ELSEIF (DRNH42S4.LE.RH) THEN
166             SCASE = 'D3'
167             CALL CALCD3              ! Only liquid          ; case D3
168          ENDIF
169       ENDIF
171 ! *** SULFATE RICH (NO ACID)
172 !     FOR SOLVING THIS CASE, NITRIC ACID IS ASSUMED A MINOR SPECIES,
173 !     THAT DOES NOT SIGNIFICANTLY PERTURB THE HSO4-SO4 EQUILIBRIUM.
174 !     SUBROUTINES CALCB? ARE CALLED, AND THEN THE NITRIC ACID IS DISSOLVED
175 !     FROM THE HNO3(G) -> (H+) + (NO3-) EQUILIBRIUM.
177       ELSEIF (1.0.LE.SULRAT .AND. SULRAT.LT.2.0) THEN
179       IF(METSTBL.EQ.1) THEN
180          SCASE = 'B4'
181          CALL CALCB4                 ! Only liquid (metastable)
182          SCASE = 'E4'
183       ELSE
185          IF (RH.LT.DRNH4HS4) THEN
186             SCASE = 'B1'
187             CALL CALCB1              ! NH4HSO4,LC,NH42SO4   ; case E1
188             SCASE = 'E1'
190          ELSEIF (DRNH4HS4.LE.RH .AND. RH.LT.DRLC) THEN
191             SCASE = 'B2'
192             CALL CALCB2              ! LC,NH42S4            ; case E2
193             SCASE = 'E2'
195          ELSEIF (DRLC.LE.RH .AND. RH.LT.DRNH42S4) THEN
196             SCASE = 'B3'
197             CALL CALCB3              ! NH42S4               ; case E3
198             SCASE = 'E3'
200          ELSEIF (DRNH42S4.LE.RH) THEN
201             SCASE = 'B4'
202             CALL CALCB4              ! Only liquid          ; case E4
203             SCASE = 'E4'
204          ENDIF
205       ENDIF
207       CALL CALCNA                 ! HNO3(g) DISSOLUTION
209 ! *** SULFATE RICH (FREE ACID)
210 !     FOR SOLVING THIS CASE, NITRIC ACID IS ASSUMED A MINOR SPECIES,
211 !     THAT DOES NOT SIGNIFICANTLY PERTURB THE HSO4-SO4 EQUILIBRIUM
212 !     SUBROUTINE CALCC? IS CALLED, AND THEN THE NITRIC ACID IS DISSOLVED
213 !     FROM THE HNO3(G) -> (H+) + (NO3-) EQUILIBRIUM.
215       ELSEIF (SULRAT.LT.1.0) THEN
217       IF(METSTBL.EQ.1) THEN
218          SCASE = 'C2'
219          CALL CALCC2                 ! Only liquid (metastable)
220          SCASE = 'F2'
221       ELSE
223          IF (RH.LT.DRNH4HS4) THEN
224             SCASE = 'C1'
225             CALL CALCC1              ! NH4HSO4              ; case F1
226             SCASE = 'F1'
228          ELSEIF (DRNH4HS4.LE.RH) THEN
229             SCASE = 'C2'
230             CALL CALCC2              ! Only liquid          ; case F2
231             SCASE = 'F2'
232          ENDIF
233       ENDIF
235       CALL CALCNA                 ! HNO3(g) DISSOLUTION
236       ENDIF
238 ! *** RETURN POINT
240       RETURN
242 ! *** END OF SUBROUTINE ISRP2F *****************************************
244       END SUBROUTINE ISRP2F                 
245 !=======================================================================
247 ! *** ISORROPIA CODE
248 ! *** SUBROUTINE ISRP3F
249 ! *** THIS SUBROUTINE IS THE DRIVER ROUTINE FOR THE FORWARD PROBLEM OF
250 !     AN AMMONIUM-SULFATE-NITRATE-CHLORIDE-SODIUM AEROSOL SYSTEM.
251 !     THE COMPOSITION REGIME IS DETERMINED BY THE SULFATE & SODIUM
252 !     RATIOS AND BY THE AMBIENT RELATIVE HUMIDITY.
254 ! *** COPYRIGHT 1996-2006, UNIVERSITY OF MIAMI, CARNEGIE MELLON UNIVERSITY,
255 ! *** GEORGIA INSTITUTE OF TECHNOLOGY
256 ! *** WRITTEN BY ATHANASIOS NENES
257 ! *** UPDATED BY CHRISTOS FOUNTOUKIS
259 !=======================================================================
261       SUBROUTINE ISRP3F (WI, RHI, TEMPI)
262       USE ISRPIA
263       IMPLICIT DOUBLE PRECISION (A-H, O-Z)
264 !      INCLUDE 'ISRPIA.INC'
265       DIMENSION WI(NCOMP)
267 ! *** ADJUST FOR TOO LITTLE AMMONIUM AND CHLORIDE ***********************
269       WI(3) = MAX (WI(3), 1.D-10)  ! NH4+ : 1e-4 umoles/m3
270       WI(5) = MAX (WI(5), 1.D-10)  ! Cl-  : 1e-4 umoles/m3
272 ! *** ADJUST FOR TOO LITTLE SODIUM, SULFATE AND NITRATE COMBINED ********
274       IF (WI(1)+WI(2)+WI(4) .LE. 1d-10) THEN
275          WI(1) = 1.D-10  ! Na+  : 1e-4 umoles/m3
276          WI(2) = 1.D-10  ! SO4- : 1e-4 umoles/m3
277       ENDIF
279 ! *** INITIALIZE ALL VARIABLES IN COMMON BLOCK **************************
281       CALL ISOINIT3 (WI, RHI, TEMPI)
283 ! *** CHECK IF TOO MUCH SODIUM ; ADJUST AND ISSUE ERROR MESSAGE *********
285       REST = 2.D0*W(2) + W(4) + W(5)
286       IF (W(1).GT.REST) THEN            ! NA > 2*SO4+CL+NO3 ?
287          W(1) = (ONE-1D-6)*REST         ! Adjust Na amount
288          CALL PUSHERR (0050, 'ISRP3F')  ! Warning error: Na adjusted
289       ENDIF
291 ! *** CALCULATE SULFATE & SODIUM RATIOS *********************************
293       SULRAT = (W(1)+W(3))/W(2)
294       SODRAT = W(1)/W(2)
296 ! *** FIND CALCULATION REGIME FROM (SULRAT,RH) **************************
298 ! *** SULFATE POOR ; SODIUM POOR
300       IF (2.0.LE.SULRAT .AND. SODRAT.LT.2.0) THEN
302       IF(METSTBL.EQ.1) THEN
303          SCASE = 'G5'
304          CALL CALCG5                 ! Only liquid (metastable)
305       ELSE
307          IF (RH.LT.DRNH4NO3) THEN
308             SCASE = 'G1'
309             CALL CALCG1              ! NH42SO4,NH4NO3,NH4CL,NA2SO4
311          ELSEIF (DRNH4NO3.LE.RH .AND. RH.LT.DRNH4CL) THEN
312             SCASE = 'G2'
313             CALL CALCG2              ! NH42SO4,NH4CL,NA2SO4
315          ELSEIF (DRNH4CL.LE.RH  .AND. RH.LT.DRNH42S4) THEN
316             SCASE = 'G3'
317             CALL CALCG3              ! NH42SO4,NA2SO4
319         ELSEIF (DRNH42S4.LE.RH  .AND. RH.LT.DRNA2SO4) THEN
320             SCASE = 'G4'
321             CALL CALCG4              ! NA2SO4
323          ELSEIF (DRNA2SO4.LE.RH) THEN
324             SCASE = 'G5'
325             CALL CALCG5              ! Only liquid
326          ENDIF
327       ENDIF
329 ! *** SULFATE POOR ; SODIUM RICH
331       ELSE IF (SULRAT.GE.2.0 .AND. SODRAT.GE.2.0) THEN
333       IF(METSTBL.EQ.1) THEN
334          SCASE = 'H6'
335          CALL CALCH6                 ! Only liquid (metastable)
336       ELSE
338          IF (RH.LT.DRNH4NO3) THEN
339             SCASE = 'H1'
340             CALL CALCH1              ! NH4NO3,NH4CL,NA2SO4,NACL,NANO3
342          ELSEIF (DRNH4NO3.LE.RH .AND. RH.LT.DRNANO3) THEN
343             SCASE = 'H2'
344             CALL CALCH2              ! NH4CL,NA2SO4,NACL,NANO3
346          ELSEIF (DRNANO3.LE.RH  .AND. RH.LT.DRNACL) THEN
347             SCASE = 'H3'
348             CALL CALCH3              ! NH4CL,NA2SO4,NACL
350          ELSEIF (DRNACL.LE.RH   .AND. RH.LT.DRNH4Cl) THEN
351             SCASE = 'H4'
352             CALL CALCH4              ! NH4CL,NA2SO4
354          ELSEIF (DRNH4Cl.LE.RH .AND. RH.LT.DRNA2SO4) THEN
355             SCASE = 'H5'
356             CALL CALCH5              ! NA2SO4
358          ELSEIF (DRNA2SO4.LE.RH) THEN
359             SCASE = 'H6'
360             CALL CALCH6              ! NO SOLID
361          ENDIF
362       ENDIF
364 ! *** SULFATE RICH (NO ACID)
366       ELSEIF (1.0.LE.SULRAT .AND. SULRAT.LT.2.0) THEN
368       IF(METSTBL.EQ.1) THEN
369          SCASE = 'I6'
370          CALL CALCI6                 ! Only liquid (metastable)
371       ELSE
373          IF (RH.LT.DRNH4HS4) THEN
374             SCASE = 'I1'
375             CALL CALCI1              ! NA2SO4,(NH4)2SO4,NAHSO4,NH4HSO4,LC
377          ELSEIF (DRNH4HS4.LE.RH .AND. RH.LT.DRNAHSO4) THEN
378             SCASE = 'I2'
379             CALL CALCI2              ! NA2SO4,(NH4)2SO4,NAHSO4,LC
381          ELSEIF (DRNAHSO4.LE.RH .AND. RH.LT.DRLC) THEN
382             SCASE = 'I3'
383             CALL CALCI3              ! NA2SO4,(NH4)2SO4,LC
385          ELSEIF (DRLC.LE.RH     .AND. RH.LT.DRNH42S4) THEN
386             SCASE = 'I4'
387             CALL CALCI4              ! NA2SO4,(NH4)2SO4
389          ELSEIF (DRNH42S4.LE.RH .AND. RH.LT.DRNA2SO4) THEN
390             SCASE = 'I5'
391             CALL CALCI5              ! NA2SO4
393          ELSEIF (DRNA2SO4.LE.RH) THEN
394             SCASE = 'I6'
395             CALL CALCI6              ! NO SOLIDS
396          ENDIF
397       ENDIF
399       CALL CALCNHA                ! MINOR SPECIES: HNO3, HCl
400       CALL CALCNH3                !                NH3
402 ! *** SULFATE RICH (FREE ACID)
404       ELSEIF (SULRAT.LT.1.0) THEN
406       IF(METSTBL.EQ.1) THEN
407          SCASE = 'J3'
408          CALL CALCJ3                 ! Only liquid (metastable)
409       ELSE
411          IF (RH.LT.DRNH4HS4) THEN
412             SCASE = 'J1'
413             CALL CALCJ1              ! NH4HSO4,NAHSO4
415          ELSEIF (DRNH4HS4.LE.RH .AND. RH.LT.DRNAHSO4) THEN
416             SCASE = 'J2'
417             CALL CALCJ2              ! NAHSO4
419          ELSEIF (DRNAHSO4.LE.RH) THEN
420             SCASE = 'J3'
421             CALL CALCJ3
422          ENDIF
423       ENDIF
425       CALL CALCNHA                ! MINOR SPECIES: HNO3, HCl
426       CALL CALCNH3                !                NH3
427       ENDIF
429 ! *** RETURN POINT
431       RETURN
433 ! *** END OF SUBROUTINE ISRP3F *****************************************
435       END SUBROUTINE ISRP3F                 
436 !=======================================================================
438 ! *** ISORROPIA CODE
439 ! *** SUBROUTINE CALCA2
440 ! *** CASE A2
442 !     THE MAIN CHARACTERISTICS OF THIS REGIME ARE:
443 !     1. SULFATE POOR (SULRAT > 2.0)
444 !     2. LIQUID AEROSOL PHASE ONLY POSSIBLE
446 !     FOR CALCULATIONS, A BISECTION IS PERFORMED TOWARDS X, THE
447 !     AMOUNT OF HYDROGEN IONS (H+) FOUND IN THE LIQUID PHASE.
448 !     FOR EACH ESTIMATION OF H+, FUNCTION FUNCB2A CALCULATES THE
449 !     CONCENTRATION OF IONS FROM THE NH3(GAS) - NH4+(LIQ) EQUILIBRIUM.
450 !     ELECTRONEUTRALITY IS USED AS THE OBJECTIVE FUNCTION.
452 ! *** COPYRIGHT 1996-2006, UNIVERSITY OF MIAMI, CARNEGIE MELLON UNIVERSITY,
453 ! *** GEORGIA INSTITUTE OF TECHNOLOGY
454 ! *** WRITTEN BY ATHANASIOS NENES
455 ! *** UPDATED BY CHRISTOS FOUNTOUKIS
457 !=======================================================================
459       SUBROUTINE CALCA2
460       USE ISRPIA
461       IMPLICIT DOUBLE PRECISION (A-H, O-Z)
462 !      INCLUDE 'ISRPIA.INC'
464 ! *** SETUP PARAMETERS ************************************************
466       CALAOU    =.TRUE.       ! Outer loop activity calculation flag
467       OMELO     = TINY        ! Low  limit: SOLUTION IS VERY BASIC
468       OMEHI     = 2.0D0*W(2)  ! High limit: FROM NH4+ -> NH3(g) + H+(aq)
470 ! *** CALCULATE WATER CONTENT *****************************************
472       MOLAL(5) = W(2)
473       MOLAL(6) = ZERO
474       CALL CALCMR
476 ! *** INITIAL VALUES FOR BISECTION ************************************
478       X1 = OMEHI
479       Y1 = FUNCA2 (X1)
480       IF (ABS(Y1).LE.EPS) RETURN
482 ! *** ROOT TRACKING ; FOR THE RANGE OF HI AND LO **********************
484       DX = (OMEHI-OMELO)/REAL(NDIV)
485       DO 10 I=1,NDIV
486          X2 = MAX(X1-DX, OMELO)
487          Y2 = FUNCA2 (X2)
488          IF (SIGN(1.d0,Y1)*SIGN(1.d0,Y2).LT.ZERO) GOTO 20  ! (Y1*Y2.LT.ZERO)
489          X1 = X2
490          Y1 = Y2
491 10    CONTINUE
492       IF (ABS(Y2).LE.EPS) THEN
493          RETURN
494       ELSE
495          CALL PUSHERR (0001, 'CALCA2')    ! WARNING ERROR: NO SOLUTION
496          RETURN
497       ENDIF
499 ! *** PERFORM BISECTION ***********************************************
501 20    DO 30 I=1,MAXIT
502          X3 = 0.5*(X1+X2)
503          Y3 = FUNCA2 (X3)
504          IF (SIGN(1.d0,Y1)*SIGN(1.d0,Y3) .LE. ZERO) THEN  ! (Y1*Y3 .LE. ZERO)
505             Y2    = Y3
506             X2    = X3
507          ELSE
508             Y1    = Y3
509             X1    = X3
510          ENDIF
511          IF (ABS(X2-X1) .LE. EPS*X1) GOTO 40
512 30    CONTINUE
513       CALL PUSHERR (0002, 'CALCA2')    ! WARNING ERROR: NO CONVERGENCE
515 ! *** CONVERGED ; RETURN **********************************************
517 40    X3 = 0.5*(X1+X2)
518       Y3 = FUNCA2 (X3)
519       RETURN
521 ! *** END OF SUBROUTINE CALCA2 ****************************************
523       END SUBROUTINE CALCA2
527 !=======================================================================
529 ! *** ISORROPIA CODE
530 ! *** FUNCTION FUNCA2
531 ! *** CASE A2
532 !     FUNCTION THAT SOLVES THE SYSTEM OF EQUATIONS FOR CASE A2 ;
533 !     AND RETURNS THE VALUE OF THE ZEROED FUNCTION IN FUNCA2.
535 !=======================================================================
537       DOUBLE PRECISION FUNCTION FUNCA2 (OMEGI)
538       USE ISRPIA
539       IMPLICIT DOUBLE PRECISION (A-H, O-Z)
540 !      INCLUDE 'ISRPIA.INC'
541       DOUBLE PRECISION LAMDA
543 ! *** SETUP PARAMETERS ************************************************
545       FRST   = .TRUE.
546       CALAIN = .TRUE.
547       PSI    = W(2)         ! INITIAL AMOUNT OF (NH4)2SO4 IN SOLUTION
549 ! *** SOLVE EQUATIONS ; WITH ITERATIONS FOR ACTIVITY COEF. ************
551       DO 10 I=1,NSWEEP
552          A1    = XK1*WATER/GAMA(7)*(GAMA(8)/GAMA(7))**2.
553          A2    = XK2*R*TEMP/XKW*(GAMA(8)/GAMA(9))**2.
554          A3    = XKW*RH*WATER*WATER
556          LAMDA = PSI/(A1/OMEGI+ONE)
557          ZETA  = A3/OMEGI
559 ! *** SPECIATION & WATER CONTENT ***************************************
561          MOLAL (1) = OMEGI                        ! HI
562          MOLAL (3) = W(3)/(ONE/A2/OMEGI + ONE)    ! NH4I
563          MOLAL (5) = MAX(PSI-LAMDA,TINY)          ! SO4I
564          MOLAL (6) = LAMDA                        ! HSO4I
565          GNH3      = MAX (W(3)-MOLAL(3), TINY)    ! NH3GI
566          COH       = ZETA                         ! OHI
568 ! *** CALCULATE ACTIVITIES OR TERMINATE INTERNAL LOOP *****************
570          IF (FRST.AND.CALAOU .OR. .NOT.FRST.AND.CALAIN) THEN
571             CALL CALCACT
572          ELSE
573             GOTO 20
574          ENDIF
575 10    CONTINUE
577 ! *** CALCULATE OBJECTIVE FUNCTION ************************************
579 20    DENOM = (2.0*MOLAL(5)+MOLAL(6))
580       FUNCA2= (MOLAL(3)/DENOM - ONE) + MOLAL(1)/DENOM
581       RETURN
583 ! *** END OF FUNCTION FUNCA2 ********************************************
585       END FUNCTION FUNCA2        
586 !=======================================================================
588 ! *** ISORROPIA CODE
589 ! *** SUBROUTINE CALCA1
590 ! *** CASE A1
592 !     THE MAIN CHARACTERISTICS OF THIS REGIME ARE:
593 !     1. SULFATE POOR (SULRAT > 2.0)
594 !     2. SOLID AEROSOL ONLY
595 !     3. SOLIDS POSSIBLE : (NH4)2SO4
597 !     A SIMPLE MATERIAL BALANCE IS PERFORMED, AND THE SOLID (NH4)2SO4
598 !     IS CALCULATED FROM THE SULFATES. THE EXCESS AMMONIA REMAINS IN
599 !     THE GAS PHASE.
601 ! *** COPYRIGHT 1996-2006, UNIVERSITY OF MIAMI, CARNEGIE MELLON UNIVERSITY,
602 ! *** GEORGIA INSTITUTE OF TECHNOLOGY
603 ! *** WRITTEN BY ATHANASIOS NENES
604 ! *** UPDATED BY CHRISTOS FOUNTOUKIS
606 !=======================================================================
608       SUBROUTINE CALCA1
609       USE ISRPIA
610       IMPLICIT DOUBLE PRECISION (A-H, O-Z)
611 !      INCLUDE 'ISRPIA.INC'
613       CNH42S4 = W(2)
614       GNH3    = MAX (W(3)-2.0*CNH42S4, ZERO)
615       RETURN
617 ! *** END OF SUBROUTINE CALCA1 ******************************************
619       END SUBROUTINE CALCA1
623 !=======================================================================
625 ! *** ISORROPIA CODE
626 ! *** SUBROUTINE CALCB4
627 ! *** CASE B4
629 !     THE MAIN CHARACTERISTICS OF THIS REGIME ARE:
630 !     1. SULFATE RICH, NO FREE ACID (1.0 <= SULRAT < 2.0)
631 !     2. LIQUID AEROSOL PHASE ONLY POSSIBLE
633 !     FOR CALCULATIONS, A BISECTION IS PERFORMED WITH RESPECT TO H+.
634 !     THE OBJECTIVE FUNCTION IS THE DIFFERENCE BETWEEN THE ESTIMATED H+
635 !     AND THAT CALCULATED FROM ELECTRONEUTRALITY.
637 ! *** COPYRIGHT 1996-2006, UNIVERSITY OF MIAMI, CARNEGIE MELLON UNIVERSITY,
638 ! *** GEORGIA INSTITUTE OF TECHNOLOGY
639 ! *** WRITTEN BY ATHANASIOS NENES
640 ! *** UPDATED BY CHRISTOS FOUNTOUKIS
642 !=======================================================================
644       SUBROUTINE CALCB4
645       USE ISRPIA
646       IMPLICIT DOUBLE PRECISION (A-H, O-Z)
647 !      INCLUDE 'ISRPIA.INC'
649 ! *** SOLVE EQUATIONS **************************************************
651       FRST       = .TRUE.
652       CALAIN     = .TRUE.
653       CALAOU     = .TRUE.
655 ! *** CALCULATE WATER CONTENT ******************************************
657       CALL CALCB1A         ! GET DRY SALT CONTENT, AND USE FOR WATER.
658       MOLALR(13) = CLC
659       MOLALR(9)  = CNH4HS4
660       MOLALR(4)  = CNH42S4
661       CLC        = ZERO
662       CNH4HS4    = ZERO
663       CNH42S4    = ZERO
664       WATER      = MOLALR(13)/M0(13)+MOLALR(9)/M0(9)+MOLALR(4)/M0(4)
666       MOLAL(3)   = W(3)   ! NH4I
668       DO 20 I=1,NSWEEP
669          AK1   = XK1*((GAMA(8)/GAMA(7))**2.)*(WATER/GAMA(7))
670          BET   = W(2)
671          GAM   = MOLAL(3)
673          BB    = BET + AK1 - GAM
674          CC    =-AK1*BET
675          DD    = BB*BB - 4.D0*CC
677 ! *** SPECIATION & WATER CONTENT ***************************************
679          MOLAL (5) = MAX(TINY,MIN(0.5*(-BB + SQRT(DD)), W(2))) ! SO4I
680          MOLAL (6) = MAX(TINY,MIN(W(2)-MOLAL(5),W(2)))         ! HSO4I
681          MOLAL (1) = MAX(TINY,MIN(AK1*MOLAL(6)/MOLAL(5),W(2))) ! HI
682          CALL CALCMR                                           ! Water content
684 ! *** CALCULATE ACTIVITIES OR TERMINATE INTERNAL LOOP *****************
686          IF (.NOT.CALAIN) GOTO 30
687          CALL CALCACT
688 20    CONTINUE
690 30    RETURN
692 ! *** END OF SUBROUTINE CALCB4 ******************************************
694       END SUBROUTINE CALCB4
695 !=======================================================================
697 ! *** ISORROPIA CODE
698 ! *** SUBROUTINE CALCB3
699 ! *** CASE B3
701 !     THE MAIN CHARACTERISTICS OF THIS REGIME ARE:
702 !     1. SULFATE RICH, NO FREE ACID (1.0 <= SULRAT < 2.0)
703 !     2. BOTH LIQUID & SOLID PHASE IS POSSIBLE
704 !     3. SOLIDS POSSIBLE: (NH4)2SO4
706 ! *** COPYRIGHT 1996-2006, UNIVERSITY OF MIAMI, CARNEGIE MELLON UNIVERSITY,
707 ! *** GEORGIA INSTITUTE OF TECHNOLOGY
708 ! *** WRITTEN BY ATHANASIOS NENES
709 ! *** UPDATED BY CHRISTOS FOUNTOUKIS
711 !=======================================================================
713       SUBROUTINE CALCB3
714       USE ISRPIA
715       IMPLICIT DOUBLE PRECISION (A-H, O-Z)
716 !      INCLUDE 'ISRPIA.INC'
718 ! *** CALCULATE EQUIVALENT AMOUNT OF HSO4 AND SO4 ***********************
720       X = MAX(2*W(2)-W(3), ZERO)   ! Equivalent NH4HSO4
721       Y = MAX(W(3)  -W(2), ZERO)   ! Equivalent NH42SO4
723 ! *** CALCULATE SPECIES ACCORDING TO RELATIVE ABUNDANCE OF HSO4 *********
725       IF (X.LT.Y) THEN             ! LC is the MIN (x,y)
726          SCASE   = 'B3 ; SUBCASE 1'
727          TLC     = X
728          TNH42S4 = Y-X
729          CALL CALCB3A (TLC,TNH42S4)      ! LC + (NH4)2SO4
730       ELSE
731          SCASE   = 'B3 ; SUBCASE 2'
732          TLC     = Y
733          TNH4HS4 = X-Y
734          CALL CALCB3B (TLC,TNH4HS4)      ! LC + NH4HSO4
735       ENDIF
737       RETURN
739 ! *** END OF SUBROUTINE CALCB3 ******************************************
741       END SUBROUTINE CALCB3
744 !=======================================================================
746 ! *** ISORROPIA CODE
747 ! *** SUBROUTINE CALCB3A
748 ! *** CASE B3 ; SUBCASE 1
750 !     THE MAIN CHARACTERISTICS OF THIS REGIME ARE:
751 !     1. SULFATE RICH (1.0 < SULRAT < 2.0)
752 !     2. BOTH LIQUID & SOLID PHASE IS POSSIBLE
753 !     3. SOLIDS POSSIBLE: (NH4)2SO4
755 !     FOR CALCULATIONS, A BISECTION IS PERFORMED TOWARDS ZETA, THE
756 !     AMOUNT OF SOLID (NH4)2SO4 DISSOLVED IN THE LIQUID PHASE.
757 !     FOR EACH ESTIMATION OF ZETA, FUNCTION FUNCB3A CALCULATES THE
758 !     AMOUNT OF H+ PRODUCED (BASED ON THE SO4 RELEASED INTO THE
759 !     SOLUTION). THE SOLUBILITY PRODUCT OF (NH4)2SO4 IS USED AS THE
760 !     OBJECTIVE FUNCTION.
762 ! *** COPYRIGHT 1996-2006, UNIVERSITY OF MIAMI, CARNEGIE MELLON UNIVERSITY,
763 ! *** GEORGIA INSTITUTE OF TECHNOLOGY
764 ! *** WRITTEN BY ATHANASIOS NENES
765 ! *** UPDATED BY CHRISTOS FOUNTOUKIS
767 !=======================================================================
769       SUBROUTINE CALCB3A (TLC, TNH42S4)
770       USE ISRPIA
771       IMPLICIT DOUBLE PRECISION (A-H, O-Z)
772 !      INCLUDE 'ISRPIA.INC'
774       CALAOU = .TRUE.         ! Outer loop activity calculation flag
775       ZLO    = ZERO           ! MIN DISSOLVED (NH4)2SO4
776       ZHI    = TNH42S4        ! MAX DISSOLVED (NH4)2SO4
778 ! *** INITIAL VALUES FOR BISECTION (DISSOLVED (NH4)2SO4 ****************
780       Z1 = ZLO
781       Y1 = FUNCB3A (Z1, TLC, TNH42S4)
782       IF (ABS(Y1).LE.EPS) RETURN
783       YLO= Y1
785 ! *** ROOT TRACKING ; FOR THE RANGE OF HI AND LO ***********************
787       DZ = (ZHI-ZLO)/REAL(NDIV)
788       DO 10 I=1,NDIV
789          Z2 = Z1+DZ
790          Y2 = FUNCB3A (Z2, TLC, TNH42S4)
791          IF (SIGN(1.d0,Y1)*SIGN(1.d0,Y2).LT.ZERO) GOTO 20  ! (Y1*Y2.LT.ZERO)
792          Z1 = Z2
793          Y1 = Y2
794 10    CONTINUE
796 ! *** NO SUBDIVISION WITH SOLUTION FOUND
798       YHI= Y1                      ! Save Y-value at HI position
799       IF (ABS(Y2) .LT. EPS) THEN   ! X2 IS A SOLUTION
800          RETURN
802 ! *** { YLO, YHI } < 0.0 THE SOLUTION IS ALWAYS UNDERSATURATED WITH LC
804       ELSE IF (YLO.LT.ZERO .AND. YHI.LT.ZERO) THEN
805          Z1 = ZHI
806          Z2 = ZHI
807          GOTO 40
809 ! *** { YLO, YHI } > 0.0 THE SOLUTION IS ALWAYS SUPERSATURATED WITH LC
811       ELSE IF (YLO.GT.ZERO .AND. YHI.GT.ZERO) THEN
812          Z1 = ZLO
813          Z2 = ZLO
814          GOTO 40
815       ELSE
816          CALL PUSHERR (0001, 'CALCB3A')    ! WARNING ERROR: NO SOLUTION
817          RETURN
818       ENDIF
820 ! *** PERFORM BISECTION ***********************************************
822 20    DO 30 I=1,MAXIT
823          Z3 = 0.5*(Z1+Z2)
824          Y3 = FUNCB3A (Z3, TLC, TNH42S4)
825          IF (SIGN(1.d0,Y1)*SIGN(1.d0,Y3) .LE. ZERO) THEN  ! (Y1*Y3 .LE. ZERO)
826             Y2    = Y3
827             Z2    = Z3
828          ELSE
829             Y1    = Y3
830             Z1    = Z3
831          ENDIF
832          IF (ABS(Z2-Z1) .LE. EPS*Z1) GOTO 40
833 30    CONTINUE
834       CALL PUSHERR (0002, 'CALCB3A')    ! WARNING ERROR: NO CONVERGENCE
836 ! *** CONVERGED ; RETURN ************************************************
838 40    ZK = 0.5*(Z1+Z2)
839       Y3 = FUNCB3A (ZK, TLC, TNH42S4)
841       RETURN
843 ! *** END OF SUBROUTINE CALCB3A ******************************************
845       END SUBROUTINE CALCB3A               
849 !=======================================================================
851 ! *** ISORROPIA CODE
852 ! *** FUNCTION FUNCB3A
853 ! *** CASE B3 ; SUBCASE 1
854 !     FUNCTION THAT SOLVES THE SYSTEM OF EQUATIONS FOR CASE B3
855 !     AND RETURNS THE VALUE OF THE ZEROED FUNCTION IN FUNCA3.
857 !=======================================================================
859       DOUBLE PRECISION FUNCTION FUNCB3A (ZK, Y, X)
860       USE ISRPIA
861       IMPLICIT DOUBLE PRECISION (A-H, O-Z)
862 !      INCLUDE 'ISRPIA.INC'
863       DOUBLE PRECISION KK
865 ! *** SOLVE EQUATIONS ; WITH ITERATIONS FOR ACTIVITY COEF. ************
867       FRST   = .TRUE.
868       CALAIN = .TRUE.
869       DO 20 I=1,NSWEEP
870          GRAT1 = XK1*WATER/GAMA(7)*(GAMA(8)/GAMA(7))**2.
871          DD    = SQRT( (ZK+GRAT1+Y)**2. + 4.0*Y*GRAT1)
872          KK    = 0.5*(-(ZK+GRAT1+Y) + DD )
874 ! *** SPECIATION & WATER CONTENT ***************************************
876          MOLAL (1) = KK                ! HI
877          MOLAL (5) = KK+ZK+Y           ! SO4I
878          MOLAL (6) = MAX (Y-KK, TINY)  ! HSO4I
879          MOLAL (3) = 3.0*Y+2*ZK        ! NH4I
880          CNH42S4   = X-ZK              ! Solid (NH4)2SO4
881          CALL CALCMR                   ! Water content
883 ! *** CALCULATE ACTIVITIES OR TERMINATE INTERNAL LOOP *****************
885          IF (FRST.AND.CALAOU .OR. .NOT.FRST.AND.CALAIN) THEN
886             CALL CALCACT
887          ELSE
888             GOTO 30
889          ENDIF
890 20    CONTINUE
892 ! *** CALCULATE OBJECTIVE FUNCTION ************************************
894 !CC30    FUNCB3A= ( SO4I*NH4I**2.0 )/( XK7*(WATER/GAMA(4))**3.0 )
895 30    FUNCB3A= MOLAL(5)*MOLAL(3)**2.0
896       FUNCB3A= FUNCB3A/(XK7*(WATER/GAMA(4))**3.0) - ONE
897       RETURN
899 ! *** END OF FUNCTION FUNCB3A ********************************************
901       END FUNCTION FUNCB3A           
905 !=======================================================================
907 ! *** ISORROPIA CODE
908 ! *** SUBROUTINE CALCB3B
909 ! *** CASE B3 ; SUBCASE 2
911 !     THE MAIN CHARACTERISTICS OF THIS REGIME ARE:
912 !     1. SULFATE RICH (1.0 < SULRAT < 2.0)
913 !     2. LIQUID PHASE ONLY IS POSSIBLE
915 !     SPECIATION CALCULATIONS IS BASED ON THE HSO4 <--> SO4 EQUILIBRIUM.
917 ! *** COPYRIGHT 1996-2006, UNIVERSITY OF MIAMI, CARNEGIE MELLON UNIVERSITY,
918 ! *** GEORGIA INSTITUTE OF TECHNOLOGY
919 ! *** WRITTEN BY ATHANASIOS NENES
920 ! *** UPDATED BY CHRISTOS FOUNTOUKIS
922 !=======================================================================
924       SUBROUTINE CALCB3B (Y, X)
925       USE ISRPIA
926       IMPLICIT DOUBLE PRECISION (A-H, O-Z)
927 !      INCLUDE 'ISRPIA.INC'
928       DOUBLE PRECISION KK
930       CALAOU = .FALSE.        ! Outer loop activity calculation flag
931       FRST   = .FALSE.
932       CALAIN = .TRUE.
934 ! *** SOLVE EQUATIONS ; WITH ITERATIONS FOR ACTIVITY COEF. ************
936       DO 20 I=1,NSWEEP
937          GRAT1 = XK1*WATER/GAMA(7)*(GAMA(8)/GAMA(7))**2.
938          DD    = SQRT( (GRAT1+Y)**2. + 4.0*(X+Y)*GRAT1)
939          KK    = 0.5*(-(GRAT1+Y) + DD )
941 ! *** SPECIATION & WATER CONTENT ***************************************
943          MOLAL (1) = KK                   ! HI
944          MOLAL (5) = Y+KK                 ! SO4I
945          MOLAL (6) = MAX (X+Y-KK, TINY)   ! HSO4I
946          MOLAL (3) = 3.0*Y+X              ! NH4I
947          CALL CALCMR                      ! Water content
949 ! *** CALCULATE ACTIVITIES OR TERMINATE INTERNAL LOOP *****************
951          IF (.NOT.CALAIN) GOTO 30
952          CALL CALCACT
953 20    CONTINUE
955 30    RETURN
957 ! *** END OF SUBROUTINE CALCB3B ******************************************
959       END SUBROUTINE CALCB3B       
960 !=======================================================================
962 ! *** ISORROPIA CODE
963 ! *** SUBROUTINE CALCB2
964 ! *** CASE B2
966 !     THE MAIN CHARACTERISTICS OF THIS REGIME ARE:
967 !     1. SULFATE RICH, NO FREE ACID (1.0 <= SULRAT < 2.0)
968 !     2. THERE IS BOTH A LIQUID & SOLID PHASE
969 !     3. SOLIDS POSSIBLE : LC, (NH4)2SO4
971 !     THERE ARE TWO POSSIBLE REGIMES HERE, DEPENDING ON THE SULFATE RATIO:
972 !     1. WHEN BOTH LC AND (NH4)2SO4 ARE POSSIBLE (SUBROUTINE CALCB2A)
973 !     2. WHEN ONLY LC IS POSSIBLE (SUBROUTINE CALCB2B).
975 ! *** COPYRIGHT 1996-2006, UNIVERSITY OF MIAMI, CARNEGIE MELLON UNIVERSITY,
976 ! *** GEORGIA INSTITUTE OF TECHNOLOGY
977 ! *** WRITTEN BY ATHANASIOS NENES
978 ! *** UPDATED BY CHRISTOS FOUNTOUKIS
980 !=======================================================================
982       SUBROUTINE CALCB2
983       USE ISRPIA
984       IMPLICIT DOUBLE PRECISION (A-H, O-Z)
985 !      INCLUDE 'ISRPIA.INC'
987 ! *** CALCULATE EQUIVALENT AMOUNT OF HSO4 AND SO4 ***********************
989       X = MAX(2*W(2)-W(3), TINY)   ! Equivalent NH4HSO4
990       Y = MAX(W(3)  -W(2), TINY)   ! Equivalent NH42SO4
992 ! *** CALCULATE SPECIES ACCORDING TO RELATIVE ABUNDANCE OF HSO4 *********
994       IF (X.LE.Y) THEN             ! LC is the MIN (x,y)
995          SCASE = 'B2 ; SUBCASE 1'
996          CALL CALCB2A (X,Y-X)      ! LC + (NH4)2SO4 POSSIBLE
997       ELSE
998          SCASE = 'B2 ; SUBCASE 2'
999          CALL CALCB2B (Y,X-Y)      ! LC ONLY POSSIBLE
1000       ENDIF
1002       RETURN
1004 ! *** END OF SUBROUTINE CALCB2 ******************************************
1006       END SUBROUTINE CALCB2
1010 !=======================================================================
1012 ! *** ISORROPIA CODE
1013 ! *** SUBROUTINE CALCB2
1014 ! *** CASE B2 ; SUBCASE A.
1016 !     THE MAIN CHARACTERISTICS OF THIS REGIME ARE:
1017 !     1. SULFATE RICH (1.0 < SULRAT < 2.0)
1018 !     2. SOLID PHASE ONLY POSSIBLE
1019 !     3. SOLIDS POSSIBLE: LC, (NH4)2SO4
1021 !     THERE ARE TWO POSSIBLE REGIMES HERE, DEPENDING ON RELATIVE HUMIDITY:
1022 !     1. WHEN RH >= MDRH ; LIQUID PHASE POSSIBLE (MDRH REGION)
1023 !     2. WHEN RH < MDRH  ; ONLY SOLID PHASE POSSIBLE
1025 !     FOR SOLID CALCULATIONS, A MATERIAL BALANCE BASED ON THE STOICHIMETRIC
1026 !     PROPORTION OF AMMONIUM AND SULFATE IS DONE TO CALCULATE THE AMOUNT
1027 !     OF LC AND (NH4)2SO4 IN THE SOLID PHASE.
1029 ! *** COPYRIGHT 1996-2006, UNIVERSITY OF MIAMI, CARNEGIE MELLON UNIVERSITY,
1030 ! *** GEORGIA INSTITUTE OF TECHNOLOGY
1031 ! *** WRITTEN BY ATHANASIOS NENES
1032 ! *** UPDATED BY CHRISTOS FOUNTOUKIS
1034 !=======================================================================
1036       SUBROUTINE CALCB2A (TLC, TNH42S4)
1037       USE ISRPIA
1038       IMPLICIT DOUBLE PRECISION (A-H, O-Z)
1039 !      INCLUDE 'ISRPIA.INC'
1041 ! *** REGIME DEPENDS UPON THE AMBIENT RELATIVE HUMIDITY *****************
1043       IF (RH.LT.DRMLCAS) THEN
1044          SCASE   = 'B2 ; SUBCASE A1'    ! SOLIDS POSSIBLE ONLY
1045          CLC     = TLC
1046          CNH42S4 = TNH42S4
1047          SCASE   = 'B2 ; SUBCASE A1'
1048       ELSE
1049          SCASE = 'B2 ; SUBCASE A2'
1050          CALL CALCB2A2 (TLC, TNH42S4)   ! LIQUID & SOLID PHASE POSSIBLE
1051          SCASE = 'B2 ; SUBCASE A2'
1052       ENDIF
1054       RETURN
1056 ! *** END OF SUBROUTINE CALCB2A *****************************************
1058       END SUBROUTINE CALCB2A               
1062 !=======================================================================
1064 ! *** ISORROPIA CODE
1065 ! *** SUBROUTINE CALCB2A2
1066 ! *** CASE B2 ; SUBCASE A2.
1068 !     THE MAIN CHARACTERISTICS OF THIS REGIME ARE:
1069 !     1. SULFATE RICH (1.0 < SULRAT < 2.0)
1070 !     2. SOLID PHASE ONLY POSSIBLE
1071 !     3. SOLIDS POSSIBLE: LC, (NH4)2SO4
1073 !     THIS IS THE CASE WHERE THE RELATIVE HUMIDITY IS IN THE MUTUAL
1074 !     DRH REGION. THE SOLUTION IS ASSUMED TO BE THE SUM OF TWO WEIGHTED
1075 !     SOLUTIONS ; THE SOLID PHASE ONLY (SUBROUTINE CALCB2A1) AND THE
1076 !     THE SOLID WITH LIQUID PHASE (SUBROUTINE CALCB3).
1078 ! *** COPYRIGHT 1996-2006, UNIVERSITY OF MIAMI, CARNEGIE MELLON UNIVERSITY,
1079 ! *** GEORGIA INSTITUTE OF TECHNOLOGY
1080 ! *** WRITTEN BY ATHANASIOS NENES
1081 ! *** UPDATED BY CHRISTOS FOUNTOUKIS
1083 !=======================================================================
1085       SUBROUTINE CALCB2A2 (TLC, TNH42S4)
1086       USE ISRPIA
1087       IMPLICIT DOUBLE PRECISION (A-H, O-Z)
1088 !      INCLUDE 'ISRPIA.INC'
1090 ! *** FIND WEIGHT FACTOR **********************************************
1092       IF (WFTYP.EQ.0) THEN
1093          WF = ZERO
1094       ELSEIF (WFTYP.EQ.1) THEN
1095          WF = 0.5D0
1096       ELSE
1097          WF = (DRLC-RH)/(DRLC-DRMLCAS)
1098       ENDIF
1099       ONEMWF  = ONE - WF
1101 ! *** FIND FIRST SECTION ; DRY ONE ************************************
1103       CLCO     = TLC                     ! FIRST (DRY) SOLUTION
1104       CNH42SO  = TNH42S4
1106 ! *** FIND SECOND SECTION ; DRY & LIQUID ******************************
1108       CLC     = ZERO
1109       CNH42S4 = ZERO
1110       CALL CALCB3                        ! SECOND (LIQUID) SOLUTION
1112 ! *** FIND SOLUTION AT MDRH BY WEIGHTING DRY & LIQUID SOLUTIONS.
1114       MOLAL(1)= ONEMWF*MOLAL(1)                                   ! H+
1115       MOLAL(3)= ONEMWF*(2.D0*(CNH42SO-CNH42S4) + 3.D0*(CLCO-CLC)) ! NH4+
1116       MOLAL(5)= ONEMWF*(CNH42SO-CNH42S4 + CLCO-CLC)               ! SO4--
1117       MOLAL(6)= ONEMWF*(CLCO-CLC)                                 ! HSO4-
1119       WATER   = ONEMWF*WATER
1121       CLC     = WF*CLCO    + ONEMWF*CLC
1122       CNH42S4 = WF*CNH42SO + ONEMWF*CNH42S4
1124       RETURN
1126 ! *** END OF SUBROUTINE CALCB2A2 ****************************************
1128       END SUBROUTINE CALCB2A2               
1132 !=======================================================================
1134 ! *** ISORROPIA CODE
1135 ! *** SUBROUTINE CALCB2
1136 ! *** CASE B2 ; SUBCASE B
1138 !     THE MAIN CHARACTERISTICS OF THIS REGIME ARE:
1139 !     1. SULFATE RICH (1.0 < SULRAT < 2.0)
1140 !     2. BOTH LIQUID & SOLID PHASE IS POSSIBLE
1141 !     3. SOLIDS POSSIBLE: LC
1143 !     FOR CALCULATIONS, A BISECTION IS PERFORMED TOWARDS ZETA, THE
1144 !     AMOUNT OF SOLID LC DISSOLVED IN THE LIQUID PHASE.
1145 !     FOR EACH ESTIMATION OF ZETA, FUNCTION FUNCB2A CALCULATES THE
1146 !     AMOUNT OF H+ PRODUCED (BASED ON THE HSO4, SO4 RELEASED INTO THE
1147 !     SOLUTION). THE SOLUBILITY PRODUCT OF LC IS USED AS THE OBJECTIVE
1148 !     FUNCTION.
1150 ! *** COPYRIGHT 1996-2006, UNIVERSITY OF MIAMI, CARNEGIE MELLON UNIVERSITY,
1151 ! *** GEORGIA INSTITUTE OF TECHNOLOGY
1152 ! *** WRITTEN BY ATHANASIOS NENES
1153 ! *** UPDATED BY CHRISTOS FOUNTOUKIS
1155 !=======================================================================
1157       SUBROUTINE CALCB2B (TLC,TNH4HS4)
1158       USE ISRPIA
1159       IMPLICIT DOUBLE PRECISION (A-H, O-Z)
1160 !      INCLUDE 'ISRPIA.INC'
1162       CALAOU = .TRUE.       ! Outer loop activity calculation flag
1163       ZLO    = ZERO
1164       ZHI    = TLC          ! High limit: all of it in liquid phase
1166 ! *** INITIAL VALUES FOR BISECTION **************************************
1168       X1 = ZHI
1169       Y1 = FUNCB2B (X1,TNH4HS4,TLC)
1170       IF (ABS(Y1).LE.EPS) RETURN
1171       YHI= Y1                        ! Save Y-value at Hi position
1173 ! *** ROOT TRACKING ; FOR THE RANGE OF HI AND LO ************************
1175       DX = (ZHI-ZLO)/NDIV
1176       DO 10 I=1,NDIV
1177          X2 = X1-DX
1178          Y2 = FUNCB2B (X2,TNH4HS4,TLC)
1179          IF (SIGN(1.d0,Y1)*SIGN(1.d0,Y2).LT.ZERO) GOTO 20  ! (Y1*Y2.LT.ZERO)
1180          X1 = X2
1181          Y1 = Y2
1182 10    CONTINUE
1184 ! *** NO SUBDIVISION WITH SOLUTION FOUND
1186       YLO= Y1                      ! Save Y-value at LO position
1187       IF (ABS(Y2) .LT. EPS) THEN   ! X2 IS A SOLUTION
1188          RETURN
1190 ! *** { YLO, YHI } < 0.0 THE SOLUTION IS ALWAYS UNDERSATURATED WITH LC
1192       ELSE IF (YLO.LT.ZERO .AND. YHI.LT.ZERO) THEN
1193          X1 = ZHI
1194          X2 = ZHI
1195          GOTO 40
1197 ! *** { YLO, YHI } > 0.0 THE SOLUTION IS ALWAYS SUPERSATURATED WITH LC
1199       ELSE IF (YLO.GT.ZERO .AND. YHI.GT.ZERO) THEN
1200          X1 = ZLO
1201          X2 = ZLO
1202          GOTO 40
1203       ELSE
1204          CALL PUSHERR (0001, 'CALCB2B')    ! WARNING ERROR: NO SOLUTION
1205          RETURN
1206       ENDIF
1208 ! *** PERFORM BISECTION *************************************************
1210 20    DO 30 I=1,MAXIT
1211          X3 = 0.5*(X1+X2)
1212          Y3 = FUNCB2B (X3,TNH4HS4,TLC)
1213          IF (SIGN(1.d0,Y1)*SIGN(1.d0,Y3) .LE. ZERO) THEN  ! (Y1*Y3 .LE. ZERO)
1214             Y2    = Y3
1215             X2    = X3
1216          ELSE
1217             Y1    = Y3
1218             X1    = X3
1219          ENDIF
1220          IF (ABS(X2-X1) .LE. EPS*X1) GOTO 40
1221 30    CONTINUE
1222       CALL PUSHERR (0002, 'CALCB2B')    ! WARNING ERROR: NO CONVERGENCE
1224 ! *** CONVERGED ; RETURN ************************************************
1226 40    X3 = 0.5*(X1+X2)
1227       Y3 = FUNCB2B (X3,TNH4HS4,TLC)
1229       RETURN
1231 ! *** END OF SUBROUTINE CALCB2B *****************************************
1233       END SUBROUTINE CALCB2B              
1237 !=======================================================================
1239 ! *** ISORROPIA CODE
1240 ! *** FUNCTION FUNCB2B
1241 ! *** CASE B2 ;
1242 !     FUNCTION THAT SOLVES THE SYSTEM OF EQUATIONS FOR CASE B2 ; SUBCASE 2
1243 !     AND RETURNS THE VALUE OF THE ZEROED FUNCTION IN FUNCB2B.
1245 !=======================================================================
1247       DOUBLE PRECISION FUNCTION FUNCB2B (X,TNH4HS4,TLC)
1248       USE ISRPIA
1249       IMPLICIT DOUBLE PRECISION (A-H, O-Z)
1250 !      INCLUDE 'ISRPIA.INC'
1252 ! *** SOLVE EQUATIONS **************************************************
1254       FRST   = .TRUE.
1255       CALAIN = .TRUE.
1256       DO 20 I=1,NSWEEP
1257          GRAT2 = XK1*WATER*(GAMA(8)/GAMA(7))**2./GAMA(7)
1258          PARM  = X+GRAT2
1259          DELTA = PARM*PARM + 4.0*(X+TNH4HS4)*GRAT2 ! Diakrinousa
1260          OMEGA = 0.5*(-PARM + SQRT(DELTA))         ! Thetiki riza (ie:H+>0)
1262 ! *** SPECIATION & WATER CONTENT ***************************************
1264          MOLAL (1) = OMEGA                         ! HI
1265          MOLAL (3) = 3.0*X+TNH4HS4                 ! NH4I
1266          MOLAL (5) = X+OMEGA                       ! SO4I
1267          MOLAL (6) = MAX (X+TNH4HS4-OMEGA, TINY)   ! HSO4I
1268          CLC       = MAX(TLC-X,ZERO)               ! Solid LC
1269          CALL CALCMR                               ! Water content
1271 ! *** CALCULATE ACTIVITIES OR TERMINATE INTERNAL LOOP ******************
1273          IF (FRST.AND.CALAOU .OR. .NOT.FRST.AND.CALAIN) THEN
1274             CALL CALCACT
1275          ELSE
1276             GOTO 30
1277          ENDIF
1278 20    CONTINUE
1280 ! *** CALCULATE OBJECTIVE FUNCTION **************************************
1282 !CC30    FUNCB2B= ( NH4I**3.*SO4I*HSO4I )/( XK13*(WATER/GAMA(13))**5. )
1283 30    FUNCB2B= (MOLAL(3)**3.)*MOLAL(5)*MOLAL(6)
1284       FUNCB2B= FUNCB2B/(XK13*(WATER/GAMA(13))**5.) - ONE
1285       RETURN
1287 ! *** END OF FUNCTION FUNCB2B *******************************************
1289       END FUNCTION FUNCB2B                
1292 !=======================================================================
1294 ! *** ISORROPIA CODE
1295 ! *** SUBROUTINE CALCB1
1296 ! *** CASE B1
1298 !     THE MAIN CHARACTERISTICS OF THIS REGIME ARE:
1299 !     1. SULFATE RICH, NO FREE ACID (1.0 <= SULRAT < 2.0)
1300 !     2. THERE IS BOTH A LIQUID & SOLID PHASE
1301 !     3. SOLIDS POSSIBLE : LC, (NH4)2SO4, NH4HSO4
1303 !     THERE ARE TWO POSSIBLE REGIMES HERE, DEPENDING ON RELATIVE HUMIDITY:
1304 !     1. WHEN RH >= MDRH ; LIQUID PHASE POSSIBLE (MDRH REGION)
1305 !     2. WHEN RH < MDRH  ; ONLY SOLID PHASE POSSIBLE (SUBROUTINE CALCB1A)
1307 ! *** COPYRIGHT 1996-2006, UNIVERSITY OF MIAMI, CARNEGIE MELLON UNIVERSITY,
1308 ! *** GEORGIA INSTITUTE OF TECHNOLOGY
1309 ! *** WRITTEN BY ATHANASIOS NENES
1310 ! *** UPDATED BY CHRISTOS FOUNTOUKIS
1312 !=======================================================================
1314       SUBROUTINE CALCB1
1315       USE ISRPIA
1316       IMPLICIT DOUBLE PRECISION (A-H, O-Z)
1317 !      INCLUDE 'ISRPIA.INC'
1319 ! *** REGIME DEPENDS UPON THE AMBIENT RELATIVE HUMIDITY *****************
1321       IF (RH.LT.DRMLCAB) THEN
1322          SCASE = 'B1 ; SUBCASE 1'
1323          CALL CALCB1A              ! SOLID PHASE ONLY POSSIBLE
1324          SCASE = 'B1 ; SUBCASE 1'
1325       ELSE
1326          SCASE = 'B1 ; SUBCASE 2'
1327          CALL CALCB1B              ! LIQUID & SOLID PHASE POSSIBLE
1328          SCASE = 'B1 ; SUBCASE 2'
1329       ENDIF
1331       RETURN
1333 ! *** END OF SUBROUTINE CALCB1 ******************************************
1335       END SUBROUTINE CALCB1
1339 !=======================================================================
1341 ! *** ISORROPIA CODE
1342 ! *** SUBROUTINE CALCB1A
1343 ! *** CASE B1 ; SUBCASE 1
1345 !     THE MAIN CHARACTERISTICS OF THIS REGIME ARE:
1346 !     1. SULFATE RICH
1347 !     2. THERE IS NO LIQUID PHASE
1348 !     3. SOLIDS POSSIBLE: LC, { (NH4)2SO4  XOR  NH4HSO4 } (ONE OF TWO
1349 !                         BUT NOT BOTH)
1351 !     A SIMPLE MATERIAL BALANCE IS PERFORMED, AND THE AMOUNT OF LC
1352 !     IS CALCULATED FROM THE (NH4)2SO4 AND NH4HSO4 WHICH IS LEAST
1353 !     ABUNDANT (STOICHIMETRICALLY). THE REMAINING EXCESS OF SALT
1354 !     IS MIXED WITH THE LC.
1356 ! *** COPYRIGHT 1996-2006, UNIVERSITY OF MIAMI, CARNEGIE MELLON UNIVERSITY,
1357 ! *** GEORGIA INSTITUTE OF TECHNOLOGY
1358 ! *** WRITTEN BY ATHANASIOS NENES
1359 ! *** UPDATED BY CHRISTOS FOUNTOUKIS
1361 !=======================================================================
1363       SUBROUTINE CALCB1A
1364       USE ISRPIA
1365       IMPLICIT DOUBLE PRECISION (A-H, O-Z)
1366 !      INCLUDE 'ISRPIA.INC'
1368 ! *** SETUP PARAMETERS ************************************************
1370       X = 2*W(2)-W(3)       ! Equivalent NH4HSO4
1371       Y = W(3)-W(2)         ! Equivalent (NH4)2SO4
1373 ! *** CALCULATE COMPOSITION *******************************************
1375       IF (X.LE.Y) THEN      ! LC is the MIN (x,y)
1376          CLC     = X        ! NH4HSO4 >= (NH4)2S04
1377          CNH4HS4 = ZERO
1378          CNH42S4 = Y-X
1379       ELSE
1380          CLC     = Y        ! NH4HSO4 <  (NH4)2S04
1381          CNH4HS4 = X-Y
1382          CNH42S4 = ZERO
1383       ENDIF
1384       RETURN
1386 ! *** END OF SUBROUTINE CALCB1 ******************************************
1388       END SUBROUTINE CALCB1A
1393 !=======================================================================
1395 ! *** ISORROPIA CODE
1396 ! *** SUBROUTINE CALCB1B
1397 ! *** CASE B1 ; SUBCASE 2
1399 !     THE MAIN CHARACTERISTICS OF THIS REGIME ARE:
1400 !     1. SULFATE RICH
1401 !     2. THERE IS BOTH A LIQUID & SOLID PHASE
1402 !     3. SOLIDS POSSIBLE: LC, { (NH4)2SO4  XOR  NH4HSO4 } (ONE OF TWO
1403 !                         BUT NOT BOTH)
1405 !     THIS IS THE CASE WHERE THE RELATIVE HUMIDITY IS IN THE MUTUAL
1406 !     DRH REGION. THE SOLUTION IS ASSUMED TO BE THE SUM OF TWO WEIGHTED
1407 !     SOLUTIONS ; THE SOLID PHASE ONLY (SUBROUTINE CALCB1A) AND THE
1408 !     THE SOLID WITH LIQUID PHASE (SUBROUTINE CALCB2).
1410 ! *** COPYRIGHT 1996-2006, UNIVERSITY OF MIAMI, CARNEGIE MELLON UNIVERSITY,
1411 ! *** GEORGIA INSTITUTE OF TECHNOLOGY
1412 ! *** WRITTEN BY ATHANASIOS NENES
1413 ! *** UPDATED BY CHRISTOS FOUNTOUKIS
1415 !=======================================================================
1417       SUBROUTINE CALCB1B
1418       USE ISRPIA
1419       IMPLICIT DOUBLE PRECISION (A-H, O-Z)
1420 !      INCLUDE 'ISRPIA.INC'
1422 ! *** FIND WEIGHT FACTOR **********************************************
1424       IF (WFTYP.EQ.0) THEN
1425          WF = ZERO
1426       ELSEIF (WFTYP.EQ.1) THEN
1427          WF = 0.5D0
1428       ELSE
1429          WF = (DRNH4HS4-RH)/(DRNH4HS4-DRMLCAB)
1430       ENDIF
1431       ONEMWF  = ONE - WF
1433 ! *** FIND FIRST SECTION ; DRY ONE ************************************
1435       CALL CALCB1A
1436       CLCO     = CLC               ! FIRST (DRY) SOLUTION
1437       CNH42SO  = CNH42S4
1438       CNH4HSO  = CNH4HS4
1440 ! *** FIND SECOND SECTION ; DRY & LIQUID ******************************
1442       CLC     = ZERO
1443       CNH42S4 = ZERO
1444       CNH4HS4 = ZERO
1445       CALL CALCB2                  ! SECOND (LIQUID) SOLUTION
1447 ! *** FIND SOLUTION AT MDRH BY WEIGHTING DRY & LIQUID SOLUTIONS.
1449       MOLAL(1)= ONEMWF*MOLAL(1)                                   ! H+
1450       MOLAL(3)= ONEMWF*(2.D0*(CNH42SO-CNH42S4) + (CNH4HSO-CNH4HS4)   &
1451                       + 3.D0*(CLCO-CLC))                          ! NH4+
1452       MOLAL(5)= ONEMWF*(CNH42SO-CNH42S4 + CLCO-CLC)               ! SO4--
1453       MOLAL(6)= ONEMWF*(CNH4HSO-CNH4HS4 + CLCO-CLC)               ! HSO4-
1455       WATER   = ONEMWF*WATER
1457       CLC     = WF*CLCO    + ONEMWF*CLC
1458       CNH42S4 = WF*CNH42SO + ONEMWF*CNH42S4
1459       CNH4HS4 = WF*CNH4HSO + ONEMWF*CNH4HS4
1461       RETURN
1463 ! *** END OF SUBROUTINE CALCB1B *****************************************
1465       END SUBROUTINE CALCB1B
1468 !=======================================================================
1470 ! *** ISORROPIA CODE
1471 ! *** SUBROUTINE CALCC2
1472 ! *** CASE C2
1474 !     THE MAIN CHARACTERISTICS OF THIS REGIME ARE:
1475 !     1. SULFATE RICH, FREE ACID (SULRAT < 1.0)
1476 !     2. THERE IS ONLY A LIQUID PHASE
1478 ! *** COPYRIGHT 1996-2006, UNIVERSITY OF MIAMI, CARNEGIE MELLON UNIVERSITY,
1479 ! *** GEORGIA INSTITUTE OF TECHNOLOGY
1480 ! *** WRITTEN BY ATHANASIOS NENES
1481 ! *** UPDATED BY CHRISTOS FOUNTOUKIS
1483 !=======================================================================
1485       SUBROUTINE CALCC2
1486       USE ISRPIA
1487       IMPLICIT DOUBLE PRECISION (A-H, O-Z)
1488 !      INCLUDE 'ISRPIA.INC'
1489       DOUBLE PRECISION LAMDA, KAPA
1491       CALAOU =.TRUE.         ! Outer loop activity calculation flag
1492       FRST   =.TRUE.
1493       CALAIN =.TRUE.
1495 ! *** SOLVE EQUATIONS **************************************************
1497       LAMDA  = W(3)           ! NH4HSO4 INITIALLY IN SOLUTION
1498       PSI    = W(2)-W(3)      ! H2SO4 IN SOLUTION
1499       DO 20 I=1,NSWEEP
1500          PARM  = WATER*XK1/GAMA(7)*(GAMA(8)/GAMA(7))**2.
1501          BB    = PSI+PARM
1502          CC    =-PARM*(LAMDA+PSI)
1503          KAPA  = 0.5*(-BB+SQRT(BB*BB-4.0*CC))
1505 ! *** SPECIATION & WATER CONTENT ***************************************
1507          MOLAL(1) = PSI+KAPA                               ! HI
1508          MOLAL(3) = LAMDA                                  ! NH4I
1509          MOLAL(5) = KAPA                                   ! SO4I
1510          MOLAL(6) = MAX(LAMDA+PSI-KAPA, TINY)              ! HSO4I
1511          CH2SO4   = MAX(MOLAL(5)+MOLAL(6)-MOLAL(3), ZERO)  ! Free H2SO4
1512          CALL CALCMR                                       ! Water content
1514 ! *** CALCULATE ACTIVITIES OR TERMINATE INTERNAL LOOP *****************
1516          IF (.NOT.CALAIN) GOTO 30
1517          CALL CALCACT
1518 20    CONTINUE
1520 30    RETURN
1522 ! *** END OF SUBROUTINE CALCC2 *****************************************
1524       END SUBROUTINE CALCC2
1528 !=======================================================================
1530 ! *** ISORROPIA CODE
1531 ! *** SUBROUTINE CALCC1
1532 ! *** CASE C1
1534 !     THE MAIN CHARACTERISTICS OF THIS REGIME ARE:
1535 !     1. SULFATE RICH, FREE ACID (SULRAT < 1.0)
1536 !     2. THERE IS BOTH A LIQUID & SOLID PHASE
1537 !     3. SOLIDS POSSIBLE: NH4HSO4
1539 ! *** COPYRIGHT 1996-2006, UNIVERSITY OF MIAMI, CARNEGIE MELLON UNIVERSITY,
1540 ! *** GEORGIA INSTITUTE OF TECHNOLOGY
1541 ! *** WRITTEN BY ATHANASIOS NENES
1542 ! *** UPDATED BY CHRISTOS FOUNTOUKIS
1544 !=======================================================================
1546       SUBROUTINE CALCC1
1547       USE ISRPIA
1548       IMPLICIT DOUBLE PRECISION (A-H, O-Z)
1549 !      INCLUDE 'ISRPIA.INC'
1550       DOUBLE PRECISION KLO, KHI
1552       CALAOU = .TRUE.    ! Outer loop activity calculation flag
1553       KLO    = TINY
1554       KHI    = W(3)
1556 ! *** INITIAL VALUES FOR BISECTION *************************************
1558       X1 = KLO
1559       Y1 = FUNCC1 (X1)
1560       IF (ABS(Y1).LE.EPS) GOTO 50
1561       YLO= Y1
1563 ! *** ROOT TRACKING ; FOR THE RANGE OF HI AND LO ***********************
1565       DX = (KHI-KLO)/REAL(NDIV)
1566       DO 10 I=1,NDIV
1567          X2 = X1+DX
1568          Y2 = FUNCC1 (X2)
1569          IF (SIGN(1.d0,Y1)*SIGN(1.d0,Y2) .LT. ZERO) GOTO 20 ! (Y1*Y2 .LT. ZERO)
1570          X1 = X2
1571          Y1 = Y2
1572 10    CONTINUE
1574 ! *** NO SUBDIVISION WITH SOLUTION FOUND
1576       YHI= Y2                 ! Save Y-value at HI position
1577       IF (ABS(Y2) .LT. EPS) THEN   ! X2 IS A SOLUTION
1578          GOTO 50
1580 ! *** { YLO, YHI } < 0.0  SOLUTION IS ALWAYS UNDERSATURATED WITH NH4HS04
1582       ELSE IF (YLO.LT.ZERO .AND. YHI.LT.ZERO) THEN
1583          GOTO 50
1585 ! *** { YLO, YHI } > 0.0 SOLUTION IS ALWAYS SUPERSATURATED WITH NH4HS04
1587       ELSE IF (YLO.GT.ZERO .AND. YHI.GT.ZERO) THEN
1588          X1 = KLO
1589          X2 = KLO
1590          GOTO 40
1591       ELSE
1592          CALL PUSHERR (0001, 'CALCC1')    ! WARNING ERROR: NO SOLUTION
1593          GOTO 50
1594       ENDIF
1596 ! *** PERFORM BISECTION OF DISSOLVED NH4HSO4 **************************
1598 20    DO 30 I=1,MAXIT
1599          X3 = 0.5*(X1+X2)
1600          Y3 = FUNCC1 (X3)
1601          IF (SIGN(1.d0,Y1)*SIGN(1.d0,Y3) .LE. ZERO) THEN  ! (Y1*Y3 .LE. ZERO)
1602             Y2    = Y3
1603             X2    = X3
1604          ELSE
1605             Y1    = Y3
1606             X1    = X3
1607          ENDIF
1608          IF (ABS(X2-X1) .LE. EPS*X1) GOTO 40
1609 30    CONTINUE
1610       CALL PUSHERR (0002, 'CALCC1')    ! WARNING ERROR: NO CONVERGENCE
1612 ! *** CONVERGED ; RETURN ***********************************************
1614 40    X3 = 0.5*(X1+X2)
1615       Y3 = FUNCC1 (X3)
1617 50    RETURN
1619 ! *** END OF SUBROUTINE CALCC1 *****************************************
1621       END SUBROUTINE CALCC1
1625 !=======================================================================
1627 ! *** ISORROPIA CODE
1628 ! *** FUNCTION FUNCC1
1629 ! *** CASE C1 ;
1630 !     FUNCTION THAT SOLVES THE SYSTEM OF EQUATIONS FOR CASE C1
1631 !     AND RETURNS THE VALUE OF THE ZEROED FUNCTION IN FUNCC1.
1633 ! *** COPYRIGHT 1996-2006, UNIVERSITY OF MIAMI, CARNEGIE MELLON UNIVERSITY,
1634 ! *** GEORGIA INSTITUTE OF TECHNOLOGY
1635 ! *** WRITTEN BY ATHANASIOS NENES
1636 ! *** UPDATED BY CHRISTOS FOUNTOUKIS
1638 !=======================================================================
1640       DOUBLE PRECISION FUNCTION FUNCC1 (KAPA)
1641       USE ISRPIA
1642       IMPLICIT DOUBLE PRECISION (A-H, O-Z)
1643 !      INCLUDE 'ISRPIA.INC'
1644       DOUBLE PRECISION KAPA, LAMDA
1646 ! *** SOLVE EQUATIONS **************************************************
1648       FRST   = .TRUE.
1649       CALAIN = .TRUE.
1651       PSI = W(2)-W(3)
1652       DO 20 I=1,NSWEEP
1653          PAR1  = XK1*WATER/GAMA(7)*(GAMA(8)/GAMA(7))**2.0
1654          PAR2  = XK12*(WATER/GAMA(9))**2.0
1655          BB    = PSI + PAR1
1656          CC    =-PAR1*(PSI+KAPA)
1657          LAMDA = 0.5*(-BB+SQRT(BB*BB-4*CC))
1659 ! *** SAVE CONCENTRATIONS IN MOLAL ARRAY *******************************
1661          MOLAL(1) = PSI+LAMDA                    ! HI
1662          MOLAL(3) = KAPA                         ! NH4I
1663          MOLAL(5) = LAMDA                        ! SO4I
1664          MOLAL(6) = MAX (ZERO, PSI+KAPA-LAMDA)   ! HSO4I
1665          CNH4HS4  = MAX(W(3)-MOLAL(3), ZERO)     ! Solid NH4HSO4
1666          CH2SO4   = MAX(PSI, ZERO)               ! Free H2SO4
1667          CALL CALCMR                             ! Water content
1669 ! *** CALCULATE ACTIVITIES OR TERMINATE INTERNAL LOOP *****************
1671          IF (FRST.AND.CALAOU .OR. .NOT.FRST.AND.CALAIN) THEN
1672             CALL CALCACT
1673          ELSE
1674             GOTO 30
1675          ENDIF
1676 20    CONTINUE
1678 ! *** CALCULATE ZERO FUNCTION *******************************************
1680 !CC30    FUNCC1= (NH4I*HSO4I/PAR2) - ONE
1681 30    FUNCC1= (MOLAL(3)*MOLAL(6)/PAR2) - ONE
1682       RETURN
1684 ! *** END OF FUNCTION FUNCC1 ********************************************
1686       END FUNCTION FUNCC1       
1688 !=======================================================================
1690 ! *** ISORROPIA CODE
1691 ! *** SUBROUTINE CALCD3
1692 ! *** CASE D3
1694 !     THE MAIN CHARACTERISTICS OF THIS REGIME ARE:
1695 !     1. SULFATE POOR (SULRAT > 2.0)
1696 !     2. THERE IS OLNY A LIQUID PHASE
1698 ! *** COPYRIGHT 1996-2006, UNIVERSITY OF MIAMI, CARNEGIE MELLON UNIVERSITY,
1699 ! *** GEORGIA INSTITUTE OF TECHNOLOGY
1700 ! *** WRITTEN BY ATHANASIOS NENES
1701 ! *** UPDATED BY CHRISTOS FOUNTOUKIS
1703 !=======================================================================
1705       SUBROUTINE CALCD3
1706       USE ISRPIA
1707       USE SOLUT
1708       IMPLICIT DOUBLE PRECISION (A-H, O-Z)
1709 !      INCLUDE 'ISRPIA.INC'
1710 !      COMMON /SOLUT/ CHI1, CHI2, CHI3, CHI4, CHI5, CHI6, CHI7, CHI8,   &
1711 !                     PSI1, PSI2, PSI3, PSI4, PSI5, PSI6, PSI7, PSI8,   &
1712 !                     A1,   A2,   A3,   A4,   A5,   A6,   A7,   A8
1714 ! *** FIND DRY COMPOSITION **********************************************
1716       CALL CALCD1A
1718 ! *** SETUP PARAMETERS ************************************************
1720       CHI1 = CNH4NO3               ! Save from CALCD1 run
1721       CHI2 = CNH42S4
1722       CHI3 = GHNO3
1723       CHI4 = GNH3
1725       PSI1 = CNH4NO3               ! ASSIGN INITIAL PSI's
1726       PSI2 = CHI2
1727       PSI3 = ZERO
1728       PSI4 = ZERO
1730       MOLAL(5) = ZERO
1731       MOLAL(6) = ZERO
1732       MOLAL(3) = PSI1
1733       MOLAL(7) = PSI1
1734       CALL CALCMR                  ! Initial water
1736       CALAOU = .TRUE.              ! Outer loop activity calculation flag
1737       PSI4LO = TINY                ! Low  limit
1738       PSI4HI = CHI4                ! High limit
1740 ! *** INITIAL VALUES FOR BISECTION ************************************
1742 60    X1 = PSI4LO
1743       Y1 = FUNCD3 (X1)
1744       IF (ABS(Y1).LE.EPS) RETURN
1745       YLO= Y1                 ! Save Y-value at HI position
1747 ! *** ROOT TRACKING ; FOR THE RANGE OF HI AND LO **********************
1749       DX = (PSI4HI-PSI4LO)/REAL(NDIV)
1750       DO 10 I=1,NDIV
1751          X2 = X1+DX
1752          Y2 = FUNCD3 (X2)
1753          IF (SIGN(1.d0,Y1)*SIGN(1.d0,Y2).LT.ZERO) GOTO 20  ! (Y1*Y2.LT.ZERO)
1754          X1 = X2
1755          Y1 = Y2
1756 10    CONTINUE
1758 ! *** NO SUBDIVISION WITH SOLUTION FOUND
1760       YHI= Y1                      ! Save Y-value at Hi position
1761       IF (ABS(Y2) .LT. EPS) THEN   ! X2 IS A SOLUTION
1762          RETURN
1764 ! *** { YLO, YHI } < 0.0 THE SOLUTION IS ALWAYS UNDERSATURATED WITH NH3
1765 ! Physically I dont know when this might happen, but I have put this
1766 ! branch in for completeness. I assume there is no solution; all NO3 goes to the
1767 ! gas phase.
1769       ELSE IF (YLO.LT.ZERO .AND. YHI.LT.ZERO) THEN
1770          P4 = TINY ! PSI4LO ! CHI4
1771          YY = FUNCD3(P4)
1772          GOTO 50
1774 ! *** { YLO, YHI } > 0.0 THE SOLUTION IS ALWAYS SUPERSATURATED WITH NH3
1775 ! This happens when Sul.Rat. = 2.0, so some NH4+ from sulfate evaporates
1776 ! and goes to the gas phase ; so I redefine the LO and HI limits of PSI4
1777 ! and proceed again with root tracking.
1779       ELSE IF (YLO.GT.ZERO .AND. YHI.GT.ZERO) THEN
1780          PSI4HI = PSI4LO
1781          PSI4LO = PSI4LO - 0.1*(PSI1+PSI2) ! No solution; some NH3 evaporates
1782          IF (PSI4LO.LT.-(PSI1+PSI2)) THEN
1783             CALL PUSHERR (0001, 'CALCD3')  ! WARNING ERROR: NO SOLUTION
1784             RETURN
1785          ELSE
1786             MOLAL(5) = ZERO
1787             MOLAL(6) = ZERO
1788             MOLAL(3) = PSI1
1789             MOLAL(7) = PSI1
1790             CALL CALCMR                  ! Initial water
1791             GOTO 60                        ! Redo root tracking
1792          ENDIF
1793       ENDIF
1795 ! *** PERFORM BISECTION ***********************************************
1797 20    DO 30 I=1,MAXIT
1798          X3 = 0.5*(X1+X2)
1799          Y3 = FUNCD3 (X3)
1800          IF (SIGN(1.d0,Y1)*SIGN(1.d0,Y3) .LE. ZERO) THEN  ! (Y1*Y3 .LE. ZERO)
1801             Y2    = Y3
1802             X2    = X3
1803          ELSE
1804             Y1    = Y3
1805             X1    = X3
1806          ENDIF
1807          IF (ABS(X2-X1) .LE. EPS*ABS(X1)) GOTO 40
1808 30    CONTINUE
1809       CALL PUSHERR (0002, 'CALCD3')    ! WARNING ERROR: NO CONVERGENCE
1811 ! *** CONVERGED ; RETURN **********************************************
1813 40    X3 = 0.5*(X1+X2)
1814       Y3 = FUNCD3 (X3)
1816 ! *** CALCULATE HSO4 SPECIATION AND RETURN *******************************
1818 50    CONTINUE
1819       IF (MOLAL(1).GT.TINY) THEN
1820          CALL CALCHS4 (MOLAL(1), MOLAL(5), ZERO, DELTA)
1821          MOLAL(1) = MOLAL(1) - DELTA                     ! H+   EFFECT
1822          MOLAL(5) = MOLAL(5) - DELTA                     ! SO4  EFFECT
1823          MOLAL(6) = DELTA                                ! HSO4 EFFECT
1824       ENDIF
1825       RETURN
1827 ! *** END OF SUBROUTINE CALCD3 ******************************************
1829       END SUBROUTINE CALCD3
1833 !=======================================================================
1835 ! *** ISORROPIA CODE
1836 ! *** FUNCTION FUNCD3
1837 ! *** CASE D3
1838 !     FUNCTION THAT SOLVES THE SYSTEM OF EQUATIONS FOR CASE D3 ;
1839 !     AND RETURNS THE VALUE OF THE ZEROED FUNCTION IN FUNCD3.
1841 !=======================================================================
1843       DOUBLE PRECISION FUNCTION FUNCD3 (P4)
1844       USE ISRPIA
1845       USE SOLUT
1846       IMPLICIT DOUBLE PRECISION (A-H, O-Z)
1847 !      INCLUDE 'ISRPIA.INC'
1848 !      COMMON /SOLUT/ CHI1, CHI2, CHI3, CHI4, CHI5, CHI6, CHI7, CHI8,   &
1849 !                     PSI1, PSI2, PSI3, PSI4, PSI5, PSI6, PSI7, PSI8,   &
1850 !                     A1,   A2,   A3,   A4,   A5,   A6,   A7,   A8
1852 ! *** SETUP PARAMETERS ************************************************
1854       FRST   = .TRUE.
1855       CALAIN = .TRUE.
1856       PSI4   = P4
1858 ! *** SOLVE EQUATIONS ; WITH ITERATIONS FOR ACTIVITY COEF. ************
1860       DO 10 I=1,NSWEEP
1861          A2   = XK7*(WATER/GAMA(4))**3.0
1862          A3   = XK4*R*TEMP*(WATER/GAMA(10))**2.0
1863          A4   = (XK2/XKW)*R*TEMP*(GAMA(10)/GAMA(5))**2.0
1864          A7   = XKW *RH*WATER*WATER
1866          PSI3 = A3*A4*CHI3*(CHI4-PSI4) - PSI1*(2.D0*PSI2+PSI1+PSI4)
1867          PSI3 = PSI3/(A3*A4*(CHI4-PSI4) + 2.D0*PSI2+PSI1+PSI4)
1868          PSI3 = MIN(MAX(PSI3, ZERO), CHI3)
1870          BB   = PSI4 - PSI3
1871 !CCOLD         AHI  = 0.5*(-BB + SQRT(BB*BB + 4.d0*A7)) ! This is correct also
1872 !CC         AHI  =2.0*A7/(BB+SQRT(BB*BB + 4.d0*A7)) ! Avoid overflow when HI->0
1873          DENM = BB+SQRT(BB*BB + 4.d0*A7)
1874          IF (DENM.LE.TINY) THEN       ! Avoid overflow when HI->0
1875             ABB  = ABS(BB)
1876             DENM = (BB+ABB) + 2.0*A7/ABB ! Taylor expansion of SQRT
1877          ENDIF
1878          AHI = 2.0*A7/DENM
1880 ! *** SPECIATION & WATER CONTENT ***************************************
1882          MOLAL (1) = AHI                             ! HI
1883          MOLAL (3) = PSI1 + PSI4 + 2.D0*PSI2         ! NH4I
1884          MOLAL (5) = PSI2                            ! SO4I
1885          MOLAL (6) = ZERO                            ! HSO4I
1886          MOLAL (7) = PSI3 + PSI1                     ! NO3I
1887          CNH42S4   = CHI2 - PSI2                     ! Solid (NH4)2SO4
1888          CNH4NO3   = ZERO                            ! Solid NH4NO3
1889          GHNO3     = CHI3 - PSI3                     ! Gas HNO3
1890          GNH3      = CHI4 - PSI4                     ! Gas NH3
1891          CALL CALCMR                                 ! Water content
1893 ! *** CALCULATE ACTIVITIES OR TERMINATE INTERNAL LOOP *****************
1895          IF (FRST.AND.CALAOU .OR. .NOT.FRST.AND.CALAIN) THEN
1896             CALL CALCACT
1897          ELSE
1898             GOTO 20
1899          ENDIF
1900 10    CONTINUE
1902 ! *** CALCULATE OBJECTIVE FUNCTION ************************************
1904 20    CONTINUE
1905 !CC      FUNCD3= NH4I/HI/MAX(GNH3,TINY)/A4 - ONE
1906       FUNCD3= MOLAL(3)/MOLAL(1)/MAX(GNH3,TINY)/A4 - ONE
1907       RETURN
1909 ! *** END OF FUNCTION FUNCD3 ********************************************
1911       END FUNCTION FUNCD3     
1912 !=======================================================================
1914 ! *** ISORROPIA CODE
1915 ! *** SUBROUTINE CALCD2
1916 ! *** CASE D2
1918 !     THE MAIN CHARACTERISTICS OF THIS REGIME ARE:
1919 !     1. SULFATE POOR (SULRAT > 2.0)
1920 !     2. THERE IS BOTH A LIQUID & SOLID PHASE
1921 !     3. SOLIDS POSSIBLE : (NH4)2SO4
1923 ! *** COPYRIGHT 1996-2006, UNIVERSITY OF MIAMI, CARNEGIE MELLON UNIVERSITY,
1924 ! *** GEORGIA INSTITUTE OF TECHNOLOGY
1925 ! *** WRITTEN BY ATHANASIOS NENES
1926 ! *** UPDATED BY CHRISTOS FOUNTOUKIS
1928 !=======================================================================
1930       SUBROUTINE CALCD2
1931       USE ISRPIA
1932       USE SOLUT
1933       IMPLICIT DOUBLE PRECISION (A-H, O-Z)
1934 !      INCLUDE 'ISRPIA.INC'
1935 !      COMMON /SOLUT/ CHI1, CHI2, CHI3, CHI4, CHI5, CHI6, CHI7, CHI8,   &
1936 !                     PSI1, PSI2, PSI3, PSI4, PSI5, PSI6, PSI7, PSI8,   &
1937 !                     A1,   A2,   A3,   A4,   A5,   A6,   A7,   A8
1939 ! *** FIND DRY COMPOSITION **********************************************
1941       CALL CALCD1A
1943 ! *** SETUP PARAMETERS ************************************************
1945       CHI1 = CNH4NO3               ! Save from CALCD1 run
1946       CHI2 = CNH42S4
1947       CHI3 = GHNO3
1948       CHI4 = GNH3
1950       PSI1 = CNH4NO3               ! ASSIGN INITIAL PSI's
1951       PSI2 = CNH42S4
1952       PSI3 = ZERO
1953       PSI4 = ZERO
1955       MOLAL(5) = ZERO
1956       MOLAL(6) = ZERO
1957       MOLAL(3) = PSI1
1958       MOLAL(7) = PSI1
1959       CALL CALCMR                  ! Initial water
1961       CALAOU = .TRUE.              ! Outer loop activity calculation flag
1962       PSI4LO = TINY                ! Low  limit
1963       PSI4HI = CHI4                ! High limit
1965 ! *** INITIAL VALUES FOR BISECTION ************************************
1967 60    X1 = PSI4LO
1968       Y1 = FUNCD2 (X1)
1969       IF (ABS(Y1).LE.EPS) RETURN
1970       YLO= Y1                 ! Save Y-value at HI position
1972 ! *** ROOT TRACKING ; FOR THE RANGE OF HI AND LO **********************
1974       DX   = (PSI4HI-PSI4LO)/REAL(NDIV)
1975       DO 10 I=1,NDIV
1976          X2 = X1+DX
1977          Y2 = FUNCD2 (X2)
1978          IF (SIGN(1.d0,Y1)*SIGN(1.d0,Y2).LT.ZERO) THEN
1980 ! This is done, in case if Y(PSI4LO)>0, but Y(PSI4LO+DX) < 0 (i.e.undersat)
1982              IF (Y1 .LE. Y2) GOTO 20  ! (Y1*Y2.LT.ZERO)
1983          ENDIF
1984          X1 = X2
1985          Y1 = Y2
1986 10    CONTINUE
1988 ! *** NO SUBDIVISION WITH SOLUTION FOUND
1990       YHI= Y1                      ! Save Y-value at Hi position
1991       IF (ABS(Y2) .LT. EPS) THEN   ! X2 IS A SOLUTION
1992          RETURN
1994 ! *** { YLO, YHI } < 0.0 THE SOLUTION IS ALWAYS UNDERSATURATED WITH NH3
1995 ! Physically I dont know when this might happen, but I have put this
1996 ! branch in for completeness. I assume there is no solution; all NO3 goes to the
1997 ! gas phase.
1999       ELSE IF (YLO.LT.ZERO .AND. YHI.LT.ZERO) THEN
2000          P4 = TINY ! PSI4LO ! CHI4
2001          YY = FUNCD2(P4)
2002          GOTO 50
2004 ! *** { YLO, YHI } > 0.0 THE SOLUTION IS ALWAYS SUPERSATURATED WITH NH3
2005 ! This happens when Sul.Rat. = 2.0, so some NH4+ from sulfate evaporates
2006 ! and goes to the gas phase ; so I redefine the LO and HI limits of PSI4
2007 ! and proceed again with root tracking.
2009       ELSE IF (YLO.GT.ZERO .AND. YHI.GT.ZERO) THEN
2010          PSI4HI = PSI4LO
2011          PSI4LO = PSI4LO - 0.1*(PSI1+PSI2) ! No solution; some NH3 evaporates
2012          IF (PSI4LO.LT.-(PSI1+PSI2)) THEN
2013             CALL PUSHERR (0001, 'CALCD2')  ! WARNING ERROR: NO SOLUTION
2014             RETURN
2015          ELSE
2016             MOLAL(5) = ZERO
2017             MOLAL(6) = ZERO
2018             MOLAL(3) = PSI1
2019             MOLAL(7) = PSI1
2020             CALL CALCMR                  ! Initial water
2021             GOTO 60                        ! Redo root tracking
2022          ENDIF
2023       ENDIF
2025 ! *** PERFORM BISECTION ***********************************************
2027 20    DO 30 I=1,MAXIT
2028          X3 = 0.5*(X1+X2)
2029          Y3 = FUNCD2 (X3)
2030          IF (SIGN(1.d0,Y1)*SIGN(1.d0,Y3) .LE. ZERO) THEN  ! (Y1*Y3 .LE. ZERO)
2031             Y2    = Y3
2032             X2    = X3
2033          ELSE
2034             Y1    = Y3
2035             X1    = X3
2036          ENDIF
2037          IF (ABS(X2-X1) .LE. EPS*ABS(X1)) GOTO 40
2038 30    CONTINUE
2039       CALL PUSHERR (0002, 'CALCD2')    ! WARNING ERROR: NO CONVERGENCE
2041 ! *** CONVERGED ; RETURN **********************************************
2043 40    X3 = MIN(X1,X2)   ! 0.5*(X1+X2)  ! Get "low" side, it's acidic soln.
2044       Y3 = FUNCD2 (X3)
2046 ! *** CALCULATE HSO4 SPECIATION AND RETURN *******************************
2048 50    CONTINUE
2049       IF (MOLAL(1).GT.TINY) THEN
2050          CALL CALCHS4 (MOLAL(1), MOLAL(5), ZERO, DELTA)
2051          MOLAL(1) = MOLAL(1) - DELTA                     ! H+   EFFECT
2052          MOLAL(5) = MOLAL(5) - DELTA                     ! SO4  EFFECT
2053          MOLAL(6) = DELTA                                ! HSO4 EFFECT
2054       ENDIF
2055       RETURN
2057 ! *** END OF SUBROUTINE CALCD2 ******************************************
2059       END SUBROUTINE CALCD2
2063 !=======================================================================
2065 ! *** ISORROPIA CODE
2066 ! *** FUNCTION FUNCD2
2067 ! *** CASE D2
2068 !     FUNCTION THAT SOLVES THE SYSTEM OF EQUATIONS FOR CASE D2 ;
2069 !     AND RETURNS THE VALUE OF THE ZEROED FUNCTION IN FUNCD2.
2071 !=======================================================================
2073       DOUBLE PRECISION FUNCTION FUNCD2 (P4)
2074       USE ISRPIA
2075       USE SOLUT
2076       IMPLICIT DOUBLE PRECISION (A-H, O-Z)
2077 !      INCLUDE 'ISRPIA.INC'
2078 !      COMMON /SOLUT/ CHI1, CHI2, CHI3, CHI4, CHI5, CHI6, CHI7, CHI8,   &
2079 !                     PSI1, PSI2, PSI3, PSI4, PSI5, PSI6, PSI7, PSI8,   &
2080 !                     A1,   A2,   A3,   A4,   A5,   A6,   A7,   A8
2082 ! *** SETUP PARAMETERS ************************************************
2084       CALL RSTGAM       ! Reset activity coefficients to 0.1
2085       FRST   = .TRUE.
2086       CALAIN = .TRUE.
2087       PSI4   = P4
2088       PSI2   = CHI2
2090 ! *** SOLVE EQUATIONS ; WITH ITERATIONS FOR ACTIVITY COEF. ************
2092       DO 10 I=1,NSWEEP
2093          A2  = XK7*(WATER/GAMA(4))**3.0
2094          A3  = XK4*R*TEMP*(WATER/GAMA(10))**2.0
2095          A4  = (XK2/XKW)*R*TEMP*(GAMA(10)/GAMA(5))**2.0
2096          A7  = XKW *RH*WATER*WATER
2098          IF (CHI2.GT.TINY .AND. WATER.GT.TINY) THEN
2099             PSI14 = PSI1+PSI4
2100             CALL POLY3 (PSI14,0.25*PSI14**2.,-A2/4.D0, PSI2, ISLV)  ! PSI2
2101             IF (ISLV.EQ.0) THEN
2102                 PSI2 = MIN (PSI2, CHI2)
2103             ELSE
2104                 PSI2 = ZERO
2105             ENDIF
2106          ENDIF
2108          PSI3  = A3*A4*CHI3*(CHI4-PSI4) - PSI1*(2.D0*PSI2+PSI1+PSI4)
2109          PSI3  = PSI3/(A3*A4*(CHI4-PSI4) + 2.D0*PSI2+PSI1+PSI4)
2110 !cc         PSI3  = MIN(MAX(PSI3, ZERO), CHI3)
2112          BB   = PSI4-PSI3 ! (BB > 0, acidic solution, <0 alkaline)
2114 ! Do not change computation scheme for H+, all others did not work well.
2116          DENM = BB+SQRT(BB*BB + 4.d0*A7)
2117          IF (DENM.LE.TINY) THEN       ! Avoid overflow when HI->0
2118             ABB  = ABS(BB)
2119             DENM = (BB+ABB) + 2.d0*A7/ABB ! Taylor expansion of SQRT
2120          ENDIF
2121          AHI = 2.d0*A7/DENM
2123 ! *** SPECIATION & WATER CONTENT ***************************************
2125          MOLAL (1) = AHI                              ! HI
2126          MOLAL (3) = PSI1 + PSI4 + 2.D0*PSI2          ! NH4
2127          MOLAL (5) = PSI2                             ! SO4
2128          MOLAL (6) = ZERO                             ! HSO4
2129          MOLAL (7) = PSI3 + PSI1                      ! NO3
2130          CNH42S4   = CHI2 - PSI2                      ! Solid (NH4)2SO4
2131          CNH4NO3   = ZERO                             ! Solid NH4NO3
2132          GHNO3     = CHI3 - PSI3                      ! Gas HNO3
2133          GNH3      = CHI4 - PSI4                      ! Gas NH3
2134          CALL CALCMR                                  ! Water content
2136 ! *** CALCULATE ACTIVITIES OR TERMINATE INTERNAL LOOP *****************
2138          IF (FRST.AND.CALAOU .OR. .NOT.FRST.AND.CALAIN) THEN
2139             CALL CALCACT
2140          ELSE
2141             GOTO 20
2142          ENDIF
2143 10    CONTINUE
2145 ! *** CALCULATE OBJECTIVE FUNCTION ************************************
2147 20    CONTINUE
2148 !CC      FUNCD2= NH4I/HI/MAX(GNH3,TINY)/A4 - ONE
2149       FUNCD2= MOLAL(3)/MOLAL(1)/MAX(GNH3,TINY)/A4 - ONE
2150       RETURN
2152 ! *** END OF FUNCTION FUNCD2 ********************************************
2154       END FUNCTION FUNCD2     
2155 !=======================================================================
2157 ! *** ISORROPIA CODE
2158 ! *** SUBROUTINE CALCD1
2159 ! *** CASE D1
2161 !     THE MAIN CHARACTERISTICS OF THIS REGIME ARE:
2162 !     1. SULFATE POOR (SULRAT > 2.0)
2163 !     2. SOLID AEROSOL ONLY
2164 !     3. SOLIDS POSSIBLE : (NH4)2SO4, NH4NO3
2166 !     THERE ARE TWO REGIMES DEFINED BY RELATIVE HUMIDITY:
2167 !     1. RH < MDRH ; ONLY SOLID PHASE POSSIBLE (SUBROUTINE CALCD1A)
2168 !     2. RH >= MDRH ; LIQUID PHASE POSSIBLE (MDRH REGION)
2170 ! *** COPYRIGHT 1996-2006, UNIVERSITY OF MIAMI, CARNEGIE MELLON UNIVERSITY,
2171 ! *** GEORGIA INSTITUTE OF TECHNOLOGY
2172 ! *** WRITTEN BY ATHANASIOS NENES
2173 ! *** UPDATED BY CHRISTOS FOUNTOUKIS
2175 !=======================================================================
2177       SUBROUTINE CALCD1
2178       USE ISRPIA
2179       IMPLICIT DOUBLE PRECISION (A-H, O-Z)
2180 !      INCLUDE 'ISRPIA.INC'
2181       EXTERNAL CALCD1A, CALCD2
2183 ! *** REGIME DEPENDS UPON THE AMBIENT RELATIVE HUMIDITY *****************
2185       IF (RH.LT.DRMASAN) THEN
2186          SCASE = 'D1 ; SUBCASE 1'   ! SOLID PHASE ONLY POSSIBLE
2187          CALL CALCD1A
2188          SCASE = 'D1 ; SUBCASE 1'
2189       ELSE
2190          SCASE = 'D1 ; SUBCASE 2'   ! LIQUID & SOLID PHASE POSSIBLE
2191          CALL CALCMDRH (RH, DRMASAN, DRNH4NO3, CALCD1A, CALCD2)
2192          SCASE = 'D1 ; SUBCASE 2'
2193       ENDIF
2195       RETURN
2197 ! *** END OF SUBROUTINE CALCD1 ******************************************
2199       END SUBROUTINE CALCD1
2202 !=======================================================================
2204 ! *** ISORROPIA CODE
2205 ! *** SUBROUTINE CALCD1A
2206 ! *** CASE D1 ; SUBCASE 1
2208 !     THE MAIN CHARACTERISTICS OF THIS REGIME ARE:
2209 !     1. SULFATE POOR (SULRAT > 2.0)
2210 !     2. SOLID AEROSOL ONLY
2211 !     3. SOLIDS POSSIBLE : (NH4)2SO4, NH4NO3
2213 !     THE SOLID (NH4)2SO4 IS CALCULATED FROM THE SULFATES, WHILE NH4NO3
2214 !     IS CALCULATED FROM NH3-HNO3 EQUILIBRIUM. 'ZE' IS THE AMOUNT OF
2215 !     NH4NO3 THAT VOLATIZES WHEN ALL POSSILBE NH4NO3 IS INITIALLY IN
2216 !     THE SOLID PHASE.
2218 ! *** COPYRIGHT 1996-2006, UNIVERSITY OF MIAMI, CARNEGIE MELLON UNIVERSITY,
2219 ! *** GEORGIA INSTITUTE OF TECHNOLOGY
2220 ! *** WRITTEN BY ATHANASIOS NENES
2221 ! *** UPDATED BY CHRISTOS FOUNTOUKIS
2223 !=======================================================================
2225       SUBROUTINE CALCD1A
2226       USE ISRPIA
2227       IMPLICIT DOUBLE PRECISION (A-H, O-Z)
2228 !      INCLUDE 'ISRPIA.INC'
2230 ! *** SETUP PARAMETERS ************************************************
2232       PARM    = XK10/(R*TEMP)/(R*TEMP)
2234 ! *** CALCULATE NH4NO3 THAT VOLATIZES *********************************
2236       CNH42S4 = W(2)
2237       X       = MAX(ZERO, MIN(W(3)-2.0*CNH42S4, W(4)))  ! MAX NH4NO3
2238       PS      = MAX(W(3) - X - 2.0*CNH42S4, ZERO)
2239       OM      = MAX(W(4) - X, ZERO)
2241       OMPS    = OM+PS
2242       DIAK    = SQRT(OMPS*OMPS + 4.0*PARM)              ! DIAKRINOUSA
2243       ZE      = MIN(X, 0.5*(-OMPS + DIAK))              ! THETIKI RIZA
2245 ! *** SPECIATION *******************************************************
2247       CNH4NO3 = X  - ZE    ! Solid NH4NO3
2248       GNH3    = PS + ZE    ! Gas NH3
2249       GHNO3   = OM + ZE    ! Gas HNO3
2251       RETURN
2253 ! *** END OF SUBROUTINE CALCD1A *****************************************
2255       END SUBROUTINE CALCD1A
2256 !=======================================================================
2258 ! *** ISORROPIA CODE
2259 ! *** SUBROUTINE CALCG5
2260 ! *** CASE G5
2262 !     THE MAIN CHARACTERISTICS OF THIS REGIME ARE:
2263 !     1. SULFATE POOR (SULRAT > 2.0) ; SODIUM POOR (SODRAT < 2.0)
2264 !     2. THERE IS BOTH A LIQUID & SOLID PHASE
2265 !     3. SOLIDS POSSIBLE : (NH4)2SO4, NH4CL, NA2SO4
2267 ! *** COPYRIGHT 1996-2006, UNIVERSITY OF MIAMI, CARNEGIE MELLON UNIVERSITY,
2268 ! *** GEORGIA INSTITUTE OF TECHNOLOGY
2269 ! *** WRITTEN BY ATHANASIOS NENES
2270 ! *** UPDATED BY CHRISTOS FOUNTOUKIS
2272 !=======================================================================
2274       SUBROUTINE CALCG5
2275       USE ISRPIA
2276       USE CASEG 
2277       IMPLICIT DOUBLE PRECISION (A-H, O-Z)
2278 !      INCLUDE 'ISRPIA.INC'
2280 !      DOUBLE PRECISION LAMDA
2281 !      COMMON /CASEG/ CHI1, CHI2, CHI3, CHI4, CHI5, CHI6, LAMDA,   &
2282 !                     PSI1, PSI2, PSI3, PSI4, PSI5, PSI6, PSI7,   &
2283 !                     A1,   A2,   A3,   A4,   A5,   A6,   A7
2285 ! *** SETUP PARAMETERS ************************************************
2287       CALAOU = .TRUE.
2288       CHI1   = 0.5*W(1)
2289       CHI2   = MAX (W(2)-CHI1, ZERO)
2290       CHI3   = ZERO
2291       CHI4   = MAX (W(3)-2.D0*CHI2, ZERO)
2292       CHI5   = W(4)
2293       CHI6   = W(5)
2295       PSI1   = CHI1
2296       PSI2   = CHI2
2297       PSI6LO = TINY
2298       PSI6HI = CHI6-TINY    ! MIN(CHI6-TINY, CHI4)
2300       WATER  = CHI2/M0(4) + CHI1/M0(2)
2302 ! *** INITIAL VALUES FOR BISECTION ************************************
2304       X1 = PSI6LO
2305       Y1 = FUNCG5A (X1)
2306       IF (CHI6.LE.TINY) GOTO 50
2307 !cc      IF (ABS(Y1).LE.EPS .OR. CHI6.LE.TINY) GOTO 50
2308 !cc      IF (WATER .LE. TINY) RETURN                    ! No water
2310 ! *** ROOT TRACKING ; FOR THE RANGE OF HI AND LO **********************
2312       DX = (PSI6HI-PSI6LO)/REAL(NDIV)
2313       DO 10 I=1,NDIV
2314          X2 = X1+DX
2315          Y2 = FUNCG5A (X2)
2316          IF (SIGN(1.d0,Y1)*SIGN(1.d0,Y2).LT.ZERO) GOTO 20  ! (Y1*Y2.LT.ZERO)
2317          X1 = X2
2318          Y1 = Y2
2319 10    CONTINUE
2321 ! *** NO SUBDIVISION WITH SOLUTION; IF ABS(Y2)<EPS SOLUTION IS ASSUMED
2323       IF (ABS(Y2) .GT. EPS) Y2 = FUNCG5A (PSI6LO)
2324       GOTO 50
2326 ! *** PERFORM BISECTION ***********************************************
2328 20    DO 30 I=1,MAXIT
2329          X3 = 0.5*(X1+X2)
2330          Y3 = FUNCG5A (X3)
2331          IF (SIGN(1.d0,Y1)*SIGN(1.d0,Y3) .LE. ZERO) THEN  ! (Y1*Y3 .LE. ZERO)
2332             Y2    = Y3
2333             X2    = X3
2334          ELSE
2335             Y1    = Y3
2336             X1    = X3
2337          ENDIF
2338          IF (ABS(X2-X1) .LE. EPS*X1) GOTO 40
2339 30    CONTINUE
2340       CALL PUSHERR (0002, 'CALCG5')    ! WARNING ERROR: NO CONVERGENCE
2342 ! *** CONVERGED ; RETURN **********************************************
2344 40    X3 = 0.5*(X1+X2)
2345       Y3 = FUNCG5A (X3)
2347 ! *** CALCULATE HSO4 SPECIATION AND RETURN *******************************
2349 50    CONTINUE
2350       IF (MOLAL(1).GT.TINY .AND. MOLAL(5).GT.TINY) THEN  ! If quadrat.called
2351          CALL CALCHS4 (MOLAL(1), MOLAL(5), ZERO, DELTA)
2352          MOLAL(1) = MOLAL(1) - DELTA                    ! H+   EFFECT
2353          MOLAL(5) = MOLAL(5) - DELTA                    ! SO4  EFFECT
2354          MOLAL(6) = DELTA                               ! HSO4 EFFECT
2355       ENDIF
2357       RETURN
2359 ! *** END OF SUBROUTINE CALCG5 *******************************************
2361       END SUBROUTINE CALCG5
2366 !=======================================================================
2368 ! *** ISORROPIA CODE
2369 ! *** SUBROUTINE FUNCG5A
2370 ! *** CASE G5
2372 !     THE MAIN CHARACTERISTICS OF THIS REGIME ARE:
2373 !     1. SULFATE POOR (SULRAT > 2.0) ; SODIUM POOR (SODRAT < 2.0)
2374 !     2. THERE IS BOTH A LIQUID & SOLID PHASE
2375 !     3. SOLIDS POSSIBLE : (NH4)2SO4, NH4CL, NA2SO4
2377 ! *** COPYRIGHT 1996-2006, UNIVERSITY OF MIAMI, CARNEGIE MELLON UNIVERSITY,
2378 ! *** GEORGIA INSTITUTE OF TECHNOLOGY
2379 ! *** WRITTEN BY ATHANASIOS NENES
2380 ! *** UPDATED BY CHRISTOS FOUNTOUKIS
2382 !=======================================================================
2384       DOUBLE PRECISION FUNCTION FUNCG5A (X)
2385       USE ISRPIA
2386       USE CASEG 
2387       IMPLICIT DOUBLE PRECISION (A-H, O-Z)
2388 !      INCLUDE 'ISRPIA.INC'
2390 !      DOUBLE PRECISION LAMDA
2391 !      COMMON /CASEG/ CHI1, CHI2, CHI3, CHI4, CHI5, CHI6, LAMDA,   &
2392 !                     PSI1, PSI2, PSI3, PSI4, PSI5, PSI6, PSI7,   &
2393 !                     A1,   A2,   A3,   A4,   A5,   A6,   A7
2395 ! *** SETUP PARAMETERS ************************************************
2397       PSI6   = X
2398       FRST   = .TRUE.
2399       CALAIN = .TRUE.
2401 ! *** SOLVE EQUATIONS ; WITH ITERATIONS FOR ACTIVITY COEF. ************
2403       DO 10 I=1,NSWEEP
2405       A1  = XK5 *(WATER/GAMA(2))**3.0
2406       A2  = XK7 *(WATER/GAMA(4))**3.0
2407       A4  = (XK2/XKW)*R*TEMP*(GAMA(10)/GAMA(5))**2.0
2408       A5  = XK4 *R*TEMP*(WATER/GAMA(10))**2.0
2409       A6  = XK3 *R*TEMP*(WATER/GAMA(11))**2.0
2410       AKK = A4*A6
2412 !  CALCULATE DISSOCIATION QUANTITIES
2414       IF (CHI5.GE.TINY) THEN
2415          PSI5 = PSI6*CHI5/(A6/A5*(CHI6-PSI6) + PSI6)
2416       ELSE
2417          PSI5 = TINY
2418       ENDIF
2420 !CC      IF(CHI4.GT.TINY) THEN
2421       IF(W(2).GT.TINY) THEN       ! Accounts for NH3 evaporation
2422          BB   =-(CHI4 + PSI6 + PSI5 + 1.d0/A4)
2423          CC   = CHI4*(PSI5+PSI6) - 2.d0*PSI2/A4
2424          DD   = MAX(BB*BB-4.d0*CC,ZERO)           ! Patch proposed by Uma Shankar, 19/11/01
2425          PSI4 =0.5d0*(-BB - SQRT(DD))
2426       ELSE
2427          PSI4 = TINY
2428       ENDIF
2430 ! *** CALCULATE SPECIATION ********************************************
2432       MOLAL (2) = 2.0D0*PSI1                          ! NAI
2433       MOLAL (3) = 2.0*PSI2 + PSI4                     ! NH4I
2434       MOLAL (4) = PSI6                                ! CLI
2435       MOLAL (5) = PSI2 + PSI1                         ! SO4I
2436       MOLAL (6) = ZERO
2437       MOLAL (7) = PSI5                                ! NO3I
2439       SMIN      = 2.d0*MOLAL(5)+MOLAL(7)+MOLAL(4)-MOLAL(2)-MOLAL(3)
2440       CALL CALCPH (SMIN, HI, OHI)
2441       MOLAL (1) = HI
2443       GNH3      = MAX(CHI4 - PSI4, TINY)              ! Gas NH3
2444       GHNO3     = MAX(CHI5 - PSI5, TINY)              ! Gas HNO3
2445       GHCL      = MAX(CHI6 - PSI6, TINY)              ! Gas HCl
2447       CNH42S4   = ZERO                                ! Solid (NH4)2SO4
2448       CNH4NO3   = ZERO                                ! Solid NH4NO3
2449       CNH4CL    = ZERO                                ! Solid NH4Cl
2451       CALL CALCMR                                     ! Water content
2453 ! *** CALCULATE ACTIVITIES OR TERMINATE INTERNAL LOOP *****************
2455       IF (FRST.AND.CALAOU .OR. .NOT.FRST.AND.CALAIN) THEN
2456          CALL CALCACT
2457       ELSE
2458          GOTO 20
2459       ENDIF
2460 10    CONTINUE
2462 ! *** CALCULATE FUNCTION VALUE FOR OUTER LOOP ***************************
2464 20    FUNCG5A = MOLAL(1)*MOLAL(4)/GHCL/A6 - ONE
2465 !CC         FUNCG5A = MOLAL(3)*MOLAL(4)/GHCL/GNH3/A6/A4 - ONE
2467       RETURN
2469 ! *** END OF FUNCTION FUNCG5A *******************************************
2471       END FUNCTION FUNCG5A    
2473 !=======================================================================
2475 ! *** ISORROPIA CODE
2476 ! *** SUBROUTINE CALCG4
2477 ! *** CASE G4
2479 !     THE MAIN CHARACTERISTICS OF THIS REGIME ARE:
2480 !     1. SULFATE POOR (SULRAT > 2.0) ; SODIUM POOR (SODRAT < 2.0)
2481 !     2. THERE IS BOTH A LIQUID & SOLID PHASE
2482 !     3. SOLIDS POSSIBLE : (NH4)2SO4, NH4CL, NA2SO4
2484 ! *** COPYRIGHT 1996-2006, UNIVERSITY OF MIAMI, CARNEGIE MELLON UNIVERSITY,
2485 ! *** GEORGIA INSTITUTE OF TECHNOLOGY
2486 ! *** WRITTEN BY ATHANASIOS NENES
2487 ! *** UPDATED BY CHRISTOS FOUNTOUKIS
2489 !=======================================================================
2491       SUBROUTINE CALCG4
2492       USE ISRPIA
2493       USE CASEG 
2494       IMPLICIT DOUBLE PRECISION (A-H, O-Z)
2495 !      INCLUDE 'ISRPIA.INC'
2497 !      DOUBLE PRECISION LAMDA
2498 !      COMMON /CASEG/ CHI1, CHI2, CHI3, CHI4, CHI5, CHI6, LAMDA,   &
2499 !                     PSI1, PSI2, PSI3, PSI4, PSI5, PSI6, PSI7,   &
2500 !                     A1,   A2,   A3,   A4,   A5,   A6,   A7
2502 ! *** SETUP PARAMETERS ************************************************
2504       CALAOU = .TRUE.
2505       CHI1   = 0.5*W(1)
2506       CHI2   = MAX (W(2)-CHI1, ZERO)
2507       CHI3   = ZERO
2508       CHI4   = MAX (W(3)-2.D0*CHI2, ZERO)
2509       CHI5   = W(4)
2510       CHI6   = W(5)
2512       PSI2   = CHI2
2513       PSI6LO = TINY
2514       PSI6HI = CHI6-TINY    ! MIN(CHI6-TINY, CHI4)
2516       WATER  = CHI2/M0(4) + CHI1/M0(2)
2518 ! *** INITIAL VALUES FOR BISECTION ************************************
2520       X1 = PSI6LO
2521       Y1 = FUNCG4A (X1)
2522       IF (CHI6.LE.TINY) GOTO 50
2523 !CC      IF (ABS(Y1).LE.EPS .OR. CHI6.LE.TINY .OR. WATER .LE. TINY) GOTO 50
2524 !CC      IF (WATER .LE. TINY) RETURN                    ! No water
2526 ! *** ROOT TRACKING ; FOR THE RANGE OF HI AND LO **********************
2528       DX = (PSI6HI-PSI6LO)/REAL(NDIV)
2529       DO 10 I=1,NDIV
2530          X2  = X1+DX
2531          Y2  = FUNCG4A (X2)
2532          IF (SIGN(1.d0,Y1)*SIGN(1.d0,Y2).LT.ZERO) GOTO 20  ! (Y1*Y2.LT.ZERO)
2533          X1  = X2
2534          Y1  = Y2
2535 10    CONTINUE
2537 ! *** NO SUBDIVISION WITH SOLUTION; IF ABS(Y2)<EPS SOLUTION IS ASSUMED
2539       IF (ABS(Y2) .GT. EPS) Y2 = FUNCG4A (PSI6LO)
2540       GOTO 50
2542 ! *** PERFORM BISECTION ***********************************************
2544 20    DO 30 I=1,MAXIT
2545          X3 = 0.5*(X1+X2)
2546          Y3 = FUNCG4A (X3)
2547          IF (SIGN(1.d0,Y1)*SIGN(1.d0,Y3) .LE. ZERO) THEN  ! (Y1*Y3 .LE. ZERO)
2548             Y2    = Y3
2549             X2    = X3
2550          ELSE
2551             Y1    = Y3
2552             X1    = X3
2553          ENDIF
2554          IF (ABS(X2-X1) .LE. EPS*X1) GOTO 40
2555 30    CONTINUE
2556       CALL PUSHERR (0002, 'CALCG4')    ! WARNING ERROR: NO CONVERGENCE
2558 ! *** CONVERGED ; RETURN **********************************************
2560 40    X3 = 0.5*(X1+X2)
2561       Y3 = FUNCG4A (X3)
2563 ! *** CALCULATE HSO4 SPECIATION AND RETURN *******************************
2565 50    CONTINUE
2566       IF (MOLAL(1).GT.TINY .AND. MOLAL(5).GT.TINY) THEN
2567          CALL CALCHS4 (MOLAL(1), MOLAL(5), ZERO, DELTA)
2568          MOLAL(1) = MOLAL(1) - DELTA                     ! H+   EFFECT
2569          MOLAL(5) = MOLAL(5) - DELTA                     ! SO4  EFFECT
2570          MOLAL(6) = DELTA                                ! HSO4 EFFECT
2571       ENDIF
2573       RETURN
2575 ! *** END OF SUBROUTINE CALCG4 *******************************************
2577       END SUBROUTINE CALCG4
2582 !=======================================================================
2584 ! *** ISORROPIA CODE
2585 ! *** SUBROUTINE FUNCG4A
2586 ! *** CASE G4
2588 !     THE MAIN CHARACTERISTICS OF THIS REGIME ARE:
2589 !     1. SULFATE POOR (SULRAT > 2.0) ; SODIUM POOR (SODRAT < 2.0)
2590 !     2. THERE IS BOTH A LIQUID & SOLID PHASE
2591 !     3. SOLIDS POSSIBLE : (NH4)2SO4, NH4CL, NA2SO4
2593 ! *** COPYRIGHT 1996-2006, UNIVERSITY OF MIAMI, CARNEGIE MELLON UNIVERSITY,
2594 ! *** GEORGIA INSTITUTE OF TECHNOLOGY
2595 ! *** WRITTEN BY ATHANASIOS NENES
2596 ! *** UPDATED BY CHRISTOS FOUNTOUKIS
2598 !=======================================================================
2600       DOUBLE PRECISION FUNCTION FUNCG4A (X)
2601       USE ISRPIA
2602       USE CASEG
2603       IMPLICIT DOUBLE PRECISION (A-H, O-Z)
2604 !      INCLUDE 'ISRPIA.INC'
2606 !      DOUBLE PRECISION LAMDA, NAI, NH4I, NO3I
2607       DOUBLE PRECISION  NAI, NH4I, NO3I
2608 !      COMMON /CASEG/ CHI1, CHI2, CHI3, CHI4, CHI5, CHI6, LAMDA,   &
2609 !                     PSI1, PSI2, PSI3, PSI4, PSI5, PSI6, PSI7,   &
2610 !                     A1,   A2,   A3,   A4,   A5,   A6,   A7
2612 ! *** SETUP PARAMETERS ************************************************
2614       PSI6   = X
2615       PSI1   = CHI1
2616       FRST   = .TRUE.
2617       CALAIN = .TRUE.
2619 ! *** SOLVE EQUATIONS ; WITH ITERATIONS FOR ACTIVITY COEF. ************
2621       DO 10 I=1,NSWEEP
2623       A1  = XK5 *(WATER/GAMA(2))**3.0
2624       A2  = XK7 *(WATER/GAMA(4))**3.0
2625       A4  = (XK2/XKW)*R*TEMP*(GAMA(10)/GAMA(5))**2.0
2626       A5  = XK4 *R*TEMP*(WATER/GAMA(10))**2.0
2627       A6  = XK3 *R*TEMP*(WATER/GAMA(11))**2.0
2629 !  CALCULATE DISSOCIATION QUANTITIES
2631       IF (CHI5.GE.TINY) THEN
2632          PSI5 = PSI6*CHI5/(A6/A5*(CHI6-PSI6) + PSI6)
2633       ELSE
2634          PSI5 = TINY
2635       ENDIF
2637 !CC      IF(CHI4.GT.TINY) THEN
2638       IF(W(2).GT.TINY) THEN       ! Accounts for NH3 evaporation
2639          BB   =-(CHI4 + PSI6 + PSI5 + 1.d0/A4)
2640          CC   = CHI4*(PSI5+PSI6) - 2.d0*PSI2/A4
2641          DD   = MAX(BB*BB-4.d0*CC,ZERO) ! Patch proposed by Uma shankar, 19/11/2001
2642          PSI4 =0.5d0*(-BB - SQRT(DD))
2643       ELSE
2644          PSI4 = TINY
2645       ENDIF
2647 !  CALCULATE CONCENTRATIONS
2649       NH4I = 2.0*PSI2 + PSI4
2650       CLI  = PSI6
2651       SO4I = PSI2 + PSI1
2652       NO3I = PSI5
2653       NAI  = 2.0D0*PSI1
2655       CALL CALCPH(2.d0*SO4I+NO3I+CLI-NAI-NH4I, HI, OHI)
2657 ! *** Na2SO4 DISSOLUTION
2659       IF (CHI1.GT.TINY .AND. WATER.GT.TINY) THEN        ! PSI1
2660          CALL POLY3 (PSI2, ZERO, -A1/4.D0, PSI1, ISLV)
2661          IF (ISLV.EQ.0) THEN
2662              PSI1 = MIN (PSI1, CHI1)
2663          ELSE
2664              PSI1 = ZERO
2665          ENDIF
2666       ELSE
2667          PSI1 = ZERO
2668       ENDIF
2670 ! *** SAVE CONCENTRATIONS IN MOLAL ARRAY ******************************
2672       MOLAL (1) = HI
2673       MOLAL (2) = NAI
2674       MOLAL (3) = NH4I
2675       MOLAL (4) = CLI
2676       MOLAL (5) = SO4I
2677       MOLAL (6) = ZERO
2678       MOLAL (7) = NO3I
2680 ! *** CALCULATE GAS / SOLID SPECIES (LIQUID IN MOLAL ALREADY) *********
2682       GNH3      = MAX(CHI4 - PSI4, TINY)
2683       GHNO3     = MAX(CHI5 - PSI5, TINY)
2684       GHCL      = MAX(CHI6 - PSI6, TINY)
2686       CNH42S4   = ZERO
2687       CNH4NO3   = ZERO
2688       CNH4CL    = ZERO
2689       CNA2SO4   = MAX(CHI1-PSI1,ZERO)
2691 ! *** CALCULATE MOLALR ARRAY, WATER AND ACTIVITIES **********************
2693       CALL CALCMR
2695 ! *** CALCULATE ACTIVITIES OR TERMINATE INTERNAL LOOP *****************
2697       IF (FRST.AND.CALAOU .OR. .NOT.FRST.AND.CALAIN) THEN
2698          CALL CALCACT
2699       ELSE
2700          GOTO 20
2701       ENDIF
2702 10    CONTINUE
2704 ! *** CALCULATE FUNCTION VALUE FOR OUTER LOOP ***************************
2706 20    FUNCG4A = MOLAL(1)*MOLAL(4)/GHCL/A6 - ONE
2707 !CC         FUNCG4A = MOLAL(3)*MOLAL(4)/GHCL/GNH3/A6/A4 - ONE
2709       RETURN
2711 ! *** END OF FUNCTION FUNCG4A *******************************************
2713       END FUNCTION FUNCG4A    
2715 !=======================================================================
2717 ! *** ISORROPIA CODE
2718 ! *** SUBROUTINE CALCG3
2719 ! *** CASE G3
2721 !     THE MAIN CHARACTERISTICS OF THIS REGIME ARE:
2722 !     1. SULFATE POOR (SULRAT > 2.0) ; SODIUM POOR (SODRAT < 2.0)
2723 !     2. LIQUID & SOLID PHASE ARE BOTH POSSIBLE
2724 !     3. SOLIDS POSSIBLE : (NH4)2SO4, NH4CL, NA2SO4
2726 ! *** COPYRIGHT 1996-2006, UNIVERSITY OF MIAMI, CARNEGIE MELLON UNIVERSITY,
2727 ! *** GEORGIA INSTITUTE OF TECHNOLOGY
2728 ! *** WRITTEN BY ATHANASIOS NENES
2729 ! *** UPDATED BY CHRISTOS FOUNTOUKIS
2731 !=======================================================================
2733       SUBROUTINE CALCG3
2734       USE ISRPIA
2735       IMPLICIT DOUBLE PRECISION (A-H, O-Z)
2736 !      INCLUDE 'ISRPIA.INC'
2737       EXTERNAL CALCG1A, CALCG4
2739 ! *** REGIME DEPENDS ON THE EXISTANCE OF WATER AND OF THE RH ************
2741       IF (W(4).GT.TINY .AND. W(5).GT.TINY) THEN ! NO3,CL EXIST, WATER POSSIBLE
2742          SCASE = 'G3 ; SUBCASE 1'
2743          CALL CALCG3A
2744          SCASE = 'G3 ; SUBCASE 1'
2745       ELSE                                      ! NO3, CL NON EXISTANT
2746          SCASE = 'G1 ; SUBCASE 1'
2747          CALL CALCG1A
2748          SCASE = 'G1 ; SUBCASE 1'
2749       ENDIF
2751       IF (WATER.LE.TINY) THEN
2752          IF (RH.LT.DRMG3) THEN        ! ONLY SOLIDS
2753             WATER = TINY
2754             DO 10 I=1,NIONS
2755                MOLAL(I) = ZERO
2756 10          CONTINUE
2757             CALL CALCG1A
2758             SCASE = 'G3 ; SUBCASE 2'
2759             RETURN
2760          ELSE
2761             SCASE = 'G3 ; SUBCASE 3'  ! MDRH REGION (NA2SO4, NH42S4)
2762             CALL CALCMDRH (RH, DRMG3, DRNH42S4, CALCG1A, CALCG4)
2763             SCASE = 'G3 ; SUBCASE 3'
2764          ENDIF
2765       ENDIF
2767       RETURN
2769 ! *** END OF SUBROUTINE CALCG3 ******************************************
2771       END SUBROUTINE CALCG3
2774 !=======================================================================
2776 ! *** ISORROPIA CODE
2777 ! *** SUBROUTINE CALCG3A
2778 ! *** CASE G3 ; SUBCASE 1
2780 !     THE MAIN CHARACTERISTICS OF THIS REGIME ARE:
2781 !     1. SULFATE POOR (SULRAT > 2.0) ; SODIUM POOR (SODRAT < 2.0)
2782 !     2. THERE IS BOTH A LIQUID & SOLID PHASE
2783 !     3. SOLIDS POSSIBLE : (NH4)2SO4, NH4CL, NA2SO4
2785 ! *** COPYRIGHT 1996-2006, UNIVERSITY OF MIAMI, CARNEGIE MELLON UNIVERSITY,
2786 ! *** GEORGIA INSTITUTE OF TECHNOLOGY
2787 ! *** WRITTEN BY ATHANASIOS NENES
2788 ! *** UPDATED BY CHRISTOS FOUNTOUKIS
2790 !=======================================================================
2792       SUBROUTINE CALCG3A
2793       USE ISRPIA
2794       USE CASEG
2795       IMPLICIT DOUBLE PRECISION (A-H, O-Z)
2796 !      INCLUDE 'ISRPIA.INC'
2798 !      DOUBLE PRECISION LAMDA
2799 !      COMMON /CASEG/ CHI1, CHI2, CHI3, CHI4, CHI5, CHI6, LAMDA,   &
2800 !                     PSI1, PSI2, PSI3, PSI4, PSI5, PSI6, PSI7,   &
2801 !                     A1,   A2,   A3,   A4,   A5,   A6,   A7
2803 ! *** SETUP PARAMETERS ************************************************
2805       CALAOU = .TRUE.
2806       CHI1   = 0.5*W(1)
2807       CHI2   = MAX (W(2)-CHI1, ZERO)
2808       CHI3   = ZERO
2809       CHI4   = MAX (W(3)-2.D0*CHI2, ZERO)
2810       CHI5   = W(4)
2811       CHI6   = W(5)
2813       PSI6LO = TINY
2814       PSI6HI = CHI6-TINY    ! MIN(CHI6-TINY, CHI4)
2816       WATER  = TINY
2818 ! *** INITIAL VALUES FOR BISECTION ************************************
2820       X1 = PSI6LO
2821       Y1 = FUNCG3A (X1)
2822       IF (CHI6.LE.TINY) GOTO 50
2823 !CC      IF (ABS(Y1).LE.EPS .OR. CHI6.LE.TINY .OR. WATER .LE. TINY) GOTO 50
2824 !CC      IF (WATER .LE. TINY) RETURN                    ! No water
2826 ! *** ROOT TRACKING ; FOR THE RANGE OF HI AND LO **********************
2828       DX = (PSI6HI-PSI6LO)/REAL(NDIV)
2829       DO 10 I=1,NDIV
2830          X2  = X1+DX
2831          Y2  = FUNCG3A (X2)
2833          IF (SIGN(1.d0,Y1)*SIGN(1.d0,Y2).LT.ZERO) GOTO 20  ! (Y1*Y2.LT.ZERO)
2834          X1  = X2
2835          Y1  = Y2
2836 10    CONTINUE
2838 ! *** NO SUBDIVISION WITH SOLUTION; IF ABS(Y2)<EPS SOLUTION IS ASSUMED
2840       IF (ABS(Y2) .GT. EPS) Y2 = FUNCG3A (PSI6LO)
2841       GOTO 50
2843 ! *** PERFORM BISECTION ***********************************************
2845 20    DO 30 I=1,MAXIT
2846          X3 = 0.5*(X1+X2)
2847          Y3 = FUNCG3A (X3)
2848          IF (SIGN(1.d0,Y1)*SIGN(1.d0,Y3) .LE. ZERO) THEN  ! (Y1*Y3 .LE. ZERO)
2849             Y2    = Y3
2850             X2    = X3
2851          ELSE
2852             Y1    = Y3
2853             X1    = X3
2854          ENDIF
2855          IF (ABS(X2-X1) .LE. EPS*X1) GOTO 40
2856 30    CONTINUE
2857       CALL PUSHERR (0002, 'CALCG3A')    ! WARNING ERROR: NO CONVERGENCE
2859 ! *** CONVERGED ; RETURN **********************************************
2861 40    X3 = 0.5*(X1+X2)
2862       Y3 = FUNCG3A (X3)
2864 ! *** FINAL CALCULATIONS *************************************************
2866 50    CONTINUE
2868 ! *** Na2SO4 DISSOLUTION
2870       IF (CHI1.GT.TINY .AND. WATER.GT.TINY) THEN        ! PSI1
2871          CALL POLY3 (PSI2, ZERO, -A1/4.D0, PSI1, ISLV)
2872          IF (ISLV.EQ.0) THEN
2873              PSI1 = MIN (PSI1, CHI1)
2874          ELSE
2875              PSI1 = ZERO
2876          ENDIF
2877       ELSE
2878          PSI1 = ZERO
2879       ENDIF
2880       MOLAL(2) = 2.0D0*PSI1               ! Na+  EFFECT
2881       MOLAL(5) = MOLAL(5) + PSI1          ! SO4  EFFECT
2882       CNA2SO4  = MAX(CHI1 - PSI1, ZERO)   ! NA2SO4(s) depletion
2884 ! *** HSO4 equilibrium
2886       IF (MOLAL(1).GT.TINY .AND. MOLAL(5).GT.TINY) THEN
2887          CALL CALCHS4 (MOLAL(1), MOLAL(5), ZERO, DELTA)
2888          MOLAL(1) = MOLAL(1) - DELTA                     ! H+   EFFECT
2889          MOLAL(5) = MOLAL(5) - DELTA                     ! SO4  EFFECT
2890          MOLAL(6) = DELTA                                ! HSO4 EFFECT
2891       ENDIF
2893       RETURN
2895 ! *** END OF SUBROUTINE CALCG3A ******************************************
2897       END SUBROUTINE CALCG3A
2902 !=======================================================================
2904 ! *** ISORROPIA CODE
2905 ! *** SUBROUTINE FUNCG3A
2906 ! *** CASE G3 ; SUBCASE 1
2908 !     THE MAIN CHARACTERISTICS OF THIS REGIME ARE:
2909 !     1. SULFATE POOR (SULRAT > 2.0) ; SODIUM POOR (SODRAT < 2.0)
2910 !     2. THERE IS BOTH A LIQUID & SOLID PHASE
2911 !     3. SOLIDS POSSIBLE : (NH4)2SO4, NH4CL, NA2SO4
2913 ! *** COPYRIGHT 1996-2006, UNIVERSITY OF MIAMI, CARNEGIE MELLON UNIVERSITY,
2914 ! *** GEORGIA INSTITUTE OF TECHNOLOGY
2915 ! *** WRITTEN BY ATHANASIOS NENES
2916 ! *** UPDATED BY CHRISTOS FOUNTOUKIS
2918 !=======================================================================
2920       DOUBLE PRECISION FUNCTION FUNCG3A (X)
2921       USE ISRPIA
2922       USE CASEG
2923       IMPLICIT DOUBLE PRECISION (A-H, O-Z)
2924 !      INCLUDE 'ISRPIA.INC'
2926 !      DOUBLE PRECISION LAMDA
2927 !      COMMON /CASEG/ CHI1, CHI2, CHI3, CHI4, CHI5, CHI6, LAMDA,   &
2928 !                     PSI1, PSI2, PSI3, PSI4, PSI5, PSI6, PSI7,   &
2929 !                     A1,   A2,   A3,   A4,   A5,   A6,   A7
2931 ! *** SETUP PARAMETERS ************************************************
2933       PSI6   = X
2934       PSI2   = CHI2
2935       FRST   = .TRUE.
2936       CALAIN = .TRUE.
2938 ! *** SOLVE EQUATIONS ; WITH ITERATIONS FOR ACTIVITY COEF. ************
2940       DO 10 I=1,NSWEEP
2942       A1  = XK5 *(WATER/GAMA(2))**3.0
2943       A2  = XK7 *(WATER/GAMA(4))**3.0
2944       A4  = (XK2/XKW)*R*TEMP*(GAMA(10)/GAMA(5))**2.0
2945       A5  = XK4 *R*TEMP*(WATER/GAMA(10))**2.0
2946       A6  = XK3 *R*TEMP*(WATER/GAMA(11))**2.0
2948 !  CALCULATE DISSOCIATION QUANTITIES
2950       IF (CHI5.GE.TINY) THEN
2951          PSI5 = PSI6*CHI5/(A6/A5*(CHI6-PSI6) + PSI6)
2952       ELSE
2953          PSI5 = TINY
2954       ENDIF
2956 !CC      IF(CHI4.GT.TINY) THEN
2957       IF(W(2).GT.TINY) THEN       ! Accounts for NH3 evaporation
2958          BB   =-(CHI4 + PSI6 + PSI5 + 1.d0/A4)
2959          CC   = CHI4*(PSI5+PSI6) - 2.d0*PSI2/A4
2960          DD   = MAX(BB*BB-4.d0*CC,ZERO)  ! Patch proposed by Uma Shankar, 19/11/01
2961          PSI4 =0.5d0*(-BB - SQRT(DD))
2962       ELSE
2963          PSI4 = TINY
2964       ENDIF
2966       IF (CHI2.GT.TINY .AND. WATER.GT.TINY) THEN
2967          CALL POLY3 (PSI4, PSI4*PSI4/4.D0, -A2/4.D0, PSI20, ISLV)
2968          IF (ISLV.EQ.0) PSI2 = MIN (PSI20, CHI2)
2969       ENDIF
2971 ! *** CALCULATE GAS / SOLID SPECIES (LIQUID IN MOLAL ALREADY) *********
2973       MOLAL (2) = ZERO                                ! Na
2974       MOLAL (3) = 2.0*PSI2 + PSI4                     ! NH4I
2975       MOLAL (4) = PSI6                                ! CLI
2976       MOLAL (5) = PSI2                                ! SO4I
2977       MOLAL (6) = ZERO                                ! HSO4
2978       MOLAL (7) = PSI5                                ! NO3I
2980       SMIN      = 2.d0*MOLAL(5)+MOLAL(7)+MOLAL(4)-MOLAL(2)-MOLAL(3)
2981       CALL CALCPH (SMIN, HI, OHI)
2982       MOLAL (1) = HI
2984       GNH3      = MAX(CHI4 - PSI4, TINY)              ! Gas NH3
2985       GHNO3     = MAX(CHI5 - PSI5, TINY)              ! Gas HNO3
2986       GHCL      = MAX(CHI6 - PSI6, TINY)              ! Gas HCl
2988       CNH42S4   = CHI2 - PSI2                         ! Solid (NH4)2SO4
2989       CNH4NO3   = ZERO                                ! Solid NH4NO3
2990       CNH4CL    = ZERO                                ! Solid NH4Cl
2992       CALL CALCMR                                     ! Water content
2994 ! *** CALCULATE ACTIVITIES OR TERMINATE INTERNAL LOOP *****************
2996       IF (FRST.AND.CALAOU .OR. .NOT.FRST.AND.CALAIN) THEN
2997          CALL CALCACT
2998       ELSE
2999          GOTO 20
3000       ENDIF
3001 10    CONTINUE
3003 ! *** CALCULATE FUNCTION VALUE FOR OUTER LOOP ***************************
3005 20    FUNCG3A = MOLAL(1)*MOLAL(4)/GHCL/A6 - ONE
3006 !CC         FUNCG3A = MOLAL(3)*MOLAL(4)/GHCL/GNH3/A6/A4 - ONE
3008       RETURN
3010 ! *** END OF FUNCTION FUNCG3A *******************************************
3012       END FUNCTION FUNCG3A    
3014 !=======================================================================
3016 ! *** ISORROPIA CODE
3017 ! *** SUBROUTINE CALCG2
3018 ! *** CASE G2
3020 !     THE MAIN CHARACTERISTICS OF THIS REGIME ARE:
3021 !     1. SULFATE POOR (SULRAT > 2.0) ; SODIUM POOR (SODRAT < 2.0)
3022 !     2. LIQUID & SOLID PHASE ARE BOTH POSSIBLE
3023 !     3. SOLIDS POSSIBLE : (NH4)2SO4, NH4CL, NA2SO4
3025 ! *** COPYRIGHT 1996-2006, UNIVERSITY OF MIAMI, CARNEGIE MELLON UNIVERSITY,
3026 ! *** GEORGIA INSTITUTE OF TECHNOLOGY
3027 ! *** WRITTEN BY ATHANASIOS NENES
3028 ! *** UPDATED BY CHRISTOS FOUNTOUKIS
3030 !=======================================================================
3032       SUBROUTINE CALCG2
3033       USE ISRPIA
3034       IMPLICIT DOUBLE PRECISION (A-H, O-Z)
3035 !      INCLUDE 'ISRPIA.INC'
3036       EXTERNAL CALCG1A, CALCG3A, CALCG4
3038 ! *** REGIME DEPENDS ON THE EXISTANCE OF NITRATES ***********************
3040       IF (W(4).GT.TINY) THEN        ! NO3 EXISTS, WATER POSSIBLE
3041          SCASE = 'G2 ; SUBCASE 1'
3042          CALL CALCG2A
3043          SCASE = 'G2 ; SUBCASE 1'
3044       ELSE                          ! NO3 NON EXISTANT, WATER NOT POSSIBLE
3045          SCASE = 'G1 ; SUBCASE 1'
3046          CALL CALCG1A
3047          SCASE = 'G1 ; SUBCASE 1'
3048       ENDIF
3050 ! *** REGIME DEPENDS ON THE EXISTANCE OF WATER AND OF THE RH ************
3052       IF (WATER.LE.TINY) THEN
3053          IF (RH.LT.DRMG2) THEN             ! ONLY SOLIDS
3054             WATER = TINY
3055             DO 10 I=1,NIONS
3056                MOLAL(I) = ZERO
3057 10          CONTINUE
3058             CALL CALCG1A
3059             SCASE = 'G2 ; SUBCASE 2'
3060          ELSE
3061             IF (W(5).GT. TINY) THEN
3062                SCASE = 'G2 ; SUBCASE 3'    ! MDRH (NH4CL, NA2SO4, NH42S4)
3063                CALL CALCMDRH (RH, DRMG2, DRNH4CL, CALCG1A, CALCG3A)
3064                SCASE = 'G2 ; SUBCASE 3'
3065             ENDIF
3066             IF (WATER.LE.TINY .AND. RH.GE.DRMG3) THEN
3067                SCASE = 'G2 ; SUBCASE 4'    ! MDRH (NA2SO4, NH42S4)
3068                CALL CALCMDRH (RH, DRMG3, DRNH42S4, CALCG1A, CALCG4)
3069                SCASE = 'G2 ; SUBCASE 4'
3070             ELSE
3071                WATER = TINY
3072                DO 20 I=1,NIONS
3073                   MOLAL(I) = ZERO
3074 20             CONTINUE
3075                CALL CALCG1A
3076                SCASE = 'G2 ; SUBCASE 2'
3077             ENDIF
3078          ENDIF
3079       ENDIF
3081       RETURN
3083 ! *** END OF SUBROUTINE CALCG2 ******************************************
3085       END SUBROUTINE CALCG2
3088 !=======================================================================
3090 ! *** ISORROPIA CODE
3091 ! *** SUBROUTINE CALCG2A
3092 ! *** CASE G2 ; SUBCASE 1
3094 !     THE MAIN CHARACTERISTICS OF THIS REGIME ARE:
3095 !     1. SULFATE POOR (SULRAT > 2.0) ; SODIUM POOR (SODRAT < 2.0)
3096 !     2. THERE IS BOTH A LIQUID & SOLID PHASE
3097 !     3. SOLIDS POSSIBLE : (NH4)2SO4, NH4CL, NA2SO4
3099 ! *** COPYRIGHT 1996-2006, UNIVERSITY OF MIAMI, CARNEGIE MELLON UNIVERSITY,
3100 ! *** GEORGIA INSTITUTE OF TECHNOLOGY
3101 ! *** WRITTEN BY ATHANASIOS NENES
3102 ! *** UPDATED BY CHRISTOS FOUNTOUKIS
3104 !=======================================================================
3106       SUBROUTINE CALCG2A
3107       USE ISRPIA
3108       USE CASEG
3109       IMPLICIT DOUBLE PRECISION (A-H, O-Z)
3110 !      INCLUDE 'ISRPIA.INC'
3112 !      DOUBLE PRECISION LAMDA
3113 !      COMMON /CASEG/ CHI1, CHI2, CHI3, CHI4, CHI5, CHI6, LAMDA,   &
3114 !                     PSI1, PSI2, PSI3, PSI4, PSI5, PSI6, PSI7,   &
3115 !                     A1,   A2,   A3,   A4,   A5,   A6,   A7
3117 ! *** SETUP PARAMETERS ************************************************
3119       CALAOU = .TRUE.
3120       CHI1   = 0.5*W(1)
3121       CHI2   = MAX (W(2)-CHI1, ZERO)
3122       CHI3   = ZERO
3123       CHI4   = MAX (W(3)-2.D0*CHI2, ZERO)
3124       CHI5   = W(4)
3125       CHI6   = W(5)
3127       PSI6LO = TINY
3128       PSI6HI = CHI6-TINY
3130       WATER  = TINY
3132 ! *** INITIAL VALUES FOR BISECTION ************************************
3134       X1 = PSI6LO
3135       Y1 = FUNCG2A (X1)
3136       IF (CHI6.LE.TINY) GOTO 50
3137 !CC      IF (ABS(Y1).LE.EPS .OR. CHI6.LE.TINY) GOTO 50
3138 !CC      IF (WATER .LE. TINY) GOTO 50               ! No water
3140 ! *** ROOT TRACKING ; FOR THE RANGE OF HI AND LO **********************
3142       DX = (PSI6HI-PSI6LO)/REAL(NDIV)
3143       DO 10 I=1,NDIV
3144          X2 = X1+DX
3145          Y2 = FUNCG2A (X2)
3146          IF (SIGN(1.d0,Y1)*SIGN(1.d0,Y2).LT.ZERO) GOTO 20  ! (Y1*Y2.LT.ZERO)
3147          X1 = X2
3148          Y1 = Y2
3149 10    CONTINUE
3151 ! *** NO SUBDIVISION WITH SOLUTION; IF ABS(Y2)<EPS SOLUTION IS ASSUMED
3153       IF (ABS(Y2) .GT. EPS) WATER = TINY
3154       GOTO 50
3156 ! *** PERFORM BISECTION ***********************************************
3158 20    DO 30 I=1,MAXIT
3159          X3 = 0.5*(X1+X2)
3160          Y3 = FUNCG2A (X3)
3161          IF (SIGN(1.d0,Y1)*SIGN(1.d0,Y3) .LE. ZERO) THEN  ! (Y1*Y3 .LE. ZERO)
3162             Y2    = Y3
3163             X2    = X3
3164          ELSE
3165             Y1    = Y3
3166             X1    = X3
3167          ENDIF
3168          IF (ABS(X2-X1) .LE. EPS*X1) GOTO 40
3169 30    CONTINUE
3170       CALL PUSHERR (0002, 'CALCG2A')    ! WARNING ERROR: NO CONVERGENCE
3172 ! *** CONVERGED ; RETURN **********************************************
3174 40    X3 = 0.5*(X1+X2)
3175       IF (X3.LE.TINY2) THEN   ! PRACTICALLY NO NITRATES, SO DRY SOLUTION
3176          WATER = TINY
3177       ELSE
3178          Y3 = FUNCG2A (X3)
3179       ENDIF
3181 ! *** CALCULATE HSO4 SPECIATION AND RETURN *******************************
3183 50    CONTINUE
3185 ! *** Na2SO4 DISSOLUTION
3187       IF (CHI1.GT.TINY .AND. WATER.GT.TINY) THEN        ! PSI1
3188          CALL POLY3 (PSI2, ZERO, -A1/4.D0, PSI1, ISLV)
3189          IF (ISLV.EQ.0) THEN
3190              PSI1 = MIN (PSI1, CHI1)
3191          ELSE
3192              PSI1 = ZERO
3193          ENDIF
3194       ELSE
3195          PSI1 = ZERO
3196       ENDIF
3197       MOLAL(2) = 2.0D0*PSI1               ! Na+  EFFECT
3198       MOLAL(5) = MOLAL(5) + PSI1          ! SO4  EFFECT
3199       CNA2SO4  = MAX(CHI1 - PSI1, ZERO)   ! NA2SO4(s) depletion
3201 ! *** HSO4 equilibrium
3203       IF (MOLAL(1).GT.TINY .AND. MOLAL(5).GT.TINY) THEN
3204          CALL CALCHS4 (MOLAL(1), MOLAL(5), ZERO, DELTA)
3205          MOLAL(1) = MOLAL(1) - DELTA     ! H+   AFFECT
3206          MOLAL(5) = MOLAL(5) - DELTA     ! SO4  AFFECT
3207          MOLAL(6) = DELTA                ! HSO4 AFFECT
3208       ENDIF
3210       RETURN
3212 ! *** END OF SUBROUTINE CALCG2A ******************************************
3214       END SUBROUTINE CALCG2A
3219 !=======================================================================
3221 ! *** ISORROPIA CODE
3222 ! *** SUBROUTINE FUNCG2A
3223 ! *** CASE G2
3225 !     THE MAIN CHARACTERISTICS OF THIS REGIME ARE:
3226 !     1. SULFATE POOR (SULRAT > 2.0) ; SODIUM POOR (SODRAT < 2.0)
3227 !     2. THERE IS BOTH A LIQUID & SOLID PHASE
3228 !     3. SOLIDS POSSIBLE : (NH4)2SO4, NH4CL, NA2SO4
3230 ! *** COPYRIGHT 1996-2006, UNIVERSITY OF MIAMI, CARNEGIE MELLON UNIVERSITY,
3231 ! *** GEORGIA INSTITUTE OF TECHNOLOGY
3232 ! *** WRITTEN BY ATHANASIOS NENES
3233 ! *** UPDATED BY CHRISTOS FOUNTOUKIS
3235 !=======================================================================
3237       DOUBLE PRECISION FUNCTION FUNCG2A (X)
3238       USE ISRPIA
3239       USE CASEG
3240       IMPLICIT DOUBLE PRECISION (A-H, O-Z)
3241 !      INCLUDE 'ISRPIA.INC'
3243 !      DOUBLE PRECISION LAMDA
3244 !      COMMON /CASEG/ CHI1, CHI2, CHI3, CHI4, CHI5, CHI6, LAMDA,   &
3245 !                     PSI1, PSI2, PSI3, PSI4, PSI5, PSI6, PSI7,   &
3246 !                     A1,   A2,   A3,   A4,   A5,   A6,   A7
3248 ! *** SETUP PARAMETERS ************************************************
3250       PSI6   = X
3251       PSI2   = CHI2
3252       PSI3   = ZERO
3253       FRST   = .TRUE.
3254       CALAIN = .TRUE.
3256 ! *** SOLVE EQUATIONS ; WITH ITERATIONS FOR ACTIVITY COEF. ************
3258       DO 10 I=1,NSWEEP
3260       A1  = XK5 *(WATER/GAMA(2))**3.0
3261       A2  = XK7 *(WATER/GAMA(4))**3.0
3262       A4  = (XK2/XKW)*R*TEMP*(GAMA(10)/GAMA(5))**2.0
3263       A5  = XK4 *R*TEMP*(WATER/GAMA(10))**2.0
3264       A6  = XK3 *R*TEMP*(WATER/GAMA(11))**2.0
3266       DENO = MAX(CHI6-PSI6-PSI3, ZERO)
3267       PSI5 = CHI5/((A6/A5)*(DENO/PSI6) + ONE)
3269       PSI4 = MIN(PSI5+PSI6,CHI4)
3271       IF (CHI2.GT.TINY .AND. WATER.GT.TINY) THEN
3272          CALL POLY3 (PSI4, PSI4*PSI4/4.D0, -A2/4.D0, PSI20, ISLV)
3273          IF (ISLV.EQ.0) PSI2 = MIN (PSI20, CHI2)
3274       ENDIF
3276 ! *** SAVE CONCENTRATIONS IN MOLAL ARRAY ******************************
3278       MOLAL (2) = ZERO                             ! NA
3279       MOLAL (3) = 2.0*PSI2 + PSI4                  ! NH4I
3280       MOLAL (4) = PSI6                             ! CLI
3281       MOLAL (5) = PSI2                             ! SO4I
3282       MOLAL (6) = ZERO                             ! HSO4
3283       MOLAL (7) = PSI5                             ! NO3I
3285 !CC      MOLAL (1) = MAX(CHI5 - PSI5, TINY)*A5/PSI5   ! HI
3286       SMIN      = 2.d0*MOLAL(5)+MOLAL(7)+MOLAL(4)-MOLAL(2)-MOLAL(3)
3287       CALL CALCPH (SMIN, HI, OHI)
3288       MOLAL (1) = HI
3290 ! *** CALCULATE GAS / SOLID SPECIES (LIQUID IN MOLAL ALREADY) *********
3292       GNH3      = MAX(CHI4 - PSI4, TINY)
3293       GHNO3     = MAX(CHI5 - PSI5, TINY)
3294       GHCL      = MAX(CHI6 - PSI6, TINY)
3296       CNH42S4   = MAX(CHI2 - PSI2, ZERO)
3297       CNH4NO3   = ZERO
3299 ! *** NH4Cl(s) calculations
3301       A3   = XK6 /(R*TEMP*R*TEMP)
3302       IF (GNH3*GHCL.GT.A3) THEN
3303          DELT = MIN(GNH3, GHCL)
3304          BB = -(GNH3+GHCL)
3305          CC = GNH3*GHCL-A3
3306          DD = BB*BB - 4.D0*CC
3307          PSI31 = 0.5D0*(-BB + SQRT(DD))
3308          PSI32 = 0.5D0*(-BB - SQRT(DD))
3309          IF (DELT-PSI31.GT.ZERO .AND. PSI31.GT.ZERO) THEN
3310             PSI3 = PSI31
3311          ELSEIF (DELT-PSI32.GT.ZERO .AND. PSI32.GT.ZERO) THEN
3312             PSI3 = PSI32
3313          ELSE
3314             PSI3 = ZERO
3315          ENDIF
3316       ELSE
3317          PSI3 = ZERO
3318       ENDIF
3320 ! *** CALCULATE GAS / SOLID SPECIES (LIQUID IN MOLAL ALREADY) *********
3322       GNH3    = MAX(GNH3 - PSI3, TINY)
3323       GHCL    = MAX(GHCL - PSI3, TINY)
3324       CNH4CL  = PSI3
3326 ! *** CALCULATE MOLALR ARRAY, WATER AND ACTIVITIES **********************
3328       CALL CALCMR
3330 ! *** CALCULATE ACTIVITIES OR TERMINATE INTERNAL LOOP *****************
3332       IF (FRST.AND.CALAOU .OR. .NOT.FRST.AND.CALAIN) THEN
3333          CALL CALCACT
3334       ELSE
3335          GOTO 20
3336       ENDIF
3337 10    CONTINUE
3339 ! *** CALCULATE FUNCTION VALUE FOR OUTER LOOP ***************************
3341 20    IF (CHI4.LE.TINY) THEN
3342          FUNCG2A = MOLAL(1)*MOLAL(4)/GHCL/A6 - ONE
3343       ELSE
3344          FUNCG2A = MOLAL(3)*MOLAL(4)/GHCL/GNH3/A6/A4 - ONE
3345       ENDIF
3347       RETURN
3349 ! *** END OF FUNCTION FUNCG2A *******************************************
3351       END FUNCTION FUNCG2A    
3353 !=======================================================================
3355 ! *** ISORROPIA CODE
3356 ! *** SUBROUTINE CALCG1
3357 ! *** CASE G1
3359 !     THE MAIN CHARACTERISTICS OF THIS REGIME ARE:
3360 !     1. SULFATE POOR (SULRAT > 2.0) ; SODIUM POOR (SODRAT < 2.0)
3361 !     2. SOLID AEROSOL ONLY
3362 !     3. SOLIDS POSSIBLE : (NH4)2SO4, NH4NO3, NH4CL, NA2SO4
3364 !     THERE ARE TWO POSSIBLE REGIMES HERE, DEPENDING ON RELATIVE HUMIDITY:
3365 !     1. WHEN RH >= MDRH ; LIQUID PHASE POSSIBLE (MDRH REGION)
3366 !     2. WHEN RH < MDRH  ; ONLY SOLID PHASE POSSIBLE (SUBROUTINE CALCG1A)
3368 ! *** COPYRIGHT 1996-2006, UNIVERSITY OF MIAMI, CARNEGIE MELLON UNIVERSITY,
3369 ! *** GEORGIA INSTITUTE OF TECHNOLOGY
3370 ! *** WRITTEN BY ATHANASIOS NENES
3371 ! *** UPDATED BY CHRISTOS FOUNTOUKIS
3373 !=======================================================================
3375       SUBROUTINE CALCG1
3376       USE ISRPIA
3377       IMPLICIT DOUBLE PRECISION (A-H, O-Z)
3378 !      INCLUDE 'ISRPIA.INC'
3379       EXTERNAL CALCG1A, CALCG2A
3381 ! *** REGIME DEPENDS UPON THE AMBIENT RELATIVE HUMIDITY *****************
3383       IF (RH.LT.DRMG1) THEN
3384          SCASE = 'G1 ; SUBCASE 1'
3385          CALL CALCG1A              ! SOLID PHASE ONLY POSSIBLE
3386          SCASE = 'G1 ; SUBCASE 1'
3387       ELSE
3388          SCASE = 'G1 ; SUBCASE 2'  ! LIQUID & SOLID PHASE POSSIBLE
3389          CALL CALCMDRH (RH, DRMG1, DRNH4NO3, CALCG1A, CALCG2A)
3390          SCASE = 'G1 ; SUBCASE 2'
3391       ENDIF
3393       RETURN
3395 ! *** END OF SUBROUTINE CALCG1 ******************************************
3397       END SUBROUTINE CALCG1
3400 !=======================================================================
3402 ! *** ISORROPIA CODE
3403 ! *** SUBROUTINE CALCG1A
3404 ! *** CASE G1 ; SUBCASE 1
3406 !     THE MAIN CHARACTERISTICS OF THIS REGIME ARE:
3407 !     1. SULFATE POOR (SULRAT > 2.0)
3408 !     2. SOLID AEROSOL ONLY
3409 !     3. SOLIDS POSSIBLE : (NH4)2SO4, NH4NO3
3411 !     SOLID (NH4)2SO4 IS CALCULATED FROM THE SULFATES, WHILE NH4NO3
3412 !     IS CALCULATED FROM NH3-HNO3 EQUILIBRIUM. 'ZE' IS THE AMOUNT OF
3413 !     NH4NO3 THAT VOLATIZES WHEN ALL POSSILBE NH4NO3 IS INITIALLY IN
3414 !     THE SOLID PHASE.
3416 ! *** COPYRIGHT 1996-2006, UNIVERSITY OF MIAMI, CARNEGIE MELLON UNIVERSITY,
3417 ! *** GEORGIA INSTITUTE OF TECHNOLOGY
3418 ! *** WRITTEN BY ATHANASIOS NENES
3419 ! *** UPDATED BY CHRISTOS FOUNTOUKIS
3421 !=======================================================================
3423       SUBROUTINE CALCG1A
3424       USE ISRPIA
3425       IMPLICIT DOUBLE PRECISION (A-H, O-Z)
3426 !      INCLUDE 'ISRPIA.INC'
3427       DOUBLE PRECISION LAMDA, LAMDA1, LAMDA2, KAPA, KAPA1, KAPA2
3429 ! *** CALCULATE NON VOLATILE SOLIDS ***********************************
3431       CNA2SO4 = 0.5*W(1)
3432       CNH42S4 = W(2) - CNA2SO4
3434 ! *** CALCULATE VOLATILE SPECIES **************************************
3436       ALF     = W(3) - 2.0*CNH42S4
3437       BET     = W(5)
3438       GAM     = W(4)
3440       RTSQ    = R*TEMP*R*TEMP
3441       A1      = XK6/RTSQ
3442       A2      = XK10/RTSQ
3444       THETA1  = GAM - BET*(A2/A1)
3445       THETA2  = A2/A1
3447 ! QUADRATIC EQUATION SOLUTION
3449       BB      = (THETA1-ALF-BET*(ONE+THETA2))/(ONE+THETA2)
3450       CC      = (ALF*BET-A1-BET*THETA1)/(ONE+THETA2)
3451       DD      = BB*BB - 4.0D0*CC
3452       IF (DD.LT.ZERO) GOTO 100   ! Solve each reaction seperately
3454 ! TWO ROOTS FOR KAPA, CHECK AND SEE IF ANY VALID
3456       SQDD    = SQRT(DD)
3457       KAPA1   = 0.5D0*(-BB+SQDD)
3458       KAPA2   = 0.5D0*(-BB-SQDD)
3459       LAMDA1  = THETA1 + THETA2*KAPA1
3460       LAMDA2  = THETA1 + THETA2*KAPA2
3462       IF (KAPA1.GE.ZERO .AND. LAMDA1.GE.ZERO) THEN
3463          IF (ALF-KAPA1-LAMDA1.GE.ZERO .AND.   &
3464              BET-KAPA1.GE.ZERO .AND. GAM-LAMDA1.GE.ZERO) THEN
3465              KAPA = KAPA1
3466              LAMDA= LAMDA1
3467              GOTO 200
3468          ENDIF
3469       ENDIF
3471       IF (KAPA2.GE.ZERO .AND. LAMDA2.GE.ZERO) THEN
3472          IF (ALF-KAPA2-LAMDA2.GE.ZERO .AND.   &
3473              BET-KAPA2.GE.ZERO .AND. GAM-LAMDA2.GE.ZERO) THEN
3474              KAPA = KAPA2
3475              LAMDA= LAMDA2
3476              GOTO 200
3477          ENDIF
3478       ENDIF
3480 ! SEPERATE SOLUTION OF NH4CL & NH4NO3 EQUILIBRIA
3482 100   KAPA  = ZERO
3483       LAMDA = ZERO
3484       DD1   = (ALF+BET)*(ALF+BET) - 4.0D0*(ALF*BET-A1)
3485       DD2   = (ALF+GAM)*(ALF+GAM) - 4.0D0*(ALF*GAM-A2)
3487 ! NH4CL EQUILIBRIUM
3489       IF (DD1.GE.ZERO) THEN
3490          SQDD1 = SQRT(DD1)
3491          KAPA1 = 0.5D0*(ALF+BET + SQDD1)
3492          KAPA2 = 0.5D0*(ALF+BET - SQDD1)
3494          IF (KAPA1.GE.ZERO .AND. KAPA1.LE.MIN(ALF,BET)) THEN
3495             KAPA = KAPA1
3496          ELSE IF (KAPA2.GE.ZERO .AND. KAPA2.LE.MIN(ALF,BET)) THEN
3497             KAPA = KAPA2
3498          ELSE
3499             KAPA = ZERO
3500          ENDIF
3501       ENDIF
3503 ! NH4NO3 EQUILIBRIUM
3505       IF (DD2.GE.ZERO) THEN
3506          SQDD2 = SQRT(DD2)
3507          LAMDA1= 0.5D0*(ALF+GAM + SQDD2)
3508          LAMDA2= 0.5D0*(ALF+GAM - SQDD2)
3510          IF (LAMDA1.GE.ZERO .AND. LAMDA1.LE.MIN(ALF,GAM)) THEN
3511             LAMDA = LAMDA1
3512          ELSE IF (LAMDA2.GE.ZERO .AND. LAMDA2.LE.MIN(ALF,GAM)) THEN
3513             LAMDA = LAMDA2
3514          ELSE
3515             LAMDA = ZERO
3516          ENDIF
3517       ENDIF
3519 ! IF BOTH KAPA, LAMDA ARE > 0, THEN APPLY EXISTANCE CRITERION
3521       IF (KAPA.GT.ZERO .AND. LAMDA.GT.ZERO) THEN
3522          IF (BET .LT. LAMDA/THETA1) THEN
3523             KAPA = ZERO
3524          ELSE
3525             LAMDA= ZERO
3526          ENDIF
3527       ENDIF
3529 ! *** CALCULATE COMPOSITION OF VOLATILE SPECIES ***********************
3531 200   CONTINUE
3532       CNH4NO3 = LAMDA
3533       CNH4CL  = KAPA
3535       GNH3    = MAX(ALF - KAPA - LAMDA, ZERO)
3536       GHNO3   = MAX(GAM - LAMDA, ZERO)
3537       GHCL    = MAX(BET - KAPA, ZERO)
3539       RETURN
3541 ! *** END OF SUBROUTINE CALCG1A *****************************************
3543       END SUBROUTINE CALCG1A
3544 !=======================================================================
3546 ! *** ISORROPIA CODE
3547 ! *** SUBROUTINE CALCH6
3548 ! *** CASE H6
3550 !     THE MAIN CHARACTERISTICS OF THIS REGIME ARE:
3551 !     1. SULFATE POOR (SULRAT > 2.0) ; SODIUM RICH (SODRAT >= 2.0)
3552 !     2. THERE IS BOTH A LIQUID & SOLID PHASE
3553 !     3. SOLIDS POSSIBLE : (NH4)2SO4, NH4CL, NA2SO4
3555 ! *** COPYRIGHT 1996-2006, UNIVERSITY OF MIAMI, CARNEGIE MELLON UNIVERSITY,
3556 ! *** GEORGIA INSTITUTE OF TECHNOLOGY
3557 ! *** WRITTEN BY ATHANASIOS NENES
3558 ! *** UPDATED BY CHRISTOS FOUNTOUKIS
3560 !=======================================================================
3562       SUBROUTINE CALCH6
3563       USE ISRPIA
3564       USE SOLUT
3565       IMPLICIT DOUBLE PRECISION (A-H, O-Z)
3566 !      INCLUDE 'ISRPIA.INC'
3568 !      COMMON /SOLUT/ CHI1, CHI2, CHI3, CHI4, CHI5, CHI6, CHI7, CHI8,   &
3569 !                     PSI1, PSI2, PSI3, PSI4, PSI5, PSI6, PSI7, PSI8,   &
3570 !                     A1,   A2,   A3,   A4,   A5,   A6,   A7,   A8
3572 ! *** SETUP PARAMETERS ************************************************
3574       CALAOU = .TRUE.
3575       CHI1   = W(2)                                ! CNA2SO4
3576       CHI2   = ZERO                                ! CNH42S4
3577       CHI3   = ZERO                                ! CNH4CL
3578       FRNA   = MAX (W(1)-2.D0*CHI1, ZERO)
3579       CHI8   = MIN (FRNA, W(4))                    ! CNANO3
3580       CHI4   = W(3)                                ! NH3(g)
3581       CHI5   = MAX (W(4)-CHI8, ZERO)               ! HNO3(g)
3582       CHI7   = MIN (MAX(FRNA-CHI8, ZERO), W(5))    ! CNACL
3583       CHI6   = MAX (W(5)-CHI7, ZERO)               ! HCL(g)
3585       PSI6LO = TINY
3586       PSI6HI = CHI6-TINY    ! MIN(CHI6-TINY, CHI4)
3588 ! *** INITIAL VALUES FOR BISECTION ************************************
3590       X1 = PSI6LO
3591       Y1 = FUNCH6A (X1)
3592       IF (ABS(Y1).LE.EPS .OR. CHI6.LE.TINY) GOTO 50
3594 ! *** ROOT TRACKING ; FOR THE RANGE OF HI AND LO **********************
3596       DX = (PSI6HI-PSI6LO)/REAL(NDIV)
3597       DO 10 I=1,NDIV
3598          X2 = X1+DX
3599          Y2 = FUNCH6A (X2)
3600          IF (SIGN(1.d0,Y1)*SIGN(1.d0,Y2).LT.ZERO) GOTO 20  ! (Y1*Y2.LT.ZERO)
3601          X1 = X2
3602          Y1 = Y2
3603 10    CONTINUE
3605 ! *** NO SUBDIVISION WITH SOLUTION; IF ABS(Y2)<EPS SOLUTION IS ASSUMED
3607       IF (ABS(Y2) .GT. EPS) Y2 = FUNCH6A (PSI6LO)
3608       GOTO 50
3610 ! *** PERFORM BISECTION ***********************************************
3612 20    DO 30 I=1,MAXIT
3613          X3 = 0.5*(X1+X2)
3614          Y3 = FUNCH6A (X3)
3615          IF (SIGN(1.d0,Y1)*SIGN(1.d0,Y3) .LE. ZERO) THEN  ! (Y1*Y3 .LE. ZERO)
3616             Y2    = Y3
3617             X2    = X3
3618          ELSE
3619             Y1    = Y3
3620             X1    = X3
3621          ENDIF
3622          IF (ABS(X2-X1) .LE. EPS*X1) GOTO 40
3623 30    CONTINUE
3624       CALL PUSHERR (0002, 'CALCH6')    ! WARNING ERROR: NO CONVERGENCE
3626 ! *** CONVERGED ; RETURN **********************************************
3628 40    X3 = 0.5*(X1+X2)
3629       Y3 = FUNCH6A (X3)
3631 ! *** CALCULATE HSO4 SPECIATION AND RETURN *******************************
3633 50    CONTINUE
3634       IF (MOLAL(1).GT.TINY .AND. MOLAL(5).GT.TINY) THEN
3635          CALL CALCHS4 (MOLAL(1), MOLAL(5), ZERO, DELTA)
3636          MOLAL(1) = MOLAL(1) - DELTA                     ! H+   EFFECT
3637          MOLAL(5) = MOLAL(5) - DELTA                     ! SO4  EFFECT
3638          MOLAL(6) = DELTA                                ! HSO4 EFFECT
3639       ENDIF
3641       RETURN
3643 ! *** END OF SUBROUTINE CALCH6 ******************************************
3645       END SUBROUTINE CALCH6
3650 !=======================================================================
3652 ! *** ISORROPIA CODE
3653 ! *** SUBROUTINE FUNCH6A
3654 ! *** CASE H6
3656 !     THE MAIN CHARACTERISTICS OF THIS REGIME ARE:
3657 !     1. SULFATE POOR (SULRAT > 2.0) ; SODIUM RICH (SODRAT >= 2.0)
3658 !     2. THERE IS BOTH A LIQUID & SOLID PHASE
3659 !     3. SOLIDS POSSIBLE : (NH4)2SO4, NH4CL, NA2SO4
3661 ! *** COPYRIGHT 1996-2006, UNIVERSITY OF MIAMI, CARNEGIE MELLON UNIVERSITY,
3662 ! *** GEORGIA INSTITUTE OF TECHNOLOGY
3663 ! *** WRITTEN BY ATHANASIOS NENES
3664 ! *** UPDATED BY CHRISTOS FOUNTOUKIS
3666 !=======================================================================
3668       DOUBLE PRECISION FUNCTION FUNCH6A (X)
3669       USE ISRPIA
3670       USE SOLUT
3671       IMPLICIT DOUBLE PRECISION (A-H, O-Z)
3672 !      INCLUDE 'ISRPIA.INC'
3674 !      COMMON /SOLUT/ CHI1, CHI2, CHI3, CHI4, CHI5, CHI6, CHI7, CHI8,   &
3675 !                     PSI1, PSI2, PSI3, PSI4, PSI5, PSI6, PSI7, PSI8,   &
3676 !                     A1,   A2,   A3,   A4,   A5,   A6,   A7,   A8
3678 ! *** SETUP PARAMETERS ************************************************
3680       PSI6   = X
3681       PSI1   = CHI1
3682       PSI2   = ZERO
3683       PSI3   = ZERO
3684       PSI7   = CHI7
3685       PSI8   = CHI8
3686       FRST   = .TRUE.
3687       CALAIN = .TRUE.
3689 ! *** SOLVE EQUATIONS ; WITH ITERATIONS FOR ACTIVITY COEF. ************
3691       DO 10 I=1,NSWEEP
3693       A1  = XK5 *(WATER/GAMA(2))**3.0
3694       A4  = (XK2/XKW)*R*TEMP*(GAMA(10)/GAMA(5))**2.0
3695       A5  = XK4 *R*TEMP*(WATER/GAMA(10))**2.0
3696       A6  = XK3 *R*TEMP*(WATER/GAMA(11))**2.0
3697       A7  = XK8 *(WATER/GAMA(1))**2.0
3698       A8  = XK9 *(WATER/GAMA(3))**2.0
3699       A9  = XK1*WATER/GAMA(7)*(GAMA(8)/GAMA(7))**2.
3701 !  CALCULATE DISSOCIATION QUANTITIES
3703       PSI5 = CHI5*(PSI6+PSI7) - A6/A5*PSI8*(CHI6-PSI6-PSI3)
3704       PSI5 = PSI5/(A6/A5*(CHI6-PSI6-PSI3) + PSI6 + PSI7)
3705       PSI5 = MAX(PSI5, TINY)
3707       IF (W(3).GT.TINY .AND. WATER.GT.TINY) THEN  ! First try 3rd order soln
3708          BB   =-(CHI4 + PSI6 + PSI5 + 1.d0/A4)
3709          CC   = CHI4*(PSI5+PSI6)
3710          DD   = BB*BB-4.d0*CC
3711          PSI4 =0.5d0*(-BB - SQRT(DD))
3712          PSI4 = MIN(PSI4,CHI4)
3713       ELSE
3714          PSI4 = TINY
3715       ENDIF
3717 ! *** CALCULATE SPECIATION ********************************************
3719       MOLAL (2) = PSI8 + PSI7 + 2.D0*PSI1               ! NAI
3720       MOLAL (3) = PSI4                                  ! NH4I
3721       MOLAL (4) = PSI6 + PSI7                           ! CLI
3722       MOLAL (5) = PSI2 + PSI1                           ! SO4I
3723       MOLAL (6) = ZERO                                  ! HSO4I
3724       MOLAL (7) = PSI5 + PSI8                           ! NO3I
3726       SMIN      = 2.d0*MOLAL(5)+MOLAL(7)+MOLAL(4)-MOLAL(2)-MOLAL(3)
3727       CALL CALCPH (SMIN, HI, OHI)
3728       MOLAL (1) = HI
3730       GNH3      = MAX(CHI4 - PSI4, TINY)
3731       GHNO3     = MAX(CHI5 - PSI5, TINY)
3732       GHCL      = MAX(CHI6 - PSI6, TINY)
3734       CNH42S4   = ZERO
3735       CNH4NO3   = ZERO
3736       CNACL     = MAX(CHI7 - PSI7, ZERO)
3737       CNANO3    = MAX(CHI8 - PSI8, ZERO)
3738       CNA2SO4   = MAX(CHI1 - PSI1, ZERO)
3740       CALL CALCMR                                    ! Water content
3742 ! *** CALCULATE ACTIVITIES OR TERMINATE INTERNAL LOOP *****************
3744       IF (FRST.AND.CALAOU .OR. .NOT.FRST.AND.CALAIN) THEN
3745          CALL CALCACT
3746       ELSE
3747          GOTO 20
3748       ENDIF
3749 10    CONTINUE
3751 ! *** CALCULATE FUNCTION VALUE FOR OUTER LOOP ***************************
3753 20    FUNCH6A = MOLAL(3)*MOLAL(4)/GHCL/GNH3/A6/A4 - ONE
3755       RETURN
3757 ! *** END OF FUNCTION FUNCH6A *******************************************
3759       END FUNCTION FUNCH6A    
3761 !=======================================================================
3763 ! *** ISORROPIA CODE
3764 ! *** SUBROUTINE CALCH5
3765 ! *** CASE H5
3767 !     THE MAIN CHARACTERISTICS OF THIS REGIME ARE:
3768 !     1. SULFATE POOR (SULRAT > 2.0) ; SODIUM RICH (SODRAT >= 2.0)
3769 !     2. THERE IS BOTH A LIQUID & SOLID PHASE
3770 !     3. SOLIDS POSSIBLE : (NH4)2SO4, NH4CL, NA2SO4
3772 ! *** COPYRIGHT 1996-2006, UNIVERSITY OF MIAMI, CARNEGIE MELLON UNIVERSITY,
3773 ! *** GEORGIA INSTITUTE OF TECHNOLOGY
3774 ! *** WRITTEN BY ATHANASIOS NENES
3775 ! *** UPDATED BY CHRISTOS FOUNTOUKIS
3777 !=======================================================================
3779       SUBROUTINE CALCH5
3780       USE ISRPIA
3781       USE SOLUT
3782       IMPLICIT DOUBLE PRECISION (A-H, O-Z)
3783 !      INCLUDE 'ISRPIA.INC'
3785 !      COMMON /SOLUT/ CHI1, CHI2, CHI3, CHI4, CHI5, CHI6, CHI7, CHI8,   &
3786 !                     PSI1, PSI2, PSI3, PSI4, PSI5, PSI6, PSI7, PSI8,   &
3787 !                     A1,   A2,   A3,   A4,   A5,   A6,   A7,   A8
3789 ! *** REGIME DEPENDS ON THE EXISTANCE OF NITRATES ***********************
3791       IF (W(4).LE.TINY .AND. W(5).LE.TINY) THEN
3792          SCASE = 'H5'
3793          CALL CALCH1A
3794          SCASE = 'H5'
3795          RETURN
3796       ENDIF
3798 ! *** SETUP PARAMETERS ************************************************
3800       CALAOU = .TRUE.
3801       CHI1   = W(2)                                ! CNA2SO4
3802       CHI2   = ZERO                                ! CNH42S4
3803       CHI3   = ZERO                                ! CNH4CL
3804       FRNA   = MAX (W(1)-2.D0*CHI1, ZERO)
3805       CHI8   = MIN (FRNA, W(4))                    ! CNANO3
3806       CHI4   = W(3)                                ! NH3(g)
3807       CHI5   = MAX (W(4)-CHI8, ZERO)               ! HNO3(g)
3808       CHI7   = MIN (MAX(FRNA-CHI8, ZERO), W(5))    ! CNACL
3809       CHI6   = MAX (W(5)-CHI7, ZERO)               ! HCL(g)
3811       PSI6LO = TINY
3812       PSI6HI = CHI6-TINY    ! MIN(CHI6-TINY, CHI4)
3814 ! *** INITIAL VALUES FOR BISECTION ************************************
3816       X1 = PSI6LO
3817       Y1 = FUNCH5A (X1)
3818       IF (ABS(Y1).LE.EPS .OR. CHI6.LE.TINY) GOTO 50
3820 ! *** ROOT TRACKING ; FOR THE RANGE OF HI AND LO **********************
3822       DX = (PSI6HI-PSI6LO)/REAL(NDIV)
3823       DO 10 I=1,NDIV
3824          X2 = X1+DX
3825          Y2 = FUNCH5A (X2)
3826          IF (SIGN(1.d0,Y1)*SIGN(1.d0,Y2).LT.ZERO) GOTO 20  ! (Y1*Y2.LT.ZERO)
3827          X1 = X2
3828          Y1 = Y2
3829 10    CONTINUE
3831 ! *** NO SUBDIVISION WITH SOLUTION; IF ABS(Y2)<EPS SOLUTION IS ASSUMED
3833       IF (ABS(Y2) .GT. EPS) Y2 = FUNCH5A (PSI6LO)
3834       GOTO 50
3836 ! *** PERFORM BISECTION ***********************************************
3838 20    DO 30 I=1,MAXIT
3839          X3 = 0.5*(X1+X2)
3840          Y3 = FUNCH5A (X3)
3841          IF (SIGN(1.d0,Y1)*SIGN(1.d0,Y3) .LE. ZERO) THEN  ! (Y1*Y3 .LE. ZERO)
3842             Y2    = Y3
3843             X2    = X3
3844          ELSE
3845             Y1    = Y3
3846             X1    = X3
3847          ENDIF
3848          IF (ABS(X2-X1) .LE. EPS*X1) GOTO 40
3849 30    CONTINUE
3850       CALL PUSHERR (0002, 'CALCH5')    ! WARNING ERROR: NO CONVERGENCE
3852 ! *** CONVERGED ; RETURN **********************************************
3854 40    X3 = 0.5*(X1+X2)
3855       Y3 = FUNCH5A (X3)
3857 ! *** CALCULATE HSO4 SPECIATION AND RETURN *******************************
3859 50    CONTINUE
3860       IF (MOLAL(1).GT.TINY .AND. MOLAL(5).GT.TINY) THEN
3861          CALL CALCHS4 (MOLAL(1), MOLAL(5), ZERO, DELTA)
3862          MOLAL(1) = MOLAL(1) - DELTA                     ! H+   EFECT
3863          MOLAL(5) = MOLAL(5) - DELTA                     ! SO4  EFFECT
3864          MOLAL(6) = DELTA                                ! HSO4 EFFECT
3865       ENDIF
3867       RETURN
3869 ! *** END OF SUBROUTINE CALCH5 ******************************************
3871       END SUBROUTINE CALCH5
3876 !=======================================================================
3878 ! *** ISORROPIA CODE
3879 ! *** SUBROUTINE FUNCH5A
3880 ! *** CASE H5
3882 !     THE MAIN CHARACTERISTICS OF THIS REGIME ARE:
3883 !     1. SULFATE POOR (SULRAT > 2.0) ; SODIUM RICH (SODRAT >= 2.0)
3884 !     2. THERE IS BOTH A LIQUID & SOLID PHASE
3885 !     3. SOLIDS POSSIBLE : NONE
3887 ! *** COPYRIGHT 1996-2006, UNIVERSITY OF MIAMI, CARNEGIE MELLON UNIVERSITY,
3888 ! *** GEORGIA INSTITUTE OF TECHNOLOGY
3889 ! *** WRITTEN BY ATHANASIOS NENES
3890 ! *** UPDATED BY CHRISTOS FOUNTOUKIS
3892 !=======================================================================
3894       DOUBLE PRECISION FUNCTION FUNCH5A (X)
3895       USE ISRPIA
3896       USE SOLUT
3897       IMPLICIT DOUBLE PRECISION (A-H, O-Z)
3898 !      INCLUDE 'ISRPIA.INC'
3900 !      COMMON /SOLUT/ CHI1, CHI2, CHI3, CHI4, CHI5, CHI6, CHI7, CHI8,   &
3901 !                     PSI1, PSI2, PSI3, PSI4, PSI5, PSI6, PSI7, PSI8,   &
3902 !                     A1,   A2,   A3,   A4,   A5,   A6,   A7,   A8
3904 ! *** SETUP PARAMETERS ************************************************
3906       PSI6   = X
3907       PSI1   = CHI1
3908       PSI2   = ZERO
3909       PSI3   = ZERO
3910       PSI7   = CHI7
3911       PSI8   = CHI8
3912       FRST   = .TRUE.
3913       CALAIN = .TRUE.
3915 ! *** SOLVE EQUATIONS ; WITH ITERATIONS FOR ACTIVITY COEF. ************
3917       DO 10 I=1,NSWEEP
3919       A1  = XK5 *(WATER/GAMA(2))**3.0
3920       A4  = (XK2/XKW)*R*TEMP*(GAMA(10)/GAMA(5))**2.0
3921       A5  = XK4 *R*TEMP*(WATER/GAMA(10))**2.0
3922       A6  = XK3 *R*TEMP*(WATER/GAMA(11))**2.0
3923       A7  = XK8 *(WATER/GAMA(1))**2.0
3924       A8  = XK9 *(WATER/GAMA(3))**2.0
3925       A9  = XK1*WATER/GAMA(7)*(GAMA(8)/GAMA(7))**2.
3927 !  CALCULATE DISSOCIATION QUANTITIES
3929       PSI5 = CHI5*(PSI6+PSI7) - A6/A5*PSI8*(CHI6-PSI6-PSI3)
3930       PSI5 = PSI5/(A6/A5*(CHI6-PSI6-PSI3) + PSI6 + PSI7)
3931       PSI5 = MAX(PSI5, TINY)
3933       IF (W(3).GT.TINY .AND. WATER.GT.TINY) THEN  ! First try 3rd order soln
3934          BB   =-(CHI4 + PSI6 + PSI5 + 1.d0/A4)
3935          CC   = CHI4*(PSI5+PSI6)
3936          DD   = BB*BB-4.d0*CC
3937          PSI4 =0.5d0*(-BB - SQRT(DD))
3938          PSI4 = MIN(PSI4,CHI4)
3939       ELSE
3940          PSI4 = TINY
3941       ENDIF
3943       IF (CHI1.GT.TINY .AND. WATER.GT.TINY) THEN     ! NA2SO4 DISSOLUTION
3944          AA = PSI7+PSI8
3945          BB = AA*AA
3946          CC =-A1/4.D0
3947          CALL POLY3 (AA, BB, CC, PSI1, ISLV)
3948          IF (ISLV.EQ.0) THEN
3949              PSI1 = MIN (PSI1, CHI1)
3950          ELSE
3951              PSI1 = ZERO
3952          ENDIF
3953       ENDIF
3955 ! *** CALCULATE SPECIATION ********************************************
3957       MOLAL (2) = PSI8 + PSI7 + 2.D0*PSI1                ! NAI
3958       MOLAL (3) = PSI4                                   ! NH4I
3959       MOLAL (4) = PSI6 + PSI7                            ! CLI
3960       MOLAL (5) = PSI2 + PSI1                            ! SO4I
3961       MOLAL (6) = ZERO
3962       MOLAL (7) = PSI5 + PSI8                            ! NO3I
3964       SMIN      = 2.d0*MOLAL(5)+MOLAL(7)+MOLAL(4)-MOLAL(2)-MOLAL(3)
3965       CALL CALCPH (SMIN, HI, OHI)
3966       MOLAL (1) = HI
3968       GNH3      = MAX(CHI4 - PSI4, TINY)
3969       GHNO3     = MAX(CHI5 - PSI5, TINY)
3970       GHCL      = MAX(CHI6 - PSI6, TINY)
3972       CNH42S4   = ZERO
3973       CNH4NO3   = ZERO
3974       CNACL     = MAX(CHI7 - PSI7, ZERO)
3975       CNANO3    = MAX(CHI8 - PSI8, ZERO)
3976       CNA2SO4   = MAX(CHI1 - PSI1, ZERO)
3978       CALL CALCMR                               ! Water content
3980 ! *** CALCULATE ACTIVITIES OR TERMINATE INTERNAL LOOP *****************
3982       IF (FRST.AND.CALAOU .OR. .NOT.FRST.AND.CALAIN) THEN
3983          CALL CALCACT
3984       ELSE
3985          GOTO 20
3986       ENDIF
3987 10    CONTINUE
3989 ! *** CALCULATE FUNCTION VALUE FOR OUTER LOOP ***************************
3991 20    FUNCH5A = MOLAL(3)*MOLAL(4)/GHCL/GNH3/A6/A4 - ONE
3993       RETURN
3995 ! *** END OF FUNCTION FUNCH5A *******************************************
3997       END FUNCTION FUNCH5A    
3999 !=======================================================================
4001 ! *** ISORROPIA CODE
4002 ! *** SUBROUTINE CALCH4
4003 ! *** CASE H4
4005 !     THE MAIN CHARACTERISTICS OF THIS REGIME ARE:
4006 !     1. SULFATE POOR (SULRAT > 2.0) ; SODIUM RICH (SODRAT >= 2.0)
4007 !     2. THERE IS BOTH A LIQUID & SOLID PHASE
4008 !     3. SOLIDS POSSIBLE : NA2SO4
4010 ! *** COPYRIGHT 1996-2006, UNIVERSITY OF MIAMI, CARNEGIE MELLON UNIVERSITY,
4011 ! *** GEORGIA INSTITUTE OF TECHNOLOGY
4012 ! *** WRITTEN BY ATHANASIOS NENES
4013 ! *** UPDATED BY CHRISTOS FOUNTOUKIS
4015 !=======================================================================
4017       SUBROUTINE CALCH4
4018       USE ISRPIA
4019       USE SOLUT
4020       IMPLICIT DOUBLE PRECISION (A-H, O-Z)
4021 !      INCLUDE 'ISRPIA.INC'
4023 !      COMMON /SOLUT/ CHI1, CHI2, CHI3, CHI4, CHI5, CHI6, CHI7, CHI8,   &
4024 !                     PSI1, PSI2, PSI3, PSI4, PSI5, PSI6, PSI7, PSI8,   &
4025 !                     A1,   A2,   A3,   A4,   A5,   A6,   A7,   A8
4027 ! *** REGIME DEPENDS ON THE EXISTANCE OF NITRATES ***********************
4029       IF (W(4).LE.TINY .AND. W(5).LE.TINY) THEN
4030          SCASE = 'H4'
4031          CALL CALCH1A
4032          SCASE = 'H4'
4033          RETURN
4034       ENDIF
4036 ! *** SETUP PARAMETERS ************************************************
4038       CALAOU = .TRUE.
4039       CHI1   = W(2)                                ! CNA2SO4
4040       CHI2   = ZERO                                ! CNH42S4
4041       CHI3   = ZERO                                ! CNH4CL
4042       FRNA   = MAX (W(1)-2.D0*CHI1, ZERO)
4043       CHI8   = MIN (FRNA, W(4))                    ! CNANO3
4044       CHI4   = W(3)                                ! NH3(g)
4045       CHI5   = MAX (W(4)-CHI8, ZERO)               ! HNO3(g)
4046       CHI7   = MIN (MAX(FRNA-CHI8, ZERO), W(5))    ! CNACL
4047       CHI6   = MAX (W(5)-CHI7, ZERO)               ! HCL(g)
4049       PSI6LO = TINY
4050       PSI6HI = CHI6-TINY    ! MIN(CHI6-TINY, CHI4)
4052 ! *** INITIAL VALUES FOR BISECTION ************************************
4054       X1 = PSI6LO
4055       Y1 = FUNCH4A (X1)
4056       IF (ABS(Y1).LE.EPS .OR. CHI6.LE.TINY) GOTO 50
4058 ! *** ROOT TRACKING ; FOR THE RANGE OF HI AND LO **********************
4060       DX = (PSI6HI-PSI6LO)/REAL(NDIV)
4061       DO 10 I=1,NDIV
4062          X2 = X1+DX
4063          Y2 = FUNCH4A (X2)
4064          IF (SIGN(1.d0,Y1)*SIGN(1.d0,Y2).LT.ZERO) GOTO 20  ! (Y1*Y2.LT.ZERO)
4065          X1 = X2
4066          Y1 = Y2
4067 10    CONTINUE
4069 ! *** NO SUBDIVISION WITH SOLUTION; IF ABS(Y2)<EPS SOLUTION IS ASSUMED
4071       IF (ABS(Y2) .GT. EPS) Y2 = FUNCH4A (PSI6LO)
4072       GOTO 50
4074 ! *** PERFORM BISECTION ***********************************************
4076 20    DO 30 I=1,MAXIT
4077          X3 = 0.5*(X1+X2)
4078          Y3 = FUNCH4A (X3)
4079          IF (SIGN(1.d0,Y1)*SIGN(1.d0,Y3) .LE. ZERO) THEN  ! (Y1*Y3 .LE. ZERO)
4080             Y2    = Y3
4081             X2    = X3
4082          ELSE
4083             Y1    = Y3
4084             X1    = X3
4085          ENDIF
4086          IF (ABS(X2-X1) .LE. EPS*X1) GOTO 40
4087 30    CONTINUE
4088       CALL PUSHERR (0002, 'CALCH4')    ! WARNING ERROR: NO CONVERGENCE
4090 ! *** CONVERGED ; RETURN **********************************************
4092 40    X3 = 0.5*(X1+X2)
4093       Y3 = FUNCH4A (X3)
4095 ! *** CALCULATE HSO4 SPECIATION AND RETURN *******************************
4097 50    CONTINUE
4098       IF (MOLAL(1).GT.TINY .AND. MOLAL(5).GT.TINY) THEN
4099          CALL CALCHS4 (MOLAL(1), MOLAL(5), ZERO, DELTA)
4100          MOLAL(1) = MOLAL(1) - DELTA                      ! H+   EFFECT
4101          MOLAL(5) = MOLAL(5) - DELTA                      ! SO4  EFFECT
4102          MOLAL(6) = DELTA                                 ! HSO4 EFFECT
4103       ENDIF
4105       RETURN
4107 ! *** END OF SUBROUTINE CALCH4 ******************************************
4109       END SUBROUTINE CALCH4
4114 !=======================================================================
4116 ! *** ISORROPIA CODE
4117 ! *** SUBROUTINE FUNCH4A
4118 ! *** CASE H4
4120 !     THE MAIN CHARACTERISTICS OF THIS REGIME ARE:
4121 !     1. SULFATE POOR (SULRAT > 2.0) ; SODIUM RICH (SODRAT >= 2.0)
4122 !     2. THERE IS BOTH A LIQUID & SOLID PHASE
4123 !     3. SOLIDS POSSIBLE : (NH4)2SO4, NH4CL, NA2SO4
4125 ! *** COPYRIGHT 1996-2006, UNIVERSITY OF MIAMI, CARNEGIE MELLON UNIVERSITY,
4126 ! *** GEORGIA INSTITUTE OF TECHNOLOGY
4127 ! *** WRITTEN BY ATHANASIOS NENES
4128 ! *** UPDATED BY CHRISTOS FOUNTOUKIS
4130 !=======================================================================
4132       DOUBLE PRECISION FUNCTION FUNCH4A (X)
4133       USE ISRPIA
4134       USE SOLUT
4135       IMPLICIT DOUBLE PRECISION (A-H, O-Z)
4136 !      INCLUDE 'ISRPIA.INC'
4138 !      COMMON /SOLUT/ CHI1, CHI2, CHI3, CHI4, CHI5, CHI6, CHI7, CHI8,   &
4139 !                     PSI1, PSI2, PSI3, PSI4, PSI5, PSI6, PSI7, PSI8,   &
4140 !                     A1,   A2,   A3,   A4,   A5,   A6,   A7,   A8
4142 ! *** SETUP PARAMETERS ************************************************
4144       PSI6   = X
4145       PSI1   = CHI1
4146       PSI2   = ZERO
4147       PSI3   = ZERO
4148       PSI7   = CHI7
4149       PSI8   = CHI8
4150       FRST   = .TRUE.
4151       CALAIN = .TRUE.
4153 ! *** SOLVE EQUATIONS ; WITH ITERATIONS FOR ACTIVITY COEF. ************
4155       DO 10 I=1,NSWEEP
4157       A1  = XK5 *(WATER/GAMA(2))**3.0
4158       A4  = (XK2/XKW)*R*TEMP*(GAMA(10)/GAMA(5))**2.0
4159       A5  = XK4 *R*TEMP*(WATER/GAMA(10))**2.0
4160       A6  = XK3 *R*TEMP*(WATER/GAMA(11))**2.0
4161       A7  = XK8 *(WATER/GAMA(1))**2.0
4162       A8  = XK9 *(WATER/GAMA(3))**2.0
4163       A9  = XK1*WATER/GAMA(7)*(GAMA(8)/GAMA(7))**2.
4165 !  CALCULATE DISSOCIATION QUANTITIES
4167       PSI5 = CHI5*(PSI6+PSI7) - A6/A5*PSI8*(CHI6-PSI6-PSI3)
4168       PSI5 = PSI5/(A6/A5*(CHI6-PSI6-PSI3) + PSI6 + PSI7)
4169       PSI5 = MAX(PSI5, TINY)
4171       IF (W(3).GT.TINY .AND. WATER.GT.TINY) THEN  ! First try 3rd order soln
4172          BB   =-(CHI4 + PSI6 + PSI5 + 1.d0/A4)
4173          CC   = CHI4*(PSI5+PSI6)
4174          DD   = BB*BB-4.d0*CC
4175          PSI4 =0.5d0*(-BB - SQRT(DD))
4176          PSI4 = MIN(PSI4,CHI4)
4177       ELSE
4178          PSI4 = TINY
4179       ENDIF
4181       IF (CHI1.GT.TINY .AND. WATER.GT.TINY) THEN     ! NA2SO4 DISSOLUTION
4182          AA = PSI7+PSI8
4183          BB = AA*AA
4184          CC =-A1/4.D0
4185          CALL POLY3 (AA, BB, CC, PSI1, ISLV)
4186          IF (ISLV.EQ.0) THEN
4187              PSI1 = MIN (PSI1, CHI1)
4188          ELSE
4189              PSI1 = ZERO
4190          ENDIF
4191       ENDIF
4193 ! *** CALCULATE SPECIATION ********************************************
4195       MOLAL (2) = PSI8 + PSI7 + 2.D0*PSI1                ! NAI
4196       MOLAL (3) = PSI4                                   ! NH4I
4197       MOLAL (4) = PSI6 + PSI7                            ! CLI
4198       MOLAL (5) = PSI2 + PSI1                            ! SO4I
4199       MOLAL (6) = ZERO
4200       MOLAL (7) = PSI5 + PSI8                            ! NO3I
4202       SMIN      = 2.d0*MOLAL(5)+MOLAL(7)+MOLAL(4)-MOLAL(2)-MOLAL(3)
4203       CALL CALCPH (SMIN, HI, OHI)
4204       MOLAL (1) = HI
4206 ! *** CALCULATE GAS / SOLID SPECIES (LIQUID IN MOLAL ALREADY) *********
4208       GNH3      = MAX(CHI4 - PSI4, TINY)
4209       GHNO3     = MAX(CHI5 - PSI5, TINY)
4210       GHCL      = MAX(CHI6 - PSI6, TINY)
4212       CNH42S4   = ZERO
4213       CNH4NO3   = ZERO
4214       CNACL     = MAX(CHI7 - PSI7, ZERO)
4215       CNANO3    = MAX(CHI8 - PSI8, ZERO)
4216       CNA2SO4   = MAX(CHI1 - PSI1, ZERO)
4218 ! *** NH4Cl(s) calculations
4220       A3   = XK6 /(R*TEMP*R*TEMP)
4221       DELT = MIN(GNH3, GHCL)
4222       BB = -(GNH3+GHCL)
4223       CC = GNH3*GHCL-A3
4224       DD = BB*BB - 4.D0*CC
4225       PSI31 = 0.5D0*(-BB + SQRT(DD))
4226       PSI32 = 0.5D0*(-BB - SQRT(DD))
4227       IF (DELT-PSI31.GT.ZERO .AND. PSI31.GT.ZERO) THEN
4228          PSI3 = PSI31
4229       ELSEIF (DELT-PSI32.GT.ZERO .AND. PSI32.GT.ZERO) THEN
4230          PSI3 = PSI32
4231       ELSE
4232          PSI3 = ZERO
4233       ENDIF
4235 ! *** CALCULATE GAS / SOLID SPECIES (LIQUID IN MOLAL ALREADY) *********
4237       GNH3    = MAX(GNH3 - PSI3, TINY)
4238       GHCL    = MAX(GHCL - PSI3, TINY)
4239       CNH4CL  = PSI3
4241       CALL CALCMR                           ! Water content
4243 ! *** CALCULATE ACTIVITIES OR TERMINATE INTERNAL LOOP *****************
4245       IF (FRST.AND.CALAOU .OR. .NOT.FRST.AND.CALAIN) THEN
4246          CALL CALCACT
4247       ELSE
4248          GOTO 20
4249       ENDIF
4250 10    CONTINUE
4252 ! *** CALCULATE FUNCTION VALUE FOR OUTER LOOP ***************************
4254 20    FUNCH4A = MOLAL(3)*MOLAL(4)/GHCL/GNH3/A6/A4 - ONE
4256       RETURN
4258 ! *** END OF FUNCTION FUNCH4A *******************************************
4260       END FUNCTION FUNCH4A    
4262 !=======================================================================
4264 ! *** ISORROPIA CODE
4265 ! *** SUBROUTINE CALCH3
4266 ! *** CASE H3
4268 !     THE MAIN CHARACTERISTICS OF THIS REGIME ARE:
4269 !     1. SULFATE POOR (SULRAT > 2.0) ; SODIUM RICH (SODRAT >= 2.0)
4270 !     2. THERE IS BOTH A LIQUID & SOLID PHASE
4271 !     3. SOLIDS POSSIBLE : NH4CL, NA2SO4
4273 ! *** COPYRIGHT 1996-2006, UNIVERSITY OF MIAMI, CARNEGIE MELLON UNIVERSITY,
4274 ! *** GEORGIA INSTITUTE OF TECHNOLOGY
4275 ! *** WRITTEN BY ATHANASIOS NENES
4276 ! *** UPDATED BY CHRISTOS FOUNTOUKIS
4278 !=======================================================================
4280       SUBROUTINE CALCH3
4281       USE ISRPIA
4282       USE SOLUT
4283       IMPLICIT DOUBLE PRECISION (A-H, O-Z)
4284 !      INCLUDE 'ISRPIA.INC'
4286 !      COMMON /SOLUT/ CHI1, CHI2, CHI3, CHI4, CHI5, CHI6, CHI7, CHI8,   &
4287 !                     PSI1, PSI2, PSI3, PSI4, PSI5, PSI6, PSI7, PSI8,   &
4288 !                     A1,   A2,   A3,   A4,   A5,   A6,   A7,   A8
4290 ! *** REGIME DEPENDS ON THE EXISTANCE OF NITRATES ***********************
4292       IF (W(4).LE.TINY) THEN        ! NO3 NOT EXIST, WATER NOT POSSIBLE
4293          SCASE = 'H3'
4294          CALL CALCH1A
4295          SCASE = 'H3'
4296          RETURN
4297       ENDIF
4299 ! *** SETUP PARAMETERS ************************************************
4301       CALAOU = .TRUE.
4302       CHI1   = W(2)                                ! CNA2SO4
4303       CHI2   = ZERO                                ! CNH42S4
4304       CHI3   = ZERO                                ! CNH4CL
4305       FRNA   = MAX (W(1)-2.D0*CHI1, ZERO)
4306       CHI8   = MIN (FRNA, W(4))                    ! CNANO3
4307       CHI4   = W(3)                                ! NH3(g)
4308       CHI5   = MAX (W(4)-CHI8, ZERO)               ! HNO3(g)
4309       CHI7   = MIN (MAX(FRNA-CHI8, ZERO), W(5))    ! CNACL
4310       CHI6   = MAX (W(5)-CHI7, ZERO)               ! HCL(g)
4312       PSI6LO = TINY
4313       PSI6HI = CHI6-TINY    ! MIN(CHI6-TINY, CHI4)
4315 ! *** INITIAL VALUES FOR BISECTION ************************************
4317       X1 = PSI6LO
4318       Y1 = FUNCH3A (X1)
4319       IF (ABS(Y1).LE.EPS .OR. CHI6.LE.TINY) GOTO 50
4321 ! *** ROOT TRACKING ; FOR THE RANGE OF HI AND LO **********************
4323       DX = (PSI6HI-PSI6LO)/REAL(NDIV)
4324       DO 10 I=1,NDIV
4325          X2 = X1+DX
4326          Y2 = FUNCH3A (X2)
4327          IF (SIGN(1.d0,Y1)*SIGN(1.d0,Y2).LT.ZERO) GOTO 20  ! (Y1*Y2.LT.ZERO)
4328          X1 = X2
4329          Y1 = Y2
4330 10    CONTINUE
4332 ! *** NO SUBDIVISION WITH SOLUTION; IF ABS(Y2)<EPS SOLUTION IS ASSUMED
4334       IF (ABS(Y2) .GT. EPS) Y2 = FUNCH3A (PSI6LO)
4335       GOTO 50
4337 ! *** PERFORM BISECTION ***********************************************
4339 20    DO 30 I=1,MAXIT
4340          X3 = 0.5*(X1+X2)
4341          Y3 = FUNCH3A (X3)
4342          IF (SIGN(1.d0,Y1)*SIGN(1.d0,Y3) .LE. ZERO) THEN  ! (Y1*Y3 .LE. ZERO)
4343             Y2    = Y3
4344             X2    = X3
4345          ELSE
4346             Y1    = Y3
4347             X1    = X3
4348          ENDIF
4349          IF (ABS(X2-X1) .LE. EPS*X1) GOTO 40
4350 30    CONTINUE
4351       CALL PUSHERR (0002, 'CALCH3')    ! WARNING ERROR: NO CONVERGENCE
4353 ! *** CONVERGED ; RETURN **********************************************
4355 40    X3 = 0.5*(X1+X2)
4356       Y3 = FUNCH3A (X3)
4358 ! *** CALCULATE HSO4 SPECIATION AND RETURN *******************************
4360 50    CONTINUE
4361       IF (MOLAL(1).GT.TINY .AND. MOLAL(5).GT.TINY) THEN
4362          CALL CALCHS4 (MOLAL(1), MOLAL(5), ZERO, DELTA)
4363          MOLAL(1) = MOLAL(1) - DELTA                     ! H+   EFFECT
4364          MOLAL(5) = MOLAL(5) - DELTA                     ! SO4  EFFECT
4365          MOLAL(6) = DELTA                                ! HSO4 EFFECT
4366       ENDIF
4368       RETURN
4370 ! *** END OF SUBROUTINE CALCH3 ******************************************
4372       END SUBROUTINE CALCH3
4377 !=======================================================================
4379 ! *** ISORROPIA CODE
4380 ! *** SUBROUTINE FUNCH3A
4381 ! *** CASE H3
4383 !     THE MAIN CHARACTERISTICS OF THIS REGIME ARE:
4384 !     1. SULFATE POOR (SULRAT > 2.0) ; SODIUM RICH (SODRAT >= 2.0)
4385 !     2. THERE IS BOTH A LIQUID & SOLID PHASE
4386 !     3. SOLIDS POSSIBLE : (NH4)2SO4, NH4CL, NA2SO4
4388 ! *** COPYRIGHT 1996-2006, UNIVERSITY OF MIAMI, CARNEGIE MELLON UNIVERSITY,
4389 ! *** GEORGIA INSTITUTE OF TECHNOLOGY
4390 ! *** WRITTEN BY ATHANASIOS NENES
4391 ! *** UPDATED BY CHRISTOS FOUNTOUKIS
4393 !=======================================================================
4395       DOUBLE PRECISION FUNCTION FUNCH3A (X)
4396       USE ISRPIA
4397       USE SOLUT
4398       IMPLICIT DOUBLE PRECISION (A-H, O-Z)
4399 !      INCLUDE 'ISRPIA.INC'
4401 !      COMMON /SOLUT/ CHI1, CHI2, CHI3, CHI4, CHI5, CHI6, CHI7, CHI8,   &
4402 !                     PSI1, PSI2, PSI3, PSI4, PSI5, PSI6, PSI7, PSI8,   &
4403 !                     A1,   A2,   A3,   A4,   A5,   A6,   A7,   A8
4405 ! *** SETUP PARAMETERS ************************************************
4407       PSI6   = X
4408       PSI1   = CHI1
4409       PSI2   = ZERO
4410       PSI3   = ZERO
4411       PSI7   = CHI7
4412       PSI8   = CHI8
4413       FRST   = .TRUE.
4414       CALAIN = .TRUE.
4416 ! *** SOLVE EQUATIONS ; WITH ITERATIONS FOR ACTIVITY COEF. ************
4418       DO 10 I=1,NSWEEP
4420       A1  = XK5 *(WATER/GAMA(2))**3.0
4421       A4  = (XK2/XKW)*R*TEMP*(GAMA(10)/GAMA(5))**2.0
4422       A5  = XK4 *R*TEMP*(WATER/GAMA(10))**2.0
4423       A6  = XK3 *R*TEMP*(WATER/GAMA(11))**2.0
4424       A7  = XK8 *(WATER/GAMA(1))**2.0
4425       A8  = XK9 *(WATER/GAMA(3))**2.0
4426       A9  = XK1*WATER/GAMA(7)*(GAMA(8)/GAMA(7))**2.
4428 !  CALCULATE DISSOCIATION QUANTITIES
4430       PSI5 = CHI5*(PSI6+PSI7) - A6/A5*PSI8*(CHI6-PSI6-PSI3)
4431       PSI5 = PSI5/(A6/A5*(CHI6-PSI6-PSI3) + PSI6 + PSI7)
4432       PSI5 = MAX(PSI5, TINY)
4434       IF (W(3).GT.TINY .AND. WATER.GT.TINY) THEN  ! First try 3rd order soln
4435          BB   =-(CHI4 + PSI6 + PSI5 + 1.d0/A4)
4436          CC   = CHI4*(PSI5+PSI6)
4437          DD   = BB*BB-4.d0*CC
4438          PSI4 =0.5d0*(-BB - SQRT(DD))
4439          PSI4 = MIN(PSI4,CHI4)
4440       ELSE
4441          PSI4 = TINY
4442       ENDIF
4444       IF (CHI7.GT.TINY .AND. WATER.GT.TINY) THEN     ! NACL DISSOLUTION
4445          DIAK = (PSI8-PSI6)**2.D0 + 4.D0*A7
4446          PSI7 = 0.5D0*( -(PSI8+PSI6) + SQRT(DIAK) )
4447          PSI7 = MAX(MIN(PSI7, CHI7), ZERO)
4448       ENDIF
4450       IF (CHI1.GT.TINY .AND. WATER.GT.TINY) THEN     ! NA2SO4 DISSOLUTION
4451          AA = PSI7+PSI8
4452          BB = AA*AA
4453          CC =-A1/4.D0
4454          CALL POLY3 (AA, BB, CC, PSI1, ISLV)
4455          IF (ISLV.EQ.0) THEN
4456              PSI1 = MIN (PSI1, CHI1)
4457          ELSE
4458              PSI1 = ZERO
4459          ENDIF
4460       ENDIF
4462 ! *** CALCULATE SPECIATION ********************************************
4464       MOLAL (2) = PSI8 + PSI7 + 2.D0*PSI1             ! NAI
4465       MOLAL (3) = PSI4                                ! NH4I
4466       MOLAL (4) = PSI6 + PSI7                         ! CLI
4467       MOLAL (5) = PSI2 + PSI1                         ! SO4I
4468       MOLAL (6) = ZERO
4469       MOLAL (7) = PSI5 + PSI8                         ! NO3I
4471       SMIN      = 2.d0*MOLAL(5)+MOLAL(7)+MOLAL(4)-MOLAL(2)-MOLAL(3)
4472       CALL CALCPH (SMIN, HI, OHI)
4473       MOLAL (1) = HI
4475 ! *** CALCULATE GAS / SOLID SPECIES (LIQUID IN MOLAL ALREADY) *********
4477       GNH3      = MAX(CHI4 - PSI4, TINY)
4478       GHNO3     = MAX(CHI5 - PSI5, TINY)
4479       GHCL      = MAX(CHI6 - PSI6, TINY)
4481       CNH42S4   = ZERO
4482       CNH4NO3   = ZERO
4483       CNACL     = MAX(CHI7 - PSI7, ZERO)
4484       CNANO3    = MAX(CHI8 - PSI8, ZERO)
4485       CNA2SO4   = MAX(CHI1 - PSI1, ZERO)
4487 ! *** NH4Cl(s) calculations
4489       A3   = XK6 /(R*TEMP*R*TEMP)
4490       DELT = MIN(GNH3, GHCL)
4491       BB = -(GNH3+GHCL)
4492       CC = GNH3*GHCL-A3
4493       DD = BB*BB - 4.D0*CC
4494       PSI31 = 0.5D0*(-BB + SQRT(DD))
4495       PSI32 = 0.5D0*(-BB - SQRT(DD))
4496       IF (DELT-PSI31.GT.ZERO .AND. PSI31.GT.ZERO) THEN
4497          PSI3 = PSI31
4498       ELSEIF (DELT-PSI32.GT.ZERO .AND. PSI32.GT.ZERO) THEN
4499          PSI3 = PSI32
4500       ELSE
4501          PSI3 = ZERO
4502       ENDIF
4504 ! *** CALCULATE GAS / SOLID SPECIES (LIQUID IN MOLAL ALREADY) *********
4506       GNH3    = MAX(GNH3 - PSI3, TINY)
4507       GHCL    = MAX(GHCL - PSI3, TINY)
4508       CNH4CL  = PSI3
4510       CALL CALCMR                                 ! Water content
4512 ! *** CALCULATE ACTIVITIES OR TERMINATE INTERNAL LOOP *****************
4514       IF (FRST.AND.CALAOU .OR. .NOT.FRST.AND.CALAIN) THEN
4515          CALL CALCACT
4516       ELSE
4517          GOTO 20
4518       ENDIF
4519 10    CONTINUE
4521 ! *** CALCULATE FUNCTION VALUE FOR OUTER LOOP ***************************
4523 20    FUNCH3A = MOLAL(3)*MOLAL(4)/GHCL/GNH3/A6/A4 - ONE
4525       RETURN
4527 ! *** END OF FUNCTION FUNCH3A *******************************************
4529       END FUNCTION FUNCH3A    
4531 !=======================================================================
4533 ! *** ISORROPIA CODE
4534 ! *** SUBROUTINE CALCH2
4535 ! *** CASE H2
4537 !     THE MAIN CHARACTERISTICS OF THIS REGIME ARE:
4538 !     1. SULFATE POOR (SULRAT > 2.0) ; SODIUM RICH (SODRAT >= 2.0)
4539 !     2. SOLID & LIQUID AEROSOL POSSIBLE
4540 !     3. SOLIDS POSSIBLE : NH4Cl, NA2SO4, NANO3, NACL
4542 !     THERE ARE THREE REGIMES IN THIS CASE:
4543 !     1. NH4NO3(s) POSSIBLE. LIQUID & SOLID AEROSOL (SUBROUTINE CALCH2A)
4544 !     2. NH4NO3(s) NOT POSSIBLE, AND RH < MDRH. SOLID AEROSOL ONLY
4545 !     3. NH4NO3(s) NOT POSSIBLE, AND RH >= MDRH. (MDRH REGION)
4547 !     REGIMES 2. AND 3. ARE CONSIDERED TO BE THE SAME AS CASES H1A, H2B
4548 !     RESPECTIVELY (BECAUSE MDRH POINTS COINCIDE).
4550 ! *** COPYRIGHT 1996-2006, UNIVERSITY OF MIAMI, CARNEGIE MELLON UNIVERSITY,
4551 ! *** GEORGIA INSTITUTE OF TECHNOLOGY
4552 ! *** WRITTEN BY ATHANASIOS NENES
4553 ! *** UPDATED BY CHRISTOS FOUNTOUKIS
4555 !=======================================================================
4557       SUBROUTINE CALCH2
4558       USE ISRPIA
4559       IMPLICIT DOUBLE PRECISION (A-H, O-Z)
4560 !      INCLUDE 'ISRPIA.INC'
4561       EXTERNAL CALCH1A, CALCH3
4563 ! *** REGIME DEPENDS ON THE EXISTANCE OF NITRATES ***********************
4565       IF (W(4).GT.TINY) THEN        ! NO3 EXISTS, WATER POSSIBLE
4566          SCASE = 'H2 ; SUBCASE 1'
4567          CALL CALCH2A
4568          SCASE = 'H2 ; SUBCASE 1'
4569       ELSE                          ! NO3 NON EXISTANT, WATER NOT POSSIBLE
4570          SCASE = 'H2 ; SUBCASE 1'
4571          CALL CALCH1A
4572          SCASE = 'H2 ; SUBCASE 1'
4573       ENDIF
4575       IF (WATER.LE.TINY .AND. RH.LT.DRMH2) THEN      ! DRY AEROSOL
4576          SCASE = 'H2 ; SUBCASE 2'
4578       ELSEIF (WATER.LE.TINY .AND. RH.GE.DRMH2) THEN  ! MDRH OF H2
4579          SCASE = 'H2 ; SUBCASE 3'
4580          CALL CALCMDRH (RH, DRMH2, DRNANO3, CALCH1A, CALCH3)
4581          SCASE = 'H2 ; SUBCASE 3'
4582       ENDIF
4584       RETURN
4586 ! *** END OF SUBROUTINE CALCH2 ******************************************
4588       END SUBROUTINE CALCH2
4593 !=======================================================================
4595 ! *** ISORROPIA CODE
4596 ! *** SUBROUTINE CALCH2A
4597 ! *** CASE H2 ; SUBCASE 1
4599 !     THE MAIN CHARACTERISTICS OF THIS REGIME ARE:
4600 !     1. SULFATE POOR (SULRAT > 2.0) ; SODIUM RICH (SODRAT >= 2.0)
4601 !     2. THERE IS BOTH A LIQUID & SOLID PHASE
4602 !     3. SOLIDS POSSIBLE : NH4CL, NA2SO4
4604 ! *** COPYRIGHT 1996-2006, UNIVERSITY OF MIAMI, CARNEGIE MELLON UNIVERSITY,
4605 ! *** GEORGIA INSTITUTE OF TECHNOLOGY
4606 ! *** WRITTEN BY ATHANASIOS NENES
4607 ! *** UPDATED BY CHRISTOS FOUNTOUKIS
4609 !=======================================================================
4611       SUBROUTINE CALCH2A
4612       USE ISRPIA
4613       USE SOLUT
4614       IMPLICIT DOUBLE PRECISION (A-H, O-Z)
4615 !      INCLUDE 'ISRPIA.INC'
4617 !      COMMON /SOLUT/ CHI1, CHI2, CHI3, CHI4, CHI5, CHI6, CHI7, CHI8,   &
4618 !                     PSI1, PSI2, PSI3, PSI4, PSI5, PSI6, PSI7, PSI8,   &
4619 !                     A1,   A2,   A3,   A4,   A5,   A6,   A7,   A8
4621 ! *** SETUP PARAMETERS ************************************************
4623       CALAOU = .TRUE.
4624       CHI1   = W(2)                                ! CNA2SO4
4625       CHI2   = ZERO                                ! CNH42S4
4626       CHI3   = ZERO                                ! CNH4CL
4627       FRNA   = MAX (W(1)-2.D0*CHI1, ZERO)
4628       CHI8   = MIN (FRNA, W(4))                    ! CNANO3
4629       CHI4   = W(3)                                ! NH3(g)
4630       CHI5   = MAX (W(4)-CHI8, ZERO)               ! HNO3(g)
4631       CHI7   = MIN (MAX(FRNA-CHI8, ZERO), W(5))    ! CNACL
4632       CHI6   = MAX (W(5)-CHI7, ZERO)               ! HCL(g)
4634       PSI6LO = TINY
4635       PSI6HI = CHI6-TINY    ! MIN(CHI6-TINY, CHI4)
4637 ! *** INITIAL VALUES FOR BISECTION ************************************
4639       X1 = PSI6LO
4640       Y1 = FUNCH2A (X1)
4641       IF (ABS(Y1).LE.EPS .OR. CHI6.LE.TINY) GOTO 50
4643 ! *** ROOT TRACKING ; FOR THE RANGE OF HI AND LO **********************
4645       DX = (PSI6HI-PSI6LO)/REAL(NDIV)
4646       DO 10 I=1,NDIV
4647          X2 = X1+DX
4648          Y2 = FUNCH2A (X2)
4649          IF (SIGN(1.d0,Y1)*SIGN(1.d0,Y2).LT.ZERO) GOTO 20  ! (Y1*Y2.LT.ZERO)
4650          X1 = X2
4651          Y1 = Y2
4652 10    CONTINUE
4654 ! *** NO SUBDIVISION WITH SOLUTION; IF ABS(Y2)<EPS SOLUTION IS ASSUMED
4656       IF (Y2 .GT. EPS) Y2 = FUNCH2A (PSI6LO)
4657       GOTO 50
4659 ! *** PERFORM BISECTION ***********************************************
4661 20    DO 30 I=1,MAXIT
4662          X3 = 0.5*(X1+X2)
4663          Y3 = FUNCH2A (X3)
4664          IF (SIGN(1.d0,Y1)*SIGN(1.d0,Y3) .LE. ZERO) THEN  ! (Y1*Y3 .LE. ZERO)
4665             Y2    = Y3
4666             X2    = X3
4667          ELSE
4668             Y1    = Y3
4669             X1    = X3
4670          ENDIF
4671          IF (ABS(X2-X1) .LE. EPS*X1) GOTO 40
4672 30    CONTINUE
4673       CALL PUSHERR (0002, 'CALCH2A')    ! WARNING ERROR: NO CONVERGENCE
4675 ! *** CONVERGED ; RETURN **********************************************
4677 40    X3 = 0.5*(X1+X2)
4678       Y3 = FUNCH2A (X3)
4680 ! *** CALCULATE HSO4 SPECIATION AND RETURN *******************************
4682 50    CONTINUE
4683       IF (MOLAL(1).GT.TINY .AND. MOLAL(5).GT.TINY) THEN
4684          CALL CALCHS4 (MOLAL(1), MOLAL(5), ZERO, DELTA)
4685          MOLAL(1) = MOLAL(1) - DELTA                    ! H+   EFFECT
4686          MOLAL(5) = MOLAL(5) - DELTA                    ! SO4  EFFECT
4687          MOLAL(6) = DELTA                               ! HSO4 EFFECT
4688       ENDIF
4690       RETURN
4692 ! *** END OF SUBROUTINE CALCH2A ******************************************
4694       END SUBROUTINE CALCH2A
4699 !=======================================================================
4701 ! *** ISORROPIA CODE
4702 ! *** SUBROUTINE FUNCH2A
4703 ! *** CASE H2 ; SUBCASE 1
4705 !     THE MAIN CHARACTERISTICS OF THIS REGIME ARE:
4706 !     1. SULFATE POOR (SULRAT > 2.0) ; SODIUM RICH (SODRAT >= 2.0)
4707 !     2. THERE IS BOTH A LIQUID & SOLID PHASE
4708 !     3. SOLIDS POSSIBLE : (NH4)2SO4, NH4CL, NA2SO4
4710 ! *** COPYRIGHT 1996-2006, UNIVERSITY OF MIAMI, CARNEGIE MELLON UNIVERSITY,
4711 ! *** GEORGIA INSTITUTE OF TECHNOLOGY
4712 ! *** WRITTEN BY ATHANASIOS NENES
4713 ! *** UPDATED BY CHRISTOS FOUNTOUKIS
4715 !=======================================================================
4717       DOUBLE PRECISION FUNCTION FUNCH2A (X)
4718       USE ISRPIA
4719       USE SOLUT
4720       IMPLICIT DOUBLE PRECISION (A-H, O-Z)
4721 !      INCLUDE 'ISRPIA.INC'
4723 !      COMMON /SOLUT/ CHI1, CHI2, CHI3, CHI4, CHI5, CHI6, CHI7, CHI8,   &
4724 !                     PSI1, PSI2, PSI3, PSI4, PSI5, PSI6, PSI7, PSI8,   &
4725 !                     A1,   A2,   A3,   A4,   A5,   A6,   A7,   A8
4727 ! *** SETUP PARAMETERS ************************************************
4729       PSI6   = X
4730       PSI1   = CHI1
4731       PSI2   = ZERO
4732       PSI3   = ZERO
4733       PSI7   = CHI7
4734       PSI8   = CHI8
4735       FRST   = .TRUE.
4736       CALAIN = .TRUE.
4738 ! *** SOLVE EQUATIONS ; WITH ITERATIONS FOR ACTIVITY COEF. ************
4740       DO 10 I=1,NSWEEP
4742       A1  = XK5 *(WATER/GAMA(2))**3.0
4743       A4  = (XK2/XKW)*R*TEMP*(GAMA(10)/GAMA(5))**2.0
4744       A5  = XK4 *R*TEMP*(WATER/GAMA(10))**2.0
4745       A6  = XK3 *R*TEMP*(WATER/GAMA(11))**2.0
4746       A7  = XK8 *(WATER/GAMA(1))**2.0
4747       A8  = XK9 *(WATER/GAMA(3))**2.0
4748       A64 = (XK3*XK2/XKW)*(GAMA(10)/GAMA(5)/GAMA(11))**2.0
4749       A64 = A64*(R*TEMP*WATER)**2.0
4750       A9  = XK1*WATER/GAMA(7)*(GAMA(8)/GAMA(7))**2.
4752 !  CALCULATE DISSOCIATION QUANTITIES
4754       PSI5 = CHI5*(PSI6+PSI7) - A6/A5*PSI8*(CHI6-PSI6-PSI3)
4755       PSI5 = PSI5/(A6/A5*(CHI6-PSI6-PSI3) + PSI6 + PSI7)
4756       PSI5 = MAX(PSI5, TINY)
4758       IF (W(3).GT.TINY .AND. WATER.GT.TINY) THEN  ! First try 3rd order soln
4759          BB   =-(CHI4 + PSI6 + PSI5 + 1.d0/A4)
4760          CC   = CHI4*(PSI5+PSI6)
4761          DD   = BB*BB-4.d0*CC
4762          PSI4 =0.5d0*(-BB - SQRT(DD))
4763          PSI4 = MIN(PSI4,CHI4)
4764       ELSE
4765          PSI4 = TINY
4766       ENDIF
4768       IF (CHI7.GT.TINY .AND. WATER.GT.TINY) THEN     ! NACL DISSOLUTION
4769          DIAK = (PSI8-PSI6)**2.D0 + 4.D0*A7
4770          PSI7 = 0.5D0*( -(PSI8+PSI6) + SQRT(DIAK) )
4771          PSI7 = MAX(MIN(PSI7, CHI7), ZERO)
4772       ENDIF
4774       IF (CHI8.GT.TINY .AND. WATER.GT.TINY) THEN     ! NANO3 DISSOLUTION
4775          DIAK = (PSI7-PSI5)**2.D0 + 4.D0*A8
4776          PSI8 = 0.5D0*( -(PSI7+PSI5) + SQRT(DIAK) )
4777          PSI8 = MAX(MIN(PSI8, CHI8), ZERO)
4778       ENDIF
4780       IF (CHI1.GT.TINY .AND. WATER.GT.TINY) THEN     ! NA2SO4 DISSOLUTION
4781          AA = PSI7+PSI8
4782          BB = AA*AA
4783          CC =-A1/4.D0
4784          CALL POLY3 (AA, BB, CC, PSI1, ISLV)
4785          IF (ISLV.EQ.0) THEN
4786              PSI1 = MIN (PSI1, CHI1)
4787          ELSE
4788              PSI1 = ZERO
4789          ENDIF
4790       ENDIF
4792 ! *** CALCULATE SPECIATION ********************************************
4794       MOLAL (2) = PSI8 + PSI7 + 2.D0*PSI1                 ! NAI
4795       MOLAL (3) = PSI4                                    ! NH4I
4796       MOLAL (4) = PSI6 + PSI7                             ! CLI
4797       MOLAL (5) = PSI2 + PSI1                             ! SO4I
4798       MOLAL (6) = ZERO                                    ! HSO4I
4799       MOLAL (7) = PSI5 + PSI8                             ! NO3I
4801       SMIN      = 2.d0*MOLAL(5)+MOLAL(7)+MOLAL(4)-MOLAL(2)-MOLAL(3)
4802       CALL CALCPH (SMIN, HI, OHI)
4803       MOLAL (1) = HI
4805       GNH3      = MAX(CHI4 - PSI4, TINY)
4806       GHNO3     = MAX(CHI5 - PSI5, TINY)
4807       GHCL      = MAX(CHI6 - PSI6, TINY)
4809       CNH42S4   = ZERO
4810       CNH4NO3   = ZERO
4811       CNACL     = MAX(CHI7 - PSI7, ZERO)
4812       CNANO3    = MAX(CHI8 - PSI8, ZERO)
4813       CNA2SO4   = MAX(CHI1 - PSI1, ZERO)
4815 ! *** NH4Cl(s) calculations
4817       A3   = XK6 /(R*TEMP*R*TEMP)
4818       DELT = MIN(GNH3, GHCL)
4819       BB = -(GNH3+GHCL)
4820       CC = GNH3*GHCL-A3
4821       DD = BB*BB - 4.D0*CC
4822       PSI31 = 0.5D0*(-BB + SQRT(DD))
4823       PSI32 = 0.5D0*(-BB - SQRT(DD))
4824       IF (DELT-PSI31.GT.ZERO .AND. PSI31.GT.ZERO) THEN
4825          PSI3 = PSI31
4826       ELSEIF (DELT-PSI32.GT.ZERO .AND. PSI32.GT.ZERO) THEN
4827          PSI3 = PSI32
4828       ELSE
4829          PSI3 = ZERO
4830       ENDIF
4832 ! *** CALCULATE GAS / SOLID SPECIES (LIQUID IN MOLAL ALREADY) *********
4834       GNH3    = MAX(GNH3 - PSI3, TINY)
4835       GHCL    = MAX(GHCL - PSI3, TINY)
4836       CNH4CL  = PSI3
4838       CALL CALCMR                        ! Water content
4840 ! *** CALCULATE ACTIVITIES OR TERMINATE INTERNAL LOOP *****************
4842       IF (FRST.AND.CALAOU .OR. .NOT.FRST.AND.CALAIN) THEN
4843          CALL CALCACT
4844       ELSE
4845          GOTO 20
4846       ENDIF
4847 10    CONTINUE
4849 ! *** CALCULATE FUNCTION VALUE FOR OUTER LOOP ***************************
4851 20    FUNCH2A = MOLAL(3)*MOLAL(4)/GHCL/GNH3/A64 - ONE
4853       RETURN
4855 ! *** END OF FUNCTION FUNCH2A *******************************************
4857       END FUNCTION FUNCH2A    
4860 !=======================================================================
4862 ! *** ISORROPIA CODE
4863 ! *** SUBROUTINE CALCH1
4864 ! *** CASE H1
4866 !     THE MAIN CHARACTERISTICS OF THIS REGIME ARE:
4867 !     1. SULFATE POOR (SULRAT > 2.0) ; SODIUM RICH (SODRAT >= 2.0)
4868 !     2. SOLID AEROSOL ONLY
4869 !     3. SOLIDS POSSIBLE : NH4NO3, NH4CL, NA2SO4
4871 !     THERE ARE TWO POSSIBLE REGIMES HERE, DEPENDING ON RELATIVE HUMIDITY:
4872 !     1. WHEN RH >= MDRH ; LIQUID PHASE POSSIBLE (MDRH REGION)
4873 !     2. WHEN RH < MDRH  ; ONLY SOLID PHASE POSSIBLE (SUBROUTINE CALCH1A)
4875 ! *** COPYRIGHT 1996-2006, UNIVERSITY OF MIAMI, CARNEGIE MELLON UNIVERSITY,
4876 ! *** GEORGIA INSTITUTE OF TECHNOLOGY
4877 ! *** WRITTEN BY ATHANASIOS NENES
4878 ! *** UPDATED BY CHRISTOS FOUNTOUKIS
4880 !=======================================================================
4882       SUBROUTINE CALCH1
4883       USE ISRPIA
4884       IMPLICIT DOUBLE PRECISION (A-H, O-Z)
4885 !      INCLUDE 'ISRPIA.INC'
4886       EXTERNAL CALCH1A, CALCH2A
4888 ! *** REGIME DEPENDS UPON THE AMBIENT RELATIVE HUMIDITY *****************
4890       IF (RH.LT.DRMH1) THEN
4891          SCASE = 'H1 ; SUBCASE 1'
4892          CALL CALCH1A              ! SOLID PHASE ONLY POSSIBLE
4893          SCASE = 'H1 ; SUBCASE 1'
4894       ELSE
4895          SCASE = 'H1 ; SUBCASE 2'  ! LIQUID & SOLID PHASE POSSIBLE
4896          CALL CALCMDRH (RH, DRMH1, DRNH4NO3, CALCH1A, CALCH2A)
4897          SCASE = 'H1 ; SUBCASE 2'
4898       ENDIF
4900       RETURN
4902 ! *** END OF SUBROUTINE CALCH1 ******************************************
4904       END SUBROUTINE CALCH1
4907 !=======================================================================
4909 ! *** ISORROPIA CODE
4910 ! *** SUBROUTINE CALCH1A
4911 ! *** CASE H1 ; SUBCASE 1
4913 !     THE MAIN CHARACTERISTICS OF THIS REGIME ARE:
4914 !     1. SULFATE POOR (SULRAT > 2.0) ; SODIUM RICH (SODRAT >= 2.0)
4915 !     2. SOLID AEROSOL ONLY
4916 !     3. SOLIDS POSSIBLE : NH4NO3, NH4CL, NANO3, NA2SO4
4918 ! *** COPYRIGHT 1996-2006, UNIVERSITY OF MIAMI, CARNEGIE MELLON UNIVERSITY,
4919 ! *** GEORGIA INSTITUTE OF TECHNOLOGY
4920 ! *** WRITTEN BY ATHANASIOS NENES
4921 ! *** UPDATED BY CHRISTOS FOUNTOUKIS
4923 !=======================================================================
4925       SUBROUTINE CALCH1A
4926       USE ISRPIA
4927       IMPLICIT DOUBLE PRECISION (A-H, O-Z)
4928 !      INCLUDE 'ISRPIA.INC'
4929       DOUBLE PRECISION LAMDA, LAMDA1, LAMDA2, KAPA, KAPA1, KAPA2, NAFR,   &
4930                        NO3FR
4932 ! *** CALCULATE NON VOLATILE SOLIDS ***********************************
4934       CNA2SO4 = W(2)
4935       CNH42S4 = ZERO
4936       NAFR    = MAX (W(1)-2*CNA2SO4, ZERO)
4937       CNANO3  = MIN (NAFR, W(4))
4938       NO3FR   = MAX (W(4)-CNANO3, ZERO)
4939       CNACL   = MIN (MAX(NAFR-CNANO3, ZERO), W(5))
4940       CLFR    = MAX (W(5)-CNACL, ZERO)
4942 ! *** CALCULATE VOLATILE SPECIES **************************************
4944       ALF     = W(3)                     ! FREE NH3
4945       BET     = CLFR                     ! FREE CL
4946       GAM     = NO3FR                    ! FREE NO3
4948       RTSQ    = R*TEMP*R*TEMP
4949       A1      = XK6/RTSQ
4950       A2      = XK10/RTSQ
4952       THETA1  = GAM - BET*(A2/A1)
4953       THETA2  = A2/A1
4955 ! QUADRATIC EQUATION SOLUTION
4957       BB      = (THETA1-ALF-BET*(ONE+THETA2))/(ONE+THETA2)
4958       CC      = (ALF*BET-A1-BET*THETA1)/(ONE+THETA2)
4959       DD      = BB*BB - 4.0D0*CC
4960       IF (DD.LT.ZERO) GOTO 100   ! Solve each reaction seperately
4962 ! TWO ROOTS FOR KAPA, CHECK AND SEE IF ANY VALID
4964       SQDD    = SQRT(DD)
4965       KAPA1   = 0.5D0*(-BB+SQDD)
4966       KAPA2   = 0.5D0*(-BB-SQDD)
4967       LAMDA1  = THETA1 + THETA2*KAPA1
4968       LAMDA2  = THETA1 + THETA2*KAPA2
4970       IF (KAPA1.GE.ZERO .AND. LAMDA1.GE.ZERO) THEN
4971          IF (ALF-KAPA1-LAMDA1.GE.ZERO .AND.   &
4972              BET-KAPA1.GE.ZERO .AND. GAM-LAMDA1.GE.ZERO) THEN
4973              KAPA = KAPA1
4974              LAMDA= LAMDA1
4975              GOTO 200
4976          ENDIF
4977       ENDIF
4979       IF (KAPA2.GE.ZERO .AND. LAMDA2.GE.ZERO) THEN
4980          IF (ALF-KAPA2-LAMDA2.GE.ZERO .AND.   &
4981              BET-KAPA2.GE.ZERO .AND. GAM-LAMDA2.GE.ZERO) THEN
4982              KAPA = KAPA2
4983              LAMDA= LAMDA2
4984              GOTO 200
4985          ENDIF
4986       ENDIF
4988 ! SEPERATE SOLUTION OF NH4CL & NH4NO3 EQUILIBRIA
4990 100   KAPA  = ZERO
4991       LAMDA = ZERO
4992       DD1   = (ALF+BET)*(ALF+BET) - 4.0D0*(ALF*BET-A1)
4993       DD2   = (ALF+GAM)*(ALF+GAM) - 4.0D0*(ALF*GAM-A2)
4995 ! NH4CL EQUILIBRIUM
4997       IF (DD1.GE.ZERO) THEN
4998          SQDD1 = SQRT(DD1)
4999          KAPA1 = 0.5D0*(ALF+BET + SQDD1)
5000          KAPA2 = 0.5D0*(ALF+BET - SQDD1)
5002          IF (KAPA1.GE.ZERO .AND. KAPA1.LE.MIN(ALF,BET)) THEN
5003             KAPA = KAPA1
5004          ELSE IF (KAPA2.GE.ZERO .AND. KAPA2.LE.MIN(ALF,BET)) THEN
5005             KAPA = KAPA2
5006          ELSE
5007             KAPA = ZERO
5008          ENDIF
5009       ENDIF
5011 ! NH4NO3 EQUILIBRIUM
5013       IF (DD2.GE.ZERO) THEN
5014          SQDD2 = SQRT(DD2)
5015          LAMDA1= 0.5D0*(ALF+GAM + SQDD2)
5016          LAMDA2= 0.5D0*(ALF+GAM - SQDD2)
5018          IF (LAMDA1.GE.ZERO .AND. LAMDA1.LE.MIN(ALF,GAM)) THEN
5019             LAMDA = LAMDA1
5020          ELSE IF (LAMDA2.GE.ZERO .AND. LAMDA2.LE.MIN(ALF,GAM)) THEN
5021             LAMDA = LAMDA2
5022          ELSE
5023             LAMDA = ZERO
5024          ENDIF
5025       ENDIF
5027 ! IF BOTH KAPA, LAMDA ARE > 0, THEN APPLY EXISTANCE CRITERION
5029       IF (KAPA.GT.ZERO .AND. LAMDA.GT.ZERO) THEN
5030          IF (BET .LT. LAMDA/THETA1) THEN
5031             KAPA = ZERO
5032          ELSE
5033             LAMDA= ZERO
5034          ENDIF
5035       ENDIF
5037 ! *** CALCULATE COMPOSITION OF VOLATILE SPECIES ***********************
5039 200   CONTINUE
5040       CNH4NO3 = LAMDA
5041       CNH4CL  = KAPA
5043       GNH3    = ALF - KAPA - LAMDA
5044       GHNO3   = GAM - LAMDA
5045       GHCL    = BET - KAPA
5047       RETURN
5049 ! *** END OF SUBROUTINE CALCH1A *****************************************
5051       END SUBROUTINE CALCH1A
5052 !=======================================================================
5054 ! *** ISORROPIA CODE
5055 ! *** SUBROUTINE CALCI6
5056 ! *** CASE I6
5058 !     THE MAIN CHARACTERISTICS OF THIS REGIME ARE:
5059 !     1. SULFATE RICH, NO FREE ACID (1.0 <= SULRAT < 2.0)
5060 !     2. SOLID & LIQUID AEROSOL POSSIBLE
5061 !     3. SOLIDS POSSIBLE : (NH4)2SO4, NA2SO4
5063 ! *** COPYRIGHT 1996-2006, UNIVERSITY OF MIAMI, CARNEGIE MELLON UNIVERSITY,
5064 ! *** GEORGIA INSTITUTE OF TECHNOLOGY
5065 ! *** WRITTEN BY ATHANASIOS NENES
5066 ! *** UPDATED BY CHRISTOS FOUNTOUKIS
5068 !=======================================================================
5070       SUBROUTINE CALCI6
5071       USE ISRPIA
5072       USE SOLUT
5073       IMPLICIT DOUBLE PRECISION (A-H, O-Z)
5074 !      INCLUDE 'ISRPIA.INC'
5075 !      COMMON /SOLUT/ CHI1, CHI2, CHI3, CHI4, CHI5, CHI6, CHI7, CHI8,   &
5076 !                     PSI1, PSI2, PSI3, PSI4, PSI5, PSI6, PSI7, PSI8,   &
5077 !                     A1,   A2,   A3,   A4,   A5,   A6,   A7,   A8
5079 ! *** FIND DRY COMPOSITION **********************************************
5081       CALL CALCI1A
5083 ! *** SETUP PARAMETERS ************************************************
5085       CHI1 = CNH4HS4               ! Save from CALCI1 run
5086       CHI2 = CLC
5087       CHI3 = CNAHSO4
5088       CHI4 = CNA2SO4
5089       CHI5 = CNH42S4
5091       PSI1 = CNH4HS4               ! ASSIGN INITIAL PSI's
5092       PSI2 = CLC
5093       PSI3 = CNAHSO4
5094       PSI4 = CNA2SO4
5095       PSI5 = CNH42S4
5097       CALAOU = .TRUE.              ! Outer loop activity calculation flag
5098       FRST   = .TRUE.
5099       CALAIN = .TRUE.
5101 ! *** SOLVE EQUATIONS ; WITH ITERATIONS FOR ACTIVITY COEF. ************
5103       DO 10 I=1,NSWEEP
5105       A6 = XK1 *WATER/GAMA(7)*(GAMA(8)/GAMA(7))**2.
5107 !  CALCULATE DISSOCIATION QUANTITIES
5109       BB   = PSI2 + PSI4 + PSI5 + A6                    ! PSI6
5110       CC   =-A6*(PSI2 + PSI3 + PSI1)
5111       DD   = BB*BB - 4.D0*CC
5112       PSI6 = 0.5D0*(-BB + SQRT(DD))
5114 ! *** CALCULATE SPECIATION ********************************************
5116       MOLAL (1) = PSI6                                    ! HI
5117       MOLAL (2) = 2.D0*PSI4 + PSI3                        ! NAI
5118       MOLAL (3) = 3.D0*PSI2 + 2.D0*PSI5 + PSI1            ! NH4I
5119       MOLAL (5) = PSI2 + PSI4 + PSI5 + PSI6               ! SO4I
5120       MOLAL (6) = PSI2 + PSI3 + PSI1 - PSI6               ! HSO4I
5121       CLC       = ZERO
5122       CNAHSO4   = ZERO
5123       CNA2SO4   = CHI4 - PSI4
5124       CNH42S4   = ZERO
5125       CNH4HS4   = ZERO
5126       CALL CALCMR                                         ! Water content
5128 ! *** CALCULATE ACTIVITIES OR TERMINATE INTERNAL LOOP *****************
5130       IF (FRST.AND.CALAOU .OR. .NOT.FRST.AND.CALAIN) THEN
5131          CALL CALCACT
5132       ELSE
5133          GOTO 20
5134       ENDIF
5135 10    CONTINUE
5137 20    RETURN
5139 ! *** END OF SUBROUTINE CALCI6 *****************************************
5141       END SUBROUTINE CALCI6
5143 !=======================================================================
5145 ! *** ISORROPIA CODE
5146 ! *** SUBROUTINE CALCI5
5147 ! *** CASE I5
5149 !     THE MAIN CHARACTERISTICS OF THIS REGIME ARE:
5150 !     1. SULFATE RICH, NO FREE ACID (1.0 <= SULRAT < 2.0)
5151 !     2. SOLID & LIQUID AEROSOL POSSIBLE
5152 !     3. SOLIDS POSSIBLE : (NH4)2SO4, NA2SO4
5154 ! *** COPYRIGHT 1996-2006, UNIVERSITY OF MIAMI, CARNEGIE MELLON UNIVERSITY,
5155 ! *** GEORGIA INSTITUTE OF TECHNOLOGY
5156 ! *** WRITTEN BY ATHANASIOS NENES
5157 ! *** UPDATED BY CHRISTOS FOUNTOUKIS
5159 !=======================================================================
5161       SUBROUTINE CALCI5
5162       USE ISRPIA
5163       USE SOLUT
5164       IMPLICIT DOUBLE PRECISION (A-H, O-Z)
5165 !      INCLUDE 'ISRPIA.INC'
5166 !      COMMON /SOLUT/ CHI1, CHI2, CHI3, CHI4, CHI5, CHI6, CHI7, CHI8,   &
5167 !                     PSI1, PSI2, PSI3, PSI4, PSI5, PSI6, PSI7, PSI8,   &
5168 !                     A1,   A2,   A3,   A4,   A5,   A6,   A7,   A8
5170 ! *** FIND DRY COMPOSITION **********************************************
5172       CALL CALCI1A
5174 ! *** SETUP PARAMETERS ************************************************
5176       CHI1 = CNH4HS4               ! Save from CALCI1 run
5177       CHI2 = CLC
5178       CHI3 = CNAHSO4
5179       CHI4 = CNA2SO4
5180       CHI5 = CNH42S4
5182       PSI1 = CNH4HS4               ! ASSIGN INITIAL PSI's
5183       PSI2 = CLC
5184       PSI3 = CNAHSO4
5185       PSI4 = ZERO
5186       PSI5 = CNH42S4
5188       CALAOU =.TRUE.               ! Outer loop activity calculation flag
5189       PSI4LO = ZERO                ! Low  limit
5190       PSI4HI = CHI4                ! High limit
5192 ! *** IF NA2SO4(S) =0, CALL FUNCI5B FOR Y4=0 ***************************
5194       IF (CHI4.LE.TINY) THEN
5195          Y1 = FUNCI5A (ZERO)
5196          GOTO 50
5197       ENDIF
5199 ! *** INITIAL VALUES FOR BISECTION ************************************
5201       X1 = PSI4HI
5202       Y1 = FUNCI5A (X1)
5203       YHI= Y1                      ! Save Y-value at HI position
5205 ! *** YHI < 0.0 THE SOLUTION IS ALWAYS UNDERSATURATED WITH NA2SO4 **
5207       IF (ABS(Y1).LE.EPS .OR. YHI.LT.ZERO) GOTO 50
5209 ! *** ROOT TRACKING ; FOR THE RANGE OF HI AND LO **********************
5211       DX = (PSI4HI-PSI4LO)/REAL(NDIV)
5212       DO 10 I=1,NDIV
5213          X2 = X1-DX
5214          Y2 = FUNCI5A (X2)
5215          IF (SIGN(1.d0,Y1)*SIGN(1.d0,Y2).LT.ZERO) GOTO 20  ! (Y1*Y2.LT.ZERO)
5216          X1 = X2
5217          Y1 = Y2
5218 10    CONTINUE
5220 ! *** { YLO, YHI } > 0.0 THE SOLUTION IS ALWAYS SUPERSATURATED WITH NH4CL
5222       YLO= Y1                      ! Save Y-value at Hi position
5223       IF (YLO.GT.ZERO .AND. YHI.GT.ZERO) THEN
5224          Y3 = FUNCI5A (ZERO)
5225          GOTO 50
5226       ELSE IF (ABS(Y2) .LT. EPS) THEN   ! X2 IS A SOLUTION
5227          GOTO 50
5228       ELSE
5229          CALL PUSHERR (0001, 'CALCI5')    ! WARNING ERROR: NO SOLUTION
5230          GOTO 50
5231       ENDIF
5233 ! *** PERFORM BISECTION ***********************************************
5235 20    DO 30 I=1,MAXIT
5236          X3 = 0.5*(X1+X2)
5237          Y3 = FUNCI5A (X3)
5238          IF (SIGN(1.d0,Y1)*SIGN(1.d0,Y3) .LE. ZERO) THEN  ! (Y1*Y3 .LE. ZERO)
5239             Y2    = Y3
5240             X2    = X3
5241          ELSE
5242             Y1    = Y3
5243             X1    = X3
5244          ENDIF
5245          IF (ABS(X2-X1) .LE. EPS*X1) GOTO 40
5246 30    CONTINUE
5247       CALL PUSHERR (0002, 'CALCI5')    ! WARNING ERROR: NO CONVERGENCE
5249 ! *** CONVERGED ; RETURN **********************************************
5251 40    X3 = 0.5*(X1+X2)
5252       Y3 = FUNCI5A (X3)
5254 50    RETURN
5256 ! *** END OF SUBROUTINE CALCI5 *****************************************
5258       END SUBROUTINE CALCI5
5263 !=======================================================================
5265 ! *** ISORROPIA CODE
5266 ! *** SUBROUTINE FUNCI5A
5267 ! *** CASE I5
5269 !     THE MAIN CHARACTERISTICS OF THIS REGIME ARE:
5270 !     1. SULFATE RICH, NO FREE ACID (1.0 <= SULRAT < 2.0)
5271 !     2. SOLID & LIQUID AEROSOL POSSIBLE
5272 !     3. SOLIDS POSSIBLE : NA2SO4
5274 ! *** COPYRIGHT 1996-2006, UNIVERSITY OF MIAMI, CARNEGIE MELLON UNIVERSITY,
5275 ! *** GEORGIA INSTITUTE OF TECHNOLOGY
5276 ! *** WRITTEN BY ATHANASIOS NENES
5277 ! *** UPDATED BY CHRISTOS FOUNTOUKIS
5279 !=======================================================================
5281       DOUBLE PRECISION FUNCTION FUNCI5A (P4)
5282       USE ISRPIA
5283       USE SOLUT
5284       IMPLICIT DOUBLE PRECISION (A-H, O-Z)
5285 !      INCLUDE 'ISRPIA.INC'
5286 !      COMMON /SOLUT/ CHI1, CHI2, CHI3, CHI4, CHI5, CHI6, CHI7, CHI8,   &
5287 !                     PSI1, PSI2, PSI3, PSI4, PSI5, PSI6, PSI7, PSI8,   &
5288 !                     A1,   A2,   A3,   A4,   A5,   A6,   A7,   A8
5290 ! *** SETUP PARAMETERS ************************************************
5292       PSI4   = P4     ! PSI3 already assigned in FUNCI5A
5293       FRST   = .TRUE.
5294       CALAIN = .TRUE.
5296 ! *** SOLVE EQUATIONS ; WITH ITERATIONS FOR ACTIVITY COEF. ************
5298       DO 10 I=1,NSWEEP
5300       A4 = XK5 *(WATER/GAMA(2))**3.0
5301       A5 = XK7 *(WATER/GAMA(4))**3.0
5302       A6 = XK1 *WATER/GAMA(7)*(GAMA(8)/GAMA(7))**2.
5304 !  CALCULATE DISSOCIATION QUANTITIES
5306       BB   = PSI2 + PSI4 + PSI5 + A6                    ! PSI6
5307       CC   =-A6*(PSI2 + PSI3 + PSI1)
5308       DD   = BB*BB - 4.D0*CC
5309       PSI6 = 0.5D0*(-BB + SQRT(DD))
5311 ! *** CALCULATE SPECIATION ********************************************
5313       MOLAL (1) = PSI6                            ! HI
5314       MOLAL (2) = 2.D0*PSI4 + PSI3                ! NAI
5315       MOLAL (3) = 3.D0*PSI2 + 2.D0*PSI5 + PSI1    ! NH4I
5316       MOLAL (5) = PSI2 + PSI4 + PSI5 + PSI6       ! SO4I
5317       MOLAL (6) = PSI2 + PSI3 + PSI1 - PSI6       ! HSO4I
5318       CLC       = ZERO
5319       CNAHSO4   = ZERO
5320       CNA2SO4   = CHI4 - PSI4
5321       CNH42S4   = ZERO
5322       CNH4HS4   = ZERO
5323       CALL CALCMR                                 ! Water content
5325 ! *** CALCULATE ACTIVITIES OR TERMINATE INTERNAL LOOP *****************
5327       IF (FRST.AND.CALAOU .OR. .NOT.FRST.AND.CALAIN) THEN
5328          CALL CALCACT
5329       ELSE
5330          GOTO 20
5331       ENDIF
5332 10    CONTINUE
5334 ! *** CALCULATE OBJECTIVE FUNCTION ************************************
5336 20    A4     = XK5 *(WATER/GAMA(2))**3.0
5337       FUNCI5A= MOLAL(5)*MOLAL(2)*MOLAL(2)/A4 - ONE
5338       RETURN
5340 ! *** END OF FUNCTION FUNCI5A ********************************************
5342       END FUNCTION FUNCI5A     
5343 !=======================================================================
5345 ! *** ISORROPIA CODE
5346 ! *** SUBROUTINE CALCI4
5347 ! *** CASE I4
5349 !     THE MAIN CHARACTERISTICS OF THIS REGIME ARE:
5350 !     1. SULFATE RICH, NO FREE ACID (1.0 <= SULRAT < 2.0)
5351 !     2. SOLID & LIQUID AEROSOL POSSIBLE
5352 !     3. SOLIDS POSSIBLE : (NH4)2SO4, NA2SO4
5354 ! *** COPYRIGHT 1996-2006, UNIVERSITY OF MIAMI, CARNEGIE MELLON UNIVERSITY,
5355 ! *** GEORGIA INSTITUTE OF TECHNOLOGY
5356 ! *** WRITTEN BY ATHANASIOS NENES
5357 ! *** UPDATED BY CHRISTOS FOUNTOUKIS
5359 !=======================================================================
5361       SUBROUTINE CALCI4
5362       USE ISRPIA
5363       USE SOLUT
5364       IMPLICIT DOUBLE PRECISION (A-H, O-Z)
5365 !      INCLUDE 'ISRPIA.INC'
5366 !      COMMON /SOLUT/ CHI1, CHI2, CHI3, CHI4, CHI5, CHI6, CHI7, CHI8,   &
5367 !                     PSI1, PSI2, PSI3, PSI4, PSI5, PSI6, PSI7, PSI8,   &
5368 !                     A1,   A2,   A3,   A4,   A5,   A6,   A7,   A8
5370 ! *** FIND DRY COMPOSITION **********************************************
5372       CALL CALCI1A
5374 ! *** SETUP PARAMETERS ************************************************
5376       CHI1 = CNH4HS4               ! Save from CALCI1 run
5377       CHI2 = CLC
5378       CHI3 = CNAHSO4
5379       CHI4 = CNA2SO4
5380       CHI5 = CNH42S4
5382       PSI1 = CNH4HS4               ! ASSIGN INITIAL PSI's
5383       PSI2 = CLC
5384       PSI3 = CNAHSO4
5385       PSI4 = ZERO
5386       PSI5 = ZERO
5388       CALAOU = .TRUE.              ! Outer loop activity calculation flag
5389       PSI4LO = ZERO                ! Low  limit
5390       PSI4HI = CHI4                ! High limit
5392 ! *** IF NA2SO4(S) =0, CALL FUNCI4B FOR Y4=0 ***************************
5394       IF (CHI4.LE.TINY) THEN
5395          Y1 = FUNCI4A (ZERO)
5396          GOTO 50
5397       ENDIF
5399 ! *** INITIAL VALUES FOR BISECTION ************************************
5401       X1 = PSI4HI
5402       Y1 = FUNCI4A (X1)
5403       YHI= Y1                      ! Save Y-value at HI position
5405 ! *** YHI < 0.0 THE SOLUTION IS ALWAYS UNDERSATURATED WITH NA2SO4 **
5407       IF (ABS(Y1).LE.EPS .OR. YHI.LT.ZERO) GOTO 50
5409 ! *** ROOT TRACKING ; FOR THE RANGE OF HI AND LO **********************
5411       DX = (PSI4HI-PSI4LO)/REAL(NDIV)
5412       DO 10 I=1,NDIV
5413          X2 = X1-DX
5414          Y2 = FUNCI4A (X2)
5415          IF (SIGN(1.d0,Y1)*SIGN(1.d0,Y2).LT.ZERO) GOTO 20  ! (Y1*Y2.LT.ZERO)
5416          X1 = X2
5417          Y1 = Y2
5418 10    CONTINUE
5420 ! *** { YLO, YHI } > 0.0 THE SOLUTION IS ALWAYS SUPERSATURATED WITH NH4CL
5422       YLO= Y1                      ! Save Y-value at Hi position
5423       IF (YLO.GT.ZERO .AND. YHI.GT.ZERO) THEN
5424          Y3 = FUNCI4A (ZERO)
5425          GOTO 50
5426       ELSE IF (ABS(Y2) .LT. EPS) THEN   ! X2 IS A SOLUTION
5427          GOTO 50
5428       ELSE
5429          CALL PUSHERR (0001, 'CALCI4')    ! WARNING ERROR: NO SOLUTION
5430          GOTO 50
5431       ENDIF
5433 ! *** PERFORM BISECTION ***********************************************
5435 20    DO 30 I=1,MAXIT
5436          X3 = 0.5*(X1+X2)
5437          Y3 = FUNCI4A (X3)
5438          IF (SIGN(1.d0,Y1)*SIGN(1.d0,Y3) .LE. ZERO) THEN  ! (Y1*Y3 .LE. ZERO)
5439             Y2    = Y3
5440             X2    = X3
5441          ELSE
5442             Y1    = Y3
5443             X1    = X3
5444          ENDIF
5445          IF (ABS(X2-X1) .LE. EPS*X1) GOTO 40
5446 30    CONTINUE
5447       CALL PUSHERR (0002, 'CALCI4')    ! WARNING ERROR: NO CONVERGENCE
5449 ! *** CONVERGED ; RETURN **********************************************
5451 40    X3 = 0.5*(X1+X2)
5452       Y3 = FUNCI4A (X3)
5454 50    RETURN
5456 ! *** END OF SUBROUTINE CALCI4 *****************************************
5458       END SUBROUTINE CALCI4
5463 !=======================================================================
5465 ! *** ISORROPIA CODE
5466 ! *** SUBROUTINE FUNCI4A
5467 ! *** CASE I4
5469 !     THE MAIN CHARACTERISTICS OF THIS REGIME ARE:
5470 !     1. SULFATE RICH, NO FREE ACID (1.0 <= SULRAT < 2.0)
5471 !     2. SOLID & LIQUID AEROSOL POSSIBLE
5472 !     3. SOLIDS POSSIBLE : (NH4)2SO4, NA2SO4
5474 ! *** COPYRIGHT 1996-2006, UNIVERSITY OF MIAMI, CARNEGIE MELLON UNIVERSITY,
5475 ! *** GEORGIA INSTITUTE OF TECHNOLOGY
5476 ! *** WRITTEN BY ATHANASIOS NENES
5477 ! *** UPDATED BY CHRISTOS FOUNTOUKIS
5479 !=======================================================================
5481       DOUBLE PRECISION FUNCTION FUNCI4A (P4)
5482       USE ISRPIA
5483       USE SOLUT
5484       IMPLICIT DOUBLE PRECISION (A-H, O-Z)
5485 !      INCLUDE 'ISRPIA.INC'
5486 !      COMMON /SOLUT/ CHI1, CHI2, CHI3, CHI4, CHI5, CHI6, CHI7, CHI8,   &
5487 !                     PSI1, PSI2, PSI3, PSI4, PSI5, PSI6, PSI7, PSI8,   &
5488 !                     A1,   A2,   A3,   A4,   A5,   A6,   A7,   A8
5490 ! *** SETUP PARAMETERS ************************************************
5492       PSI4   = P4     ! PSI3 already assigned in FUNCI4A
5493       FRST   = .TRUE.
5494       CALAIN = .TRUE.
5496 ! *** SOLVE EQUATIONS ; WITH ITERATIONS FOR ACTIVITY COEF. ************
5498       DO 10 I=1,NSWEEP
5500       A4 = XK5 *(WATER/GAMA(2))**3.0
5501       A5 = XK7 *(WATER/GAMA(4))**3.0
5502       A6 = XK1 *WATER/GAMA(7)*(GAMA(8)/GAMA(7))**2.
5503       A7 = SQRT(A4/A5)
5505 !  CALCULATE DISSOCIATION QUANTITIES
5507       BB   = PSI2 + PSI4 + PSI5 + A6                    ! PSI6
5508       CC   =-A6*(PSI2 + PSI3 + PSI1)
5509       DD   = BB*BB - 4.D0*CC
5510       PSI6 = 0.5D0*(-BB + SQRT(DD))
5512       PSI5 = (PSI3 + 2.D0*PSI4 - A7*(3.D0*PSI2 + PSI1))/2.D0/A7
5513       PSI5 = MIN (PSI5, CHI5)
5515 ! *** CALCULATE SPECIATION ********************************************
5517       MOLAL (1) = PSI6                            ! HI
5518       MOLAL (2) = 2.D0*PSI4 + PSI3                ! NAI
5519       MOLAL (3) = 3.D0*PSI2 + 2.D0*PSI5 + PSI1    ! NH4I
5520       MOLAL (5) = PSI2 + PSI4 + PSI5 + PSI6       ! SO4I
5521       MOLAL (6) = PSI2 + PSI3 + PSI1 - PSI6       ! HSO4I
5522       CLC       = ZERO
5523       CNAHSO4   = ZERO
5524       CNA2SO4   = CHI4 - PSI4
5525       CNH42S4   = CHI5 - PSI5
5526       CNH4HS4   = ZERO
5527       CALL CALCMR                                 ! Water content
5529 ! *** CALCULATE ACTIVITIES OR TERMINATE INTERNAL LOOP *****************
5531       IF (FRST.AND.CALAOU .OR. .NOT.FRST.AND.CALAIN) THEN
5532          CALL CALCACT
5533       ELSE
5534          GOTO 20
5535       ENDIF
5536 10    CONTINUE
5538 ! *** CALCULATE OBJECTIVE FUNCTION ************************************
5540 20    A4     = XK5 *(WATER/GAMA(2))**3.0
5541       FUNCI4A= MOLAL(5)*MOLAL(2)*MOLAL(2)/A4 - ONE
5542       RETURN
5544 ! *** END OF FUNCTION FUNCI4A ********************************************
5546       END FUNCTION FUNCI4A     
5547 !=======================================================================
5549 ! *** ISORROPIA CODE
5550 ! *** SUBROUTINE CALCI3
5551 ! *** CASE I3
5553 !     THE MAIN CHARACTERISTICS OF THIS REGIME ARE:
5554 !     1. SULFATE RICH, NO FREE ACID (1.0 <= SULRAT < 2.0)
5555 !     2. SOLID & LIQUID AEROSOL POSSIBLE
5556 !     3. SOLIDS POSSIBLE : (NH4)2SO4, NA2SO4, NAHSO4, LC
5558 !     THERE ARE THREE REGIMES IN THIS CASE:
5559 !     1.(NA,NH4)HSO4(s) POSSIBLE. LIQUID & SOLID AEROSOL (SUBROUTINE CALCI3A)
5560 !     2.(NA,NH4)HSO4(s) NOT POSSIBLE, AND RH < MDRH. SOLID AEROSOL ONLY
5561 !     3.(NA,NH4)HSO4(s) NOT POSSIBLE, AND RH >= MDRH. SOLID & LIQUID AEROSOL
5563 !     REGIMES 2. AND 3. ARE CONSIDERED TO BE THE SAME AS CASES I1A, I2B
5564 !     RESPECTIVELY
5566 ! *** COPYRIGHT 1996-2006, UNIVERSITY OF MIAMI, CARNEGIE MELLON UNIVERSITY,
5567 ! *** GEORGIA INSTITUTE OF TECHNOLOGY
5568 ! *** WRITTEN BY ATHANASIOS NENES
5569 ! *** UPDATED BY CHRISTOS FOUNTOUKIS
5571 !=======================================================================
5573       SUBROUTINE CALCI3
5574       USE ISRPIA
5575       IMPLICIT DOUBLE PRECISION (A-H, O-Z)
5576 !      INCLUDE 'ISRPIA.INC'
5577       EXTERNAL CALCI1A, CALCI4
5579 ! *** FIND DRY COMPOSITION **********************************************
5581       CALL CALCI1A
5583 ! *** REGIME DEPENDS UPON THE POSSIBLE SOLIDS & RH **********************
5585       IF (CNH4HS4.GT.TINY .OR. CNAHSO4.GT.TINY) THEN
5586          SCASE = 'I3 ; SUBCASE 1'
5587          CALL CALCI3A                     ! FULL SOLUTION
5588          SCASE = 'I3 ; SUBCASE 1'
5589       ENDIF
5591       IF (WATER.LE.TINY) THEN
5592          IF (RH.LT.DRMI3) THEN         ! SOLID SOLUTION
5593             WATER = TINY
5594             DO 10 I=1,NIONS
5595                MOLAL(I) = ZERO
5596 10          CONTINUE
5597             CALL CALCI1A
5598             SCASE = 'I3 ; SUBCASE 2'
5600          ELSEIF (RH.GE.DRMI3) THEN     ! MDRH OF I3
5601             SCASE = 'I3 ; SUBCASE 3'
5602             CALL CALCMDRH (RH, DRMI3, DRLC, CALCI1A, CALCI4)
5603             SCASE = 'I3 ; SUBCASE 3'
5604          ENDIF
5605       ENDIF
5607       RETURN
5609 ! *** END OF SUBROUTINE CALCI3 ******************************************
5611       END SUBROUTINE CALCI3
5615 !=======================================================================
5617 ! *** ISORROPIA CODE
5618 ! *** SUBROUTINE CALCI3A
5619 ! *** CASE I3 ; SUBCASE 1
5621 !     THE MAIN CHARACTERISTICS OF THIS REGIME ARE:
5622 !     1. SULFATE RICH, NO FREE ACID (1.0 <= SULRAT < 2.0)
5623 !     2. SOLID & LIQUID AEROSOL POSSIBLE
5624 !     3. SOLIDS POSSIBLE : (NH4)2SO4, NA2SO4, LC
5626 ! *** COPYRIGHT 1996-2006, UNIVERSITY OF MIAMI, CARNEGIE MELLON UNIVERSITY,
5627 ! *** GEORGIA INSTITUTE OF TECHNOLOGY
5628 ! *** WRITTEN BY ATHANASIOS NENES
5629 ! *** UPDATED BY CHRISTOS FOUNTOUKIS
5631 !=======================================================================
5633       SUBROUTINE CALCI3A
5634       USE ISRPIA
5635       USE SOLUT
5636       IMPLICIT DOUBLE PRECISION (A-H, O-Z)
5637 !      INCLUDE 'ISRPIA.INC'
5638 !      COMMON /SOLUT/ CHI1, CHI2, CHI3, CHI4, CHI5, CHI6, CHI7, CHI8,   &
5639 !                     PSI1, PSI2, PSI3, PSI4, PSI5, PSI6, PSI7, PSI8,   &
5640 !                     A1,   A2,   A3,   A4,   A5,   A6,   A7,   A8
5642 ! *** FIND DRY COMPOSITION **********************************************
5644       CALL CALCI1A         ! Needed when called from CALCMDRH
5646 ! *** SETUP PARAMETERS ************************************************
5648       CHI1 = CNH4HS4               ! Save from CALCI1 run
5649       CHI2 = CLC
5650       CHI3 = CNAHSO4
5651       CHI4 = CNA2SO4
5652       CHI5 = CNH42S4
5654       PSI1 = CNH4HS4               ! ASSIGN INITIAL PSI's
5655       PSI2 = ZERO
5656       PSI3 = CNAHSO4
5657       PSI4 = ZERO
5658       PSI5 = ZERO
5660       CALAOU = .TRUE.              ! Outer loop activity calculation flag
5661       PSI2LO = ZERO                ! Low  limit
5662       PSI2HI = CHI2                ! High limit
5664 ! *** INITIAL VALUES FOR BISECTION ************************************
5666       X1 = PSI2HI
5667       Y1 = FUNCI3A (X1)
5668       YHI= Y1                      ! Save Y-value at HI position
5670 ! *** YHI < 0.0 THE SOLUTION IS ALWAYS UNDERSATURATED WITH LC *********
5672       IF (YHI.LT.EPS) GOTO 50
5674 ! *** ROOT TRACKING ; FOR THE RANGE OF HI AND LO **********************
5676       DX = (PSI2HI-PSI2LO)/REAL(NDIV)
5677       DO 10 I=1,NDIV
5678          X2 = MAX(X1-DX, PSI2LO)
5679          Y2 = FUNCI3A (X2)
5680          IF (SIGN(1.d0,Y1)*SIGN(1.d0,Y2).LT.ZERO) GOTO 20  ! (Y1*Y2.LT.ZERO)
5681          X1 = X2
5682          Y1 = Y2
5683 10    CONTINUE
5685 ! *** { YLO, YHI } > 0.0 THE SOLUTION IS ALWAYS SUPERSATURATED WITH LC
5687       IF (Y2.GT.EPS) Y2 = FUNCI3A (ZERO)
5688       GOTO 50
5690 ! *** PERFORM BISECTION ***********************************************
5692 20    DO 30 I=1,MAXIT
5693          X3 = 0.5*(X1+X2)
5694          Y3 = FUNCI3A (X3)
5695          IF (SIGN(1.d0,Y1)*SIGN(1.d0,Y3) .LE. ZERO) THEN  ! (Y1*Y3 .LE. ZERO)
5696             Y2    = Y3
5697             X2    = X3
5698          ELSE
5699             Y1    = Y3
5700             X1    = X3
5701          ENDIF
5702          IF (ABS(X2-X1) .LE. EPS*X1) GOTO 40
5703 30    CONTINUE
5704       CALL PUSHERR (0002, 'CALCI3A')    ! WARNING ERROR: NO CONVERGENCE
5706 ! *** CONVERGED ; RETURN **********************************************
5708 40    X3 = 0.5*(X1+X2)
5709       Y3 = FUNCI3A (X3)
5711 50    RETURN
5713 ! *** END OF SUBROUTINE CALCI3A *****************************************
5715       END SUBROUTINE CALCI3A
5717 !=======================================================================
5719 ! *** ISORROPIA CODE
5720 ! *** SUBROUTINE FUNCI3A
5721 ! *** CASE I3 ; SUBCASE 1
5723 !     THE MAIN CHARACTERISTICS OF THIS REGIME ARE:
5724 !     1. SULFATE RICH, NO FREE ACID (1.0 <= SULRAT < 2.0)
5725 !     2. SOLID & LIQUID AEROSOL POSSIBLE
5726 !     3. SOLIDS POSSIBLE : (NH4)2SO4, NA2SO4, LC
5728 ! *** COPYRIGHT 1996-2006, UNIVERSITY OF MIAMI, CARNEGIE MELLON UNIVERSITY,
5729 ! *** GEORGIA INSTITUTE OF TECHNOLOGY
5730 ! *** WRITTEN BY ATHANASIOS NENES
5731 ! *** UPDATED BY CHRISTOS FOUNTOUKIS
5733 !=======================================================================
5735       DOUBLE PRECISION FUNCTION FUNCI3A (P2)
5736       USE ISRPIA
5737       USE SOLUT
5738       IMPLICIT DOUBLE PRECISION (A-H, O-Z)
5739 !      INCLUDE 'ISRPIA.INC'
5740 !      COMMON /SOLUT/ CHI1, CHI2, CHI3, CHI4, CHI5, CHI6, CHI7, CHI8,   &
5741 !                     PSI1, PSI2, PSI3, PSI4, PSI5, PSI6, PSI7, PSI8,   &
5742 !                     A1,   A2,   A3,   A4,   A5,   A6,   A7,   A8
5744 ! *** SETUP PARAMETERS ************************************************
5746       PSI2   = P2                  ! Save PSI2 in COMMON BLOCK
5747       PSI4LO = ZERO                ! Low  limit for PSI4
5748       PSI4HI = CHI4                ! High limit for PSI4
5750 ! *** IF NH3 =0, CALL FUNCI3B FOR Y4=0 ********************************
5752       IF (CHI4.LE.TINY) THEN
5753          FUNCI3A = FUNCI3B (ZERO)
5754          GOTO 50
5755       ENDIF
5757 ! *** INITIAL VALUES FOR BISECTION ************************************
5759       X1 = PSI4HI
5760       Y1 = FUNCI3B (X1)
5761       IF (ABS(Y1).LE.EPS) GOTO 50
5762       YHI= Y1                      ! Save Y-value at HI position
5764 ! *** YHI < 0.0 THE SOLUTION IS ALWAYS UNDERSATURATED WITH NA2SO4 *****
5766       IF (YHI.LT.ZERO) GOTO 50
5768 ! *** ROOT TRACKING ; FOR THE RANGE OF HI AND LO **********************
5770       DX = (PSI4HI-PSI4LO)/REAL(NDIV)
5771       DO 10 I=1,NDIV
5772          X2 = MAX(X1-DX, PSI4LO)
5773          Y2 = FUNCI3B (X2)
5774          IF (SIGN(1.d0,Y1)*SIGN(1.d0,Y2).LT.ZERO) GOTO 20  ! (Y1*Y2.LT.ZERO)
5775          X1 = X2
5776          Y1 = Y2
5777 10    CONTINUE
5779 ! *** { YLO, YHI } > 0.0 THE SOLUTION IS ALWAYS SUPERSATURATED WITH NA2SO4
5781       IF (Y2.GT.EPS) Y2 = FUNCI3B (PSI4LO)
5782       GOTO 50
5784 ! *** PERFORM BISECTION ***********************************************
5786 20    DO 30 I=1,MAXIT
5787          X3 = 0.5*(X1+X2)
5788          Y3 = FUNCI3B (X3)
5789          IF (SIGN(1.d0,Y1)*SIGN(1.d0,Y3) .LE. ZERO) THEN  ! (Y1*Y3 .LE. ZERO)
5790             Y2    = Y3
5791             X2    = X3
5792          ELSE
5793             Y1    = Y3
5794             X1    = X3
5795          ENDIF
5796          IF (ABS(X2-X1) .LE. EPS*X1) GOTO 40
5797 30    CONTINUE
5798       CALL PUSHERR (0004, 'FUNCI3A')    ! WARNING ERROR: NO CONVERGENCE
5800 ! *** INNER LOOP CONVERGED **********************************************
5802 40    X3 = 0.5*(X1+X2)
5803       Y3 = FUNCI3B (X3)
5805 ! *** CALCULATE FUNCTION VALUE FOR OUTER LOOP ***************************
5807 50    A2      = XK13*(WATER/GAMA(13))**5.0
5808       FUNCI3A = MOLAL(5)*MOLAL(6)*MOLAL(3)**3.D0/A2 - ONE
5809       RETURN
5811 ! *** END OF FUNCTION FUNCI3A *******************************************
5813       END FUNCTION FUNCI3A     
5817 !=======================================================================
5819 ! *** ISORROPIA CODE
5820 ! *** FUNCTION FUNCI3B
5821 ! *** CASE I3 ; SUBCASE 2
5823 !     THE MAIN CHARACTERISTICS OF THIS REGIME ARE:
5824 !     1. SULFATE RICH, NO FREE ACID (1.0 <= SULRAT < 2.0)
5825 !     2. SOLID & LIQUID AEROSOL POSSIBLE
5826 !     3. SOLIDS POSSIBLE : (NH4)2SO4, NA2SO4, LC
5828 !     SOLUTION IS SAVED IN COMMON BLOCK /CASE/
5830 !=======================================================================
5832       DOUBLE PRECISION FUNCTION FUNCI3B (P4)
5833       USE ISRPIA
5834       USE SOLUT
5835       IMPLICIT DOUBLE PRECISION (A-H, O-Z)
5836 !      INCLUDE 'ISRPIA.INC'
5838 !      COMMON /SOLUT/ CHI1, CHI2, CHI3, CHI4, CHI5, CHI6, CHI7, CHI8,   &
5839 !                     PSI1, PSI2, PSI3, PSI4, PSI5, PSI6, PSI7, PSI8,   &
5840 !                     A1,   A2,   A3,   A4,   A5,   A6,   A7,   A8
5842 ! *** SETUP PARAMETERS ************************************************
5844       PSI4   = P4
5846 ! *** SETUP PARAMETERS ************************************************
5848       FRST   = .TRUE.
5849       CALAIN = .TRUE.
5851 ! *** SOLVE EQUATIONS ; WITH ITERATIONS FOR ACTIVITY COEF. ************
5853       DO 10 I=1,NSWEEP
5855       A4 = XK5*(WATER/GAMA(2))**3.0
5856       A5 = XK7*(WATER/GAMA(4))**3.0
5857       A6 = XK1*WATER/GAMA(7)*(GAMA(8)/GAMA(7))**2.
5858       A7 = SQRT(A4/A5)
5860 !  CALCULATE DISSOCIATION QUANTITIES
5862       BB   = PSI2 + PSI4 + PSI5 + A6                    ! PSI6
5863       CC   =-A6*(PSI2 + PSI3 + PSI1)
5864       DD   = BB*BB - 4.D0*CC
5865       PSI6 = 0.5D0*(-BB + SQRT(DD))
5867       PSI5 = (PSI3 + 2.D0*PSI4 - A7*(3.D0*PSI2 + PSI1))/2.D0/A7
5868       PSI5 = MIN (PSI5, CHI5)
5870 ! *** CALCULATE SPECIATION ********************************************
5872       MOLAL(1) = PSI6                                  ! HI
5873       MOLAL(2) = 2.D0*PSI4 + PSI3                      ! NAI
5874       MOLAL(3) = 3.D0*PSI2 + 2.D0*PSI5 + PSI1          ! NH4I
5875       MOLAL(5) = PSI2 + PSI4 + PSI5 + PSI6             ! SO4I
5876       MOLAL(6) = MAX(PSI2 + PSI3 + PSI1 - PSI6, TINY)  ! HSO4I
5877       CLC      = MAX(CHI2 - PSI2, ZERO)
5878       CNAHSO4  = ZERO
5879       CNA2SO4  = MAX(CHI4 - PSI4, ZERO)
5880       CNH42S4  = MAX(CHI5 - PSI5, ZERO)
5881       CNH4HS4  = ZERO
5882       CALL CALCMR                                       ! Water content
5884 ! *** CALCULATE ACTIVITIES OR TERMINATE INTERNAL LOOP *****************
5886       IF (FRST.AND.CALAOU .OR. .NOT.FRST.AND.CALAIN) THEN
5887          CALL CALCACT
5888       ELSE
5889          GOTO 20
5890       ENDIF
5891 10    CONTINUE
5893 ! *** CALCULATE OBJECTIVE FUNCTION ************************************
5895 20    A4     = XK5 *(WATER/GAMA(2))**3.0
5896       FUNCI3B= MOLAL(5)*MOLAL(2)*MOLAL(2)/A4 - ONE
5897       RETURN
5899 ! *** END OF FUNCTION FUNCI3B ********************************************
5901       END FUNCTION FUNCI3B     
5902 !=======================================================================
5904 ! *** ISORROPIA CODE
5905 ! *** SUBROUTINE CALCI2
5906 ! *** CASE I2
5908 !     THE MAIN CHARACTERISTICS OF THIS REGIME ARE:
5909 !     1. SULFATE RICH, NO FREE ACID (1.0 <= SULRAT < 2.0)
5910 !     2. SOLID & LIQUID AEROSOL POSSIBLE
5911 !     3. SOLIDS POSSIBLE : (NH4)2SO4, NA2SO4, NAHSO4, LC
5913 !     THERE ARE THREE REGIMES IN THIS CASE:
5914 !     1. NH4HSO4(s) POSSIBLE. LIQUID & SOLID AEROSOL (SUBROUTINE CALCI2A)
5915 !     2. NH4HSO4(s) NOT POSSIBLE, AND RH < MDRH. SOLID AEROSOL ONLY
5916 !     3. NH4HSO4(s) NOT POSSIBLE, AND RH >= MDRH. SOLID & LIQUID AEROSOL
5918 !     REGIMES 2. AND 3. ARE CONSIDERED TO BE THE SAME AS CASES I1A, I2B
5919 !     RESPECTIVELY
5921 ! *** COPYRIGHT 1996-2006, UNIVERSITY OF MIAMI, CARNEGIE MELLON UNIVERSITY,
5922 ! *** GEORGIA INSTITUTE OF TECHNOLOGY
5923 ! *** WRITTEN BY ATHANASIOS NENES
5924 ! *** UPDATED BY CHRISTOS FOUNTOUKIS
5926 !=======================================================================
5928       SUBROUTINE CALCI2
5929       USE ISRPIA
5930       IMPLICIT DOUBLE PRECISION (A-H, O-Z)
5931 !      INCLUDE 'ISRPIA.INC'
5932       EXTERNAL CALCI1A, CALCI3A
5934 ! *** FIND DRY COMPOSITION **********************************************
5936       CALL CALCI1A
5938 ! *** REGIME DEPENDS UPON THE POSSIBLE SOLIDS & RH **********************
5940       IF (CNH4HS4.GT.TINY) THEN
5941          SCASE = 'I2 ; SUBCASE 1'
5942          CALL CALCI2A
5943          SCASE = 'I2 ; SUBCASE 1'
5944       ENDIF
5946       IF (WATER.LE.TINY) THEN
5947          IF (RH.LT.DRMI2) THEN         ! SOLID SOLUTION ONLY
5948             WATER = TINY
5949             DO 10 I=1,NIONS
5950                MOLAL(I) = ZERO
5951 10          CONTINUE
5952             CALL CALCI1A
5953             SCASE = 'I2 ; SUBCASE 2'
5955          ELSEIF (RH.GE.DRMI2) THEN     ! MDRH OF I2
5956             SCASE = 'I2 ; SUBCASE 3'
5957             CALL CALCMDRH (RH, DRMI2, DRNAHSO4, CALCI1A, CALCI3A)
5958             SCASE = 'I2 ; SUBCASE 3'
5959          ENDIF
5960       ENDIF
5962       RETURN
5964 ! *** END OF SUBROUTINE CALCI2 ******************************************
5966       END SUBROUTINE CALCI2
5969 !=======================================================================
5971 ! *** ISORROPIA CODE
5972 ! *** SUBROUTINE CALCI2A
5973 ! *** CASE I2 ; SUBCASE A
5975 !     THE MAIN CHARACTERISTICS OF THIS REGIME ARE:
5976 !     1. SULFATE RICH, NO FREE ACID (1.0 <= SULRAT < 2.0)
5977 !     2. SOLID & LIQUID AEROSOL POSSIBLE
5978 !     3. SOLIDS POSSIBLE : (NH4)2SO4, NA2SO4, NAHSO4, LC
5980 ! *** COPYRIGHT 1996-2006, UNIVERSITY OF MIAMI, CARNEGIE MELLON UNIVERSITY,
5981 ! *** GEORGIA INSTITUTE OF TECHNOLOGY
5982 ! *** WRITTEN BY ATHANASIOS NENES
5983 ! *** UPDATED BY CHRISTOS FOUNTOUKIS
5985 !=======================================================================
5987       SUBROUTINE CALCI2A
5988       USE ISRPIA
5989       USE SOLUT
5990       IMPLICIT DOUBLE PRECISION (A-H, O-Z)
5991 !      INCLUDE 'ISRPIA.INC'
5992 !      COMMON /SOLUT/ CHI1, CHI2, CHI3, CHI4, CHI5, CHI6, CHI7, CHI8,   &
5993 !                     PSI1, PSI2, PSI3, PSI4, PSI5, PSI6, PSI7, PSI8,   &
5994 !                     A1,   A2,   A3,   A4,   A5,   A6,   A7,   A8
5996 ! *** FIND DRY COMPOSITION **********************************************
5998       CALL CALCI1A    ! Needed when called from CALCMDRH
6000 ! *** SETUP PARAMETERS ************************************************
6002       CHI1 = CNH4HS4               ! Save from CALCI1 run
6003       CHI2 = CLC
6004       CHI3 = CNAHSO4
6005       CHI4 = CNA2SO4
6006       CHI5 = CNH42S4
6008       PSI1 = CNH4HS4               ! ASSIGN INITIAL PSI's
6009       PSI2 = ZERO
6010       PSI3 = ZERO
6011       PSI4 = ZERO
6012       PSI5 = ZERO
6014       CALAOU = .TRUE.              ! Outer loop activity calculation flag
6015       PSI2LO = ZERO                ! Low  limit
6016       PSI2HI = CHI2                ! High limit
6018 ! *** INITIAL VALUES FOR BISECTION ************************************
6020       X1 = PSI2HI
6021       Y1 = FUNCI2A (X1)
6022       YHI= Y1                      ! Save Y-value at HI position
6024 ! *** YHI < 0.0 THE SOLUTION IS ALWAYS UNDERSATURATED WITH LC *********
6026       IF (YHI.LT.EPS) GOTO 50
6028 ! *** ROOT TRACKING ; FOR THE RANGE OF HI AND LO **********************
6030       DX = (PSI2HI-PSI2LO)/REAL(NDIV)
6031       DO 10 I=1,NDIV
6032          X2 = MAX(X1-DX, PSI2LO)
6033          Y2 = FUNCI2A (X2)
6034          IF (SIGN(1.d0,Y1)*SIGN(1.d0,Y2).LT.ZERO) GOTO 20  ! (Y1*Y2.LT.ZERO)
6035          X1 = X2
6036          Y1 = Y2
6037 10    CONTINUE
6039 ! *** { YLO, YHI } > 0.0 THE SOLUTION IS ALWAYS SUPERSATURATED WITH LC
6041       IF (Y2.GT.EPS) Y2 = FUNCI3A (ZERO)
6042       GOTO 50
6044 ! *** PERFORM BISECTION ***********************************************
6046 20    DO 30 I=1,MAXIT
6047          X3 = 0.5*(X1+X2)
6048          Y3 = FUNCI2A (X3)
6049          IF (SIGN(1.d0,Y1)*SIGN(1.d0,Y3) .LE. ZERO) THEN  ! (Y1*Y3 .LE. ZERO)
6050             Y2    = Y3
6051             X2    = X3
6052          ELSE
6053             Y1    = Y3
6054             X1    = X3
6055          ENDIF
6056          IF (ABS(X2-X1) .LE. EPS*X1) GOTO 40
6057 30    CONTINUE
6058       CALL PUSHERR (0002, 'CALCI2A')    ! WARNING ERROR: NO CONVERGENCE
6060 ! *** CONVERGED ; RETURN **********************************************
6062 40    X3 = 0.5*(X1+X2)
6063       Y3 = FUNCI2A (X3)
6065 50    RETURN
6067 ! *** END OF SUBROUTINE CALCI2A *****************************************
6069       END SUBROUTINE CALCI2A
6074 !=======================================================================
6076 ! *** ISORROPIA CODE
6077 ! *** SUBROUTINE FUNCI2A
6078 ! *** CASE I2 ; SUBCASE 1
6080 !     THE MAIN CHARACTERISTICS OF THIS REGIME ARE:
6081 !     1. SULFATE RICH, NO FREE ACID (1.0 <= SULRAT < 2.0)
6082 !     2. SOLID & LIQUID AEROSOL POSSIBLE
6083 !     3. SOLIDS POSSIBLE : (NH4)2SO4, NA2SO4, NAHSO4, LC
6085 ! *** COPYRIGHT 1996-2006, UNIVERSITY OF MIAMI, CARNEGIE MELLON UNIVERSITY,
6086 ! *** GEORGIA INSTITUTE OF TECHNOLOGY
6087 ! *** WRITTEN BY ATHANASIOS NENES
6088 ! *** UPDATED BY CHRISTOS FOUNTOUKIS
6090 !=======================================================================
6092       DOUBLE PRECISION FUNCTION FUNCI2A (P2)
6093       USE ISRPIA
6094       USE SOLUT
6095       IMPLICIT DOUBLE PRECISION (A-H, O-Z)
6096 !      INCLUDE 'ISRPIA.INC'
6097 !      COMMON /SOLUT/ CHI1, CHI2, CHI3, CHI4, CHI5, CHI6, CHI7, CHI8,   &
6098 !                     PSI1, PSI2, PSI3, PSI4, PSI5, PSI6, PSI7, PSI8,   &
6099 !                     A1,   A2,   A3,   A4,   A5,   A6,   A7,   A8
6101 ! *** SETUP PARAMETERS ************************************************
6103       FRST   = .TRUE.
6104       CALAIN = .TRUE.
6105       PSI2   = P2                  ! Save PSI2 in COMMON BLOCK
6106       PSI3   = CHI3
6107       PSI4   = CHI4
6108       PSI5   = CHI5
6109       PSI6   = ZERO
6111 ! *** SOLVE EQUATIONS ; WITH ITERATIONS FOR ACTIVITY COEF. ************
6113       DO 10 I=1,NSWEEP
6115       A3 = XK11*(WATER/GAMA(12))**2.0
6116       A4 = XK5 *(WATER/GAMA(2))**3.0
6117       A5 = XK7 *(WATER/GAMA(4))**3.0
6118       A6 = XK1 *WATER/GAMA(7)*(GAMA(8)/GAMA(7))**2.
6119       A7 = SQRT(A4/A5)
6121 !  CALCULATE DISSOCIATION QUANTITIES
6123       IF (CHI5.GT.TINY .AND. WATER.GT.TINY) THEN
6124          PSI5 = (PSI3 + 2.D0*PSI4 - A7*(3.D0*PSI2 + PSI1))/2.D0/A7
6125          PSI5 = MAX(MIN (PSI5, CHI5), TINY)
6126       ENDIF
6128       IF (CHI4.GT.TINY .AND. WATER.GT.TINY) THEN
6129          AA   = PSI2+PSI5+PSI6+PSI3
6130          BB   = PSI3*AA
6131          CC   = 0.25D0*(PSI3*PSI3*(PSI2+PSI5+PSI6)-A4)
6132          CALL POLY3 (AA, BB, CC, PSI4, ISLV)
6133          IF (ISLV.EQ.0) THEN
6134             PSI4 = MIN (PSI4, CHI4)
6135          ELSE
6136             PSI4 = ZERO
6137          ENDIF
6138       ENDIF
6140       IF (CHI3.GT.TINY .AND. WATER.GT.TINY) THEN
6141          AA   = 2.D0*PSI4 + PSI2 + PSI1 - PSI6
6142          BB   = 2.D0*PSI4*(PSI2 + PSI1 - PSI6) - A3
6143          CC   = ZERO
6144          CALL POLY3 (AA, BB, CC, PSI3, ISLV)
6145          IF (ISLV.EQ.0) THEN
6146             PSI3 = MIN (PSI3, CHI3)
6147          ELSE
6148             PSI3 = ZERO
6149          ENDIF
6150       ENDIF
6152       BB   = PSI2 + PSI4 + PSI5 + A6                    ! PSI6
6153       CC   =-A6*(PSI2 + PSI3 + PSI1)
6154       DD   = BB*BB - 4.D0*CC
6155       PSI6 = 0.5D0*(-BB + SQRT(DD))
6157 ! *** CALCULATE SPECIATION ********************************************
6159       MOLAL (1) = PSI6                           ! HI
6160       MOLAL (2) = 2.D0*PSI4 + PSI3               ! NAI
6161       MOLAL (3) = 3.D0*PSI2 + 2.D0*PSI5 + PSI1   ! NH4I
6162       MOLAL (5) = PSI2 + PSI4 + PSI5 + PSI6      ! SO4I
6163       MOLAL (6) = PSI2 + PSI3 + PSI1 - PSI6      ! HSO4I
6164       CLC       = CHI2 - PSI2
6165       CNAHSO4   = CHI3 - PSI3
6166       CNA2SO4   = CHI4 - PSI4
6167       CNH42S4   = CHI5 - PSI5
6168       CNH4HS4   = ZERO
6169       CALL CALCMR                                ! Water content
6171 ! *** CALCULATE ACTIVITIES OR TERMINATE INTERNAL LOOP *****************
6173       IF (FRST.AND.CALAOU .OR. .NOT.FRST.AND.CALAIN) THEN
6174          CALL CALCACT
6175       ELSE
6176          GOTO 20
6177       ENDIF
6178 10    CONTINUE
6180 ! *** CALCULATE FUNCTION VALUE FOR OUTER LOOP ***************************
6182 20    A2      = XK13*(WATER/GAMA(13))**5.0
6183       FUNCI2A = MOLAL(5)*MOLAL(6)*MOLAL(3)**3.D0/A2 - ONE
6184       RETURN
6186 ! *** END OF FUNCTION FUNCI2A *******************************************
6188       END FUNCTION FUNCI2A     
6190 !=======================================================================
6192 ! *** ISORROPIA CODE
6193 ! *** SUBROUTINE CALCI1
6194 ! *** CASE I1
6196 !     THE MAIN CHARACTERISTICS OF THIS REGIME ARE:
6197 !     1. SULFATE RICH, NO FREE ACID (1.0 <= SULRAT < 2.0)
6198 !     2. SOLID AEROSOL ONLY
6199 !     3. SOLIDS POSSIBLE : NH4NO3, NH4CL, NA2SO4
6201 !     THERE ARE TWO POSSIBLE REGIMES HERE, DEPENDING ON RELATIVE HUMIDITY:
6202 !     1. WHEN RH >= MDRH ; LIQUID PHASE POSSIBLE (MDRH REGION)
6203 !     2. WHEN RH < MDRH  ; ONLY SOLID PHASE POSSIBLE (SUBROUTINE CALCI1A)
6205 ! *** COPYRIGHT 1996-2006, UNIVERSITY OF MIAMI, CARNEGIE MELLON UNIVERSITY,
6206 ! *** GEORGIA INSTITUTE OF TECHNOLOGY
6207 ! *** WRITTEN BY ATHANASIOS NENES
6208 ! *** UPDATED BY CHRISTOS FOUNTOUKIS
6210 !=======================================================================
6212       SUBROUTINE CALCI1
6213       USE ISRPIA
6214       IMPLICIT DOUBLE PRECISION (A-H, O-Z)
6215 !      INCLUDE 'ISRPIA.INC'
6216       EXTERNAL CALCI1A, CALCI2A
6218 ! *** REGIME DEPENDS UPON THE AMBIENT RELATIVE HUMIDITY *****************
6220       IF (RH.LT.DRMI1) THEN
6221          SCASE = 'I1 ; SUBCASE 1'
6222          CALL CALCI1A              ! SOLID PHASE ONLY POSSIBLE
6223          SCASE = 'I1 ; SUBCASE 1'
6224       ELSE
6225          SCASE = 'I1 ; SUBCASE 2'  ! LIQUID & SOLID PHASE POSSIBLE
6226          CALL CALCMDRH (RH, DRMI1, DRNH4HS4, CALCI1A, CALCI2A)
6227          SCASE = 'I1 ; SUBCASE 2'
6228       ENDIF
6230 ! *** AMMONIA IN GAS PHASE **********************************************
6232 !      CALL CALCNH3
6234       RETURN
6236 ! *** END OF SUBROUTINE CALCI1 ******************************************
6238       END SUBROUTINE CALCI1
6241 !=======================================================================
6243 ! *** ISORROPIA CODE
6244 ! *** SUBROUTINE CALCI1A
6245 ! *** CASE I1 ; SUBCASE 1
6247 !     THE MAIN CHARACTERISTICS OF THIS REGIME ARE:
6248 !     1. SULFATE RICH, NO FREE ACID (1.0 <= SULRAT < 2.0)
6249 !     2. SOLID AEROSOL ONLY
6250 !     3. SOLIDS POSSIBLE : NH4HSO4, NAHSO4, (NH4)2SO4, NA2SO4, LC
6252 ! *** COPYRIGHT 1996-2006, UNIVERSITY OF MIAMI, CARNEGIE MELLON UNIVERSITY,
6253 ! *** GEORGIA INSTITUTE OF TECHNOLOGY
6254 ! *** WRITTEN BY ATHANASIOS NENES
6255 ! *** UPDATED BY CHRISTOS FOUNTOUKIS
6257 !=======================================================================
6259       SUBROUTINE CALCI1A
6260       USE ISRPIA
6261       IMPLICIT DOUBLE PRECISION (A-H, O-Z)
6262 !      INCLUDE 'ISRPIA.INC'
6264 ! *** CALCULATE NON VOLATILE SOLIDS ***********************************
6266       CNA2SO4 = 0.5D0*W(1)
6267       CNH4HS4 = ZERO
6268       CNAHSO4 = ZERO
6269       CNH42S4 = ZERO
6270       FRSO4   = MAX(W(2)-CNA2SO4, ZERO)
6272       CLC     = MIN(W(3)/3.D0, FRSO4/2.D0)
6273       FRSO4   = MAX(FRSO4-2.D0*CLC, ZERO)
6274       FRNH4   = MAX(W(3)-3.D0*CLC,  ZERO)
6276       IF (FRSO4.LE.TINY) THEN
6277          CLC     = MAX(CLC - FRNH4, ZERO)
6278          CNH42S4 = 2.D0*FRNH4
6280       ELSEIF (FRNH4.LE.TINY) THEN
6281          CNH4HS4 = 3.D0*MIN(FRSO4, CLC)
6282          CLC     = MAX(CLC-FRSO4, ZERO)
6283          IF (CNA2SO4.GT.TINY) THEN
6284             FRSO4   = MAX(FRSO4-CNH4HS4/3.D0, ZERO)
6285             CNAHSO4 = 2.D0*FRSO4
6286             CNA2SO4 = MAX(CNA2SO4-FRSO4, ZERO)
6287          ENDIF
6288       ENDIF
6290 ! *** CALCULATE GAS SPECIES *********************************************
6292       GHNO3 = W(4)
6293       GHCL  = W(5)
6294       GNH3  = ZERO
6296       RETURN
6298 ! *** END OF SUBROUTINE CALCI1A *****************************************
6300       END SUBROUTINE CALCI1A
6301 !=======================================================================
6303 ! *** ISORROPIA CODE
6304 ! *** SUBROUTINE CALCJ3
6305 ! *** CASE J3
6307 !     THE MAIN CHARACTERISTICS OF THIS REGIME ARE:
6308 !     1. SULFATE RICH, FREE ACID (SULRAT < 1.0)
6309 !     2. THERE IS ONLY A LIQUID PHASE
6311 ! *** COPYRIGHT 1996-2006, UNIVERSITY OF MIAMI, CARNEGIE MELLON UNIVERSITY,
6312 ! *** GEORGIA INSTITUTE OF TECHNOLOGY
6313 ! *** WRITTEN BY ATHANASIOS NENES
6314 ! *** UPDATED BY CHRISTOS FOUNTOUKIS
6316 !=======================================================================
6318       SUBROUTINE CALCJ3
6319       USE ISRPIA
6320       IMPLICIT DOUBLE PRECISION (A-H, O-Z)
6321 !      INCLUDE 'ISRPIA.INC'
6323       DOUBLE PRECISION LAMDA, KAPA
6325 ! *** SETUP PARAMETERS ************************************************
6327       CALAOU = .TRUE.              ! Outer loop activity calculation flag
6328       FRST   = .TRUE.
6329       CALAIN = .TRUE.
6331       LAMDA  = MAX(W(2) - W(3) - W(1), TINY)  ! FREE H2SO4
6332       CHI1   = W(1)                           ! NA TOTAL as NaHSO4
6333       CHI2   = W(3)                           ! NH4 TOTAL as NH4HSO4
6334       PSI1   = CHI1
6335       PSI2   = CHI2                           ! ALL NH4HSO4 DELIQUESCED
6337 ! *** SOLVE EQUATIONS ; WITH ITERATIONS FOR ACTIVITY COEF. ************
6339       DO 10 I=1,NSWEEP
6341       A3 = XK1  *WATER/GAMA(7)*(GAMA(8)/GAMA(7))**2.0
6343 !  CALCULATE DISSOCIATION QUANTITIES
6345       BB   = A3+LAMDA                        ! KAPA
6346       CC   =-A3*(LAMDA + PSI1 + PSI2)
6347       DD   = BB*BB-4.D0*CC
6348       KAPA = 0.5D0*(-BB+SQRT(DD))
6350 ! *** CALCULATE SPECIATION ********************************************
6352       MOLAL (1) = LAMDA + KAPA                 ! HI
6353       MOLAL (2) = PSI1                         ! NAI
6354       MOLAL (3) = PSI2                         ! NH4I
6355       MOLAL (4) = ZERO                         ! CLI
6356       MOLAL (5) = KAPA                         ! SO4I
6357       MOLAL (6) = LAMDA + PSI1 + PSI2 - KAPA   ! HSO4I
6358       MOLAL (7) = ZERO                         ! NO3I
6360       CNAHSO4   = ZERO
6361       CNH4HS4   = ZERO
6363       CALL CALCMR                              ! Water content
6365 ! *** CALCULATE ACTIVITIES OR TERMINATE INTERNAL LOOP *****************
6367       IF (FRST.AND.CALAOU .OR. .NOT.FRST.AND.CALAIN) THEN
6368          CALL CALCACT
6369       ELSE
6370          GOTO 50
6371       ENDIF
6372 10    CONTINUE
6374 50    RETURN
6376 ! *** END OF SUBROUTINE CALCJ3 ******************************************
6378       END SUBROUTINE CALCJ3
6379 !=======================================================================
6381 ! *** ISORROPIA CODE
6382 ! *** SUBROUTINE CALCJ2
6383 ! *** CASE J2
6385 !     THE MAIN CHARACTERISTICS OF THIS REGIME ARE:
6386 !     1. SULFATE RICH, FREE ACID (SULRAT < 1.0)
6387 !     2. THERE IS BOTH A LIQUID & SOLID PHASE
6388 !     3. SOLIDS POSSIBLE : NAHSO4
6390 ! *** COPYRIGHT 1996-2006, UNIVERSITY OF MIAMI, CARNEGIE MELLON UNIVERSITY,
6391 ! *** GEORGIA INSTITUTE OF TECHNOLOGY
6392 ! *** WRITTEN BY ATHANASIOS NENES
6393 ! *** UPDATED BY CHRISTOS FOUNTOUKIS
6395 !=======================================================================
6397       SUBROUTINE CALCJ2
6398       USE ISRPIA
6399       USE CASEJ
6400       IMPLICIT DOUBLE PRECISION (A-H, O-Z)
6401 !      INCLUDE 'ISRPIA.INC'
6403 !      DOUBLE PRECISION LAMDA, KAPA
6404 !      COMMON /CASEJ/ CHI1, CHI2, CHI3, LAMDA, KAPA, PSI1, PSI2, PSI3,   &
6405 !                     A1,   A2,   A3
6407 ! *** SETUP PARAMETERS ************************************************
6409       CALAOU = .TRUE.              ! Outer loop activity calculation flag
6410       CHI1   = W(1)                ! NA TOTAL
6411       CHI2   = W(3)                ! NH4 TOTAL
6412       PSI1LO = TINY                ! Low  limit
6413       PSI1HI = CHI1                ! High limit
6415 ! *** INITIAL VALUES FOR BISECTION ************************************
6417       X1 = PSI1HI
6418       Y1 = FUNCJ2 (X1)
6419       YHI= Y1                      ! Save Y-value at HI position
6421 ! *** YHI < 0.0 THE SOLUTION IS ALWAYS UNDERSATURATED WITH NH42SO4 ****
6423       IF (ABS(Y1).LE.EPS .OR. YHI.LT.ZERO) GOTO 50
6425 ! *** ROOT TRACKING ; FOR THE RANGE OF HI AND LO **********************
6427       DX = (PSI1HI-PSI1LO)/REAL(NDIV)
6428       DO 10 I=1,NDIV
6429          X2 = X1-DX
6430          Y2 = FUNCJ2 (X2)
6431          IF (SIGN(1.d0,Y1)*SIGN(1.d0,Y2).LT.ZERO) GOTO 20  ! (Y1*Y2.LT.ZERO)
6432          X1 = X2
6433          Y1 = Y2
6434 10    CONTINUE
6436 ! *** { YLO, YHI } > 0.0 THE SOLUTION IS ALWAYS SUPERSATURATED WITH NH42SO4
6438       YLO= Y1                      ! Save Y-value at Hi position
6439       IF (YLO.GT.ZERO .AND. YHI.GT.ZERO) THEN
6440          Y3 = FUNCJ2 (ZERO)
6441          GOTO 50
6442       ELSE IF (ABS(Y2) .LT. EPS) THEN   ! X2 IS A SOLUTION
6443          GOTO 50
6444       ELSE
6445          CALL PUSHERR (0001, 'CALCJ2')    ! WARNING ERROR: NO SOLUTION
6446          GOTO 50
6447       ENDIF
6449 ! *** PERFORM BISECTION ***********************************************
6451 20    DO 30 I=1,MAXIT
6452          X3 = 0.5*(X1+X2)
6453          Y3 = FUNCJ2 (X3)
6454          IF (SIGN(1.d0,Y1)*SIGN(1.d0,Y3) .LE. ZERO) THEN  ! (Y1*Y3 .LE. ZERO)
6455             Y2    = Y3
6456             X2    = X3
6457          ELSE
6458             Y1    = Y3
6459             X1    = X3
6460          ENDIF
6461          IF (ABS(X2-X1) .LE. EPS*X1) GOTO 40
6462 30    CONTINUE
6463       CALL PUSHERR (0002, 'CALCJ2')    ! WARNING ERROR: NO CONVERGENCE
6465 ! *** CONVERGED ; RETURN **********************************************
6467 40    X3 = 0.5*(X1+X2)
6468       Y3 = FUNCJ2 (X3)
6470 50    RETURN
6472 ! *** END OF SUBROUTINE CALCJ2 ******************************************
6474       END SUBROUTINE CALCJ2
6479 !=======================================================================
6481 ! *** ISORROPIA CODE
6482 ! *** SUBROUTINE FUNCJ2
6483 ! *** CASE J2
6485 !     THE MAIN CHARACTERISTICS OF THIS REGIME ARE:
6486 !     1. SULFATE RICH, FREE ACID (SULRAT < 1.0)
6487 !     2. THERE IS BOTH A LIQUID & SOLID PHASE
6488 !     3. SOLIDS POSSIBLE : (NH4)2SO4, NH4CL, NA2SO4
6490 ! *** COPYRIGHT 1996-2006, UNIVERSITY OF MIAMI, CARNEGIE MELLON UNIVERSITY,
6491 ! *** GEORGIA INSTITUTE OF TECHNOLOGY
6492 ! *** WRITTEN BY ATHANASIOS NENES
6493 ! *** UPDATED BY CHRISTOS FOUNTOUKIS
6495 !=======================================================================
6497       DOUBLE PRECISION FUNCTION FUNCJ2 (P1)
6498       USE CASEJ
6499       USE ISRPIA
6500       IMPLICIT DOUBLE PRECISION (A-H, O-Z)
6501 !      INCLUDE 'ISRPIA.INC'
6503 !      DOUBLE PRECISION LAMDA, KAPA
6504 !      COMMON /CASEJ/ CHI1, CHI2, CHI3, LAMDA, KAPA, PSI1, PSI2, PSI3,   &
6505 !                     A1,   A2,   A3
6507 ! *** SETUP PARAMETERS ************************************************
6509       FRST   = .TRUE.
6510       CALAIN = .TRUE.
6512       LAMDA  = MAX(W(2) - W(3) - W(1), TINY)  ! FREE H2SO4
6513       PSI1   = P1
6514       PSI2   = CHI2                           ! ALL NH4HSO4 DELIQUESCED
6516 ! *** SOLVE EQUATIONS ; WITH ITERATIONS FOR ACTIVITY COEF. ************
6518       DO 10 I=1,NSWEEP
6520       A1 = XK11 *(WATER/GAMA(12))**2.0
6521       A3 = XK1  *WATER/GAMA(7)*(GAMA(8)/GAMA(7))**2.0
6523 !  CALCULATE DISSOCIATION QUANTITIES
6525       BB   = A3+LAMDA                        ! KAPA
6526       CC   =-A3*(LAMDA + PSI1 + PSI2)
6527       DD   = BB*BB-4.D0*CC
6528       KAPA = 0.5D0*(-BB+SQRT(DD))
6530 ! *** CALCULATE SPECIATION ********************************************
6532       MOLAL (1) = LAMDA + KAPA                  ! HI
6533       MOLAL (2) = PSI1                          ! NAI
6534       MOLAL (3) = PSI2                          ! NH4I
6535       MOLAL (4) = ZERO                          ! CLI
6536       MOLAL (5) = KAPA                          ! SO4I
6537       MOLAL (6) = LAMDA + PSI1 + PSI2 - KAPA    ! HSO4I
6538       MOLAL (7) = ZERO                          ! NO3I
6540       CNAHSO4   = MAX(CHI1-PSI1,ZERO)
6541       CNH4HS4   = ZERO
6543       CALL CALCMR                               ! Water content
6545 ! *** CALCULATE ACTIVITIES OR TERMINATE INTERNAL LOOP *****************
6547       IF (FRST.AND.CALAOU .OR. .NOT.FRST.AND.CALAIN) THEN
6548          CALL CALCACT
6549       ELSE
6550          GOTO 20
6551       ENDIF
6552 10    CONTINUE
6554 ! *** CALCULATE OBJECTIVE FUNCTION ************************************
6556 20    FUNCJ2 = MOLAL(2)*MOLAL(6)/A1 - ONE
6558 ! *** END OF FUNCTION FUNCJ2 *******************************************
6560       END FUNCTION FUNCJ2     
6562 !=======================================================================
6564 ! *** ISORROPIA CODE
6565 ! *** SUBROUTINE CALCJ1
6566 ! *** CASE J1
6568 !     THE MAIN CHARACTERISTICS OF THIS REGIME ARE:
6569 !     1. SULFATE RICH, FREE ACID (SULRAT < 1.0)
6570 !     2. THERE IS BOTH A LIQUID & SOLID PHASE
6571 !     3. SOLIDS POSSIBLE : NH4HSO4, NAHSO4
6573 ! *** COPYRIGHT 1996-2006, UNIVERSITY OF MIAMI, CARNEGIE MELLON UNIVERSITY,
6574 ! *** GEORGIA INSTITUTE OF TECHNOLOGY
6575 ! *** WRITTEN BY ATHANASIOS NENES
6576 ! *** UPDATED BY CHRISTOS FOUNTOUKIS
6578 !=======================================================================
6580       SUBROUTINE CALCJ1
6581       USE ISRPIA
6582       USE CASEJ
6583       IMPLICIT DOUBLE PRECISION (A-H, O-Z)
6584 !      INCLUDE 'ISRPIA.INC'
6586 !      DOUBLE PRECISION LAMDA, KAPA
6587 !      COMMON /CASEJ/ CHI1, CHI2, CHI3, LAMDA, KAPA, PSI1, PSI2, PSI3,   &
6588 !                     A1,   A2,   A3
6590 ! *** SETUP PARAMETERS ************************************************
6592       CALAOU =.TRUE.               ! Outer loop activity calculation flag
6593       CHI1   = W(1)                ! Total NA initially as NaHSO4
6594       CHI2   = W(3)                ! Total NH4 initially as NH4HSO4
6596       PSI1LO = TINY                ! Low  limit
6597       PSI1HI = CHI1                ! High limit
6599 ! *** INITIAL VALUES FOR BISECTION ************************************
6601       X1 = PSI1HI
6602       Y1 = FUNCJ1 (X1)
6603       YHI= Y1                      ! Save Y-value at HI position
6605 ! *** YHI < 0.0 THE SOLUTION IS ALWAYS UNDERSATURATED WITH NH42SO4 ****
6607       IF (ABS(Y1).LE.EPS .OR. YHI.LT.ZERO) GOTO 50
6609 ! *** ROOT TRACKING ; FOR THE RANGE OF HI AND LO **********************
6611       DX = (PSI1HI-PSI1LO)/REAL(NDIV)
6612       DO 10 I=1,NDIV
6613          X2 = X1-DX
6614          Y2 = FUNCJ1 (X2)
6615          IF (SIGN(1.d0,Y1)*SIGN(1.d0,Y2).LT.ZERO) GOTO 20  ! (Y1*Y2.LT.ZERO)
6616          X1 = X2
6617          Y1 = Y2
6618 10    CONTINUE
6620 ! *** { YLO, YHI } > 0.0 THE SOLUTION IS ALWAYS SUPERSATURATED WITH NH42SO4
6622       YLO= Y1                      ! Save Y-value at Hi position
6623       IF (YLO.GT.ZERO .AND. YHI.GT.ZERO) THEN
6624          Y3 = FUNCJ1 (ZERO)
6625          GOTO 50
6626       ELSE IF (ABS(Y2) .LT. EPS) THEN   ! X2 IS A SOLUTION
6627          GOTO 50
6628       ELSE
6629          CALL PUSHERR (0001, 'CALCJ1')    ! WARNING ERROR: NO SOLUTION
6630          GOTO 50
6631       ENDIF
6633 ! *** PERFORM BISECTION ***********************************************
6635 20    DO 30 I=1,MAXIT
6636          X3 = 0.5*(X1+X2)
6637          Y3 = FUNCJ1 (X3)
6638          IF (SIGN(1.d0,Y1)*SIGN(1.d0,Y3) .LE. ZERO) THEN  ! (Y1*Y3 .LE. ZERO)
6639             Y2    = Y3
6640             X2    = X3
6641          ELSE
6642             Y1    = Y3
6643             X1    = X3
6644          ENDIF
6645          IF (ABS(X2-X1) .LE. EPS*X1) GOTO 40
6646 30    CONTINUE
6647       CALL PUSHERR (0002, 'CALCJ1')    ! WARNING ERROR: NO CONVERGENCE
6649 ! *** CONVERGED ; RETURN **********************************************
6651 40    X3 = 0.5*(X1+X2)
6652       Y3 = FUNCJ1 (X3)
6654 50    RETURN
6656 ! *** END OF SUBROUTINE CALCJ1 ******************************************
6658       END SUBROUTINE CALCJ1
6663 !=======================================================================
6665 ! *** ISORROPIA CODE
6666 ! *** SUBROUTINE FUNCJ1
6667 ! *** CASE J1
6669 !     THE MAIN CHARACTERISTICS OF THIS REGIME ARE:
6670 !     1. SULFATE RICH, FREE ACID (SULRAT < 1.0)
6671 !     2. THERE IS BOTH A LIQUID & SOLID PHASE
6672 !     3. SOLIDS POSSIBLE : (NH4)2SO4, NH4CL, NA2SO4
6674 ! *** COPYRIGHT 1996-2006, UNIVERSITY OF MIAMI, CARNEGIE MELLON UNIVERSITY,
6675 ! *** GEORGIA INSTITUTE OF TECHNOLOGY
6676 ! *** WRITTEN BY ATHANASIOS NENES
6677 ! *** UPDATED BY CHRISTOS FOUNTOUKIS
6679 !=======================================================================
6681       DOUBLE PRECISION FUNCTION FUNCJ1 (P1)
6682       USE ISRPIA
6683       USE CASEJ
6684       IMPLICIT DOUBLE PRECISION (A-H, O-Z)
6685 !      INCLUDE 'ISRPIA.INC'
6686 !      DOUBLE PRECISION LAMDA, KAPA
6687 !      COMMON /CASEJ/ CHI1, CHI2, CHI3, LAMDA, KAPA, PSI1, PSI2, PSI3,   &
6688 !                     A1,   A2,   A3
6690 ! *** SETUP PARAMETERS ************************************************
6692       FRST   = .TRUE.
6693       CALAIN = .TRUE.
6695       LAMDA  = MAX(W(2) - W(3) - W(1), TINY)  ! FREE H2SO4
6696       PSI1   = P1
6698 ! *** SOLVE EQUATIONS ; WITH ITERATIONS FOR ACTIVITY COEF. ************
6700       DO 10 I=1,NSWEEP
6702       A1 = XK11 *(WATER/GAMA(12))**2.0
6703       A2 = XK12 *(WATER/GAMA(09))**2.0
6704       A3 = XK1  *WATER/GAMA(7)*(GAMA(8)/GAMA(7))**2.0
6706       PSI2 = 0.5*(-(LAMDA+PSI1) + SQRT((LAMDA+PSI1)**2.D0+4.D0*A2))  ! PSI2
6707       PSI2 = MIN (PSI2, CHI2)
6709       BB   = A3+LAMDA                        ! KAPA
6710       CC   =-A3*(LAMDA + PSI2 + PSI1)
6711       DD   = BB*BB-4.D0*CC
6712       KAPA = 0.5D0*(-BB+SQRT(DD))
6714 ! *** SAVE CONCENTRATIONS IN MOLAL ARRAY ******************************
6716       MOLAL (1) = LAMDA + KAPA                  ! HI
6717       MOLAL (2) = PSI1                          ! NAI
6718       MOLAL (3) = PSI2                          ! NH4I
6719       MOLAL (4) = ZERO
6720       MOLAL (5) = KAPA                          ! SO4I
6721       MOLAL (6) = LAMDA + PSI1 + PSI2 - KAPA    ! HSO4I
6722       MOLAL (7) = ZERO
6724       CNAHSO4   = MAX(CHI1-PSI1,ZERO)
6725       CNH4HS4   = MAX(CHI2-PSI2,ZERO)
6727       CALL CALCMR                               ! Water content
6729 ! *** CALCULATE ACTIVITIES OR TERMINATE INTERNAL LOOP *****************
6731       IF (FRST.AND.CALAOU .OR. .NOT.FRST.AND.CALAIN) THEN
6732          CALL CALCACT
6733       ELSE
6734          GOTO 20
6735       ENDIF
6736 10    CONTINUE
6738 ! *** CALCULATE OBJECTIVE FUNCTION ************************************
6740 20    FUNCJ1 = MOLAL(2)*MOLAL(6)/A1 - ONE
6742 ! *** END OF FUNCTION FUNCJ1 *******************************************
6744       END FUNCTION FUNCJ1