1 !=======================================================================
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
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)
24 IMPLICIT DOUBLE PRECISION (A-H, O-Z)
25 ! INCLUDE 'ISRPIA.INC'
28 ! *** INITIALIZE ALL VARIABLES IN COMMON BLOCK **************************
30 CALL INIT1 (WI, RHI, TEMPI)
32 ! *** CALCULATE SULFATE RATIO *******************************************
36 ! *** FIND CALCULATION REGIME FROM (SULRAT,RH) **************************
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)
46 CALL CALCA2 ! Only liquid (metastable)
49 IF (RH.LT.DRNH42S4) THEN
51 CALL CALCA1 ! NH42SO4 ; case A1
53 ELSEIF (DRNH42S4.LE.RH) THEN
55 CALL CALCA2 ! Only liquid ; case A2
59 ! *** SULFATE RICH (NO ACID)
61 ELSEIF (1.0.LE.SULRAT .AND. SULRAT.LT.2.0) THEN
65 CALL CALCB4 ! Only liquid (metastable)
68 IF (RH.LT.DRNH4HS4) THEN
70 CALL CALCB1 ! NH4HSO4,LC,NH42SO4 ; case B1
72 ELSEIF (DRNH4HS4.LE.RH .AND. RH.LT.DRLC) THEN
74 CALL CALCB2 ! LC,NH42S4 ; case B2
76 ELSEIF (DRLC.LE.RH .AND. RH.LT.DRNH42S4) THEN
78 CALL CALCB3 ! NH42S4 ; case B3
80 ELSEIF (DRNH42S4.LE.RH) THEN
82 CALL CALCB4 ! Only liquid ; case B4
87 ! *** SULFATE RICH (FREE ACID)
89 ELSEIF (SULRAT.LT.1.0) THEN
93 CALL CALCC2 ! Only liquid (metastable)
96 IF (RH.LT.DRNH4HS4) THEN
98 CALL CALCC1 ! NH4HSO4 ; case C1
100 ELSEIF (DRNH4HS4.LE.RH) THEN
102 CALL CALCC2 ! Only liquid ; case C2
113 ! *** END OF SUBROUTINE ISRP1F *****************************************
115 END SUBROUTINE ISRP1F
116 !=======================================================================
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)
134 IMPLICIT DOUBLE PRECISION (A-H, O-Z)
135 ! INCLUDE 'ISRPIA.INC'
138 ! *** INITIALIZE ALL VARIABLES IN COMMON BLOCK **************************
140 CALL INIT2 (WI, RHI, TEMPI)
142 ! *** CALCULATE SULFATE RATIO *******************************************
146 ! *** FIND CALCULATION REGIME FROM (SULRAT,RH) **************************
150 IF (2.0.LE.SULRAT) THEN
152 IF(METSTBL.EQ.1) THEN
154 CALL CALCD3 ! Only liquid (metastable)
157 IF (RH.LT.DRNH4NO3) THEN
159 CALL CALCD1 ! NH42SO4,NH4NO3 ; case D1
161 ELSEIF (DRNH4NO3.LE.RH .AND. RH.LT.DRNH42S4) THEN
163 CALL CALCD2 ! NH42S4 ; case D2
165 ELSEIF (DRNH42S4.LE.RH) THEN
167 CALL CALCD3 ! Only liquid ; case D3
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
181 CALL CALCB4 ! Only liquid (metastable)
185 IF (RH.LT.DRNH4HS4) THEN
187 CALL CALCB1 ! NH4HSO4,LC,NH42SO4 ; case E1
190 ELSEIF (DRNH4HS4.LE.RH .AND. RH.LT.DRLC) THEN
192 CALL CALCB2 ! LC,NH42S4 ; case E2
195 ELSEIF (DRLC.LE.RH .AND. RH.LT.DRNH42S4) THEN
197 CALL CALCB3 ! NH42S4 ; case E3
200 ELSEIF (DRNH42S4.LE.RH) THEN
202 CALL CALCB4 ! Only liquid ; case E4
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
219 CALL CALCC2 ! Only liquid (metastable)
223 IF (RH.LT.DRNH4HS4) THEN
225 CALL CALCC1 ! NH4HSO4 ; case F1
228 ELSEIF (DRNH4HS4.LE.RH) THEN
230 CALL CALCC2 ! Only liquid ; case F2
235 CALL CALCNA ! HNO3(g) DISSOLUTION
242 ! *** END OF SUBROUTINE ISRP2F *****************************************
244 END SUBROUTINE ISRP2F
245 !=======================================================================
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)
263 IMPLICIT DOUBLE PRECISION (A-H, O-Z)
264 ! INCLUDE 'ISRPIA.INC'
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
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
291 ! *** CALCULATE SULFATE & SODIUM RATIOS *********************************
293 SULRAT = (W(1)+W(3))/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
304 CALL CALCG5 ! Only liquid (metastable)
307 IF (RH.LT.DRNH4NO3) THEN
309 CALL CALCG1 ! NH42SO4,NH4NO3,NH4CL,NA2SO4
311 ELSEIF (DRNH4NO3.LE.RH .AND. RH.LT.DRNH4CL) THEN
313 CALL CALCG2 ! NH42SO4,NH4CL,NA2SO4
315 ELSEIF (DRNH4CL.LE.RH .AND. RH.LT.DRNH42S4) THEN
317 CALL CALCG3 ! NH42SO4,NA2SO4
319 ELSEIF (DRNH42S4.LE.RH .AND. RH.LT.DRNA2SO4) THEN
323 ELSEIF (DRNA2SO4.LE.RH) THEN
325 CALL CALCG5 ! Only liquid
329 ! *** SULFATE POOR ; SODIUM RICH
331 ELSE IF (SULRAT.GE.2.0 .AND. SODRAT.GE.2.0) THEN
333 IF(METSTBL.EQ.1) THEN
335 CALL CALCH6 ! Only liquid (metastable)
338 IF (RH.LT.DRNH4NO3) THEN
340 CALL CALCH1 ! NH4NO3,NH4CL,NA2SO4,NACL,NANO3
342 ELSEIF (DRNH4NO3.LE.RH .AND. RH.LT.DRNANO3) THEN
344 CALL CALCH2 ! NH4CL,NA2SO4,NACL,NANO3
346 ELSEIF (DRNANO3.LE.RH .AND. RH.LT.DRNACL) THEN
348 CALL CALCH3 ! NH4CL,NA2SO4,NACL
350 ELSEIF (DRNACL.LE.RH .AND. RH.LT.DRNH4Cl) THEN
352 CALL CALCH4 ! NH4CL,NA2SO4
354 ELSEIF (DRNH4Cl.LE.RH .AND. RH.LT.DRNA2SO4) THEN
358 ELSEIF (DRNA2SO4.LE.RH) THEN
360 CALL CALCH6 ! NO SOLID
364 ! *** SULFATE RICH (NO ACID)
366 ELSEIF (1.0.LE.SULRAT .AND. SULRAT.LT.2.0) THEN
368 IF(METSTBL.EQ.1) THEN
370 CALL CALCI6 ! Only liquid (metastable)
373 IF (RH.LT.DRNH4HS4) THEN
375 CALL CALCI1 ! NA2SO4,(NH4)2SO4,NAHSO4,NH4HSO4,LC
377 ELSEIF (DRNH4HS4.LE.RH .AND. RH.LT.DRNAHSO4) THEN
379 CALL CALCI2 ! NA2SO4,(NH4)2SO4,NAHSO4,LC
381 ELSEIF (DRNAHSO4.LE.RH .AND. RH.LT.DRLC) THEN
383 CALL CALCI3 ! NA2SO4,(NH4)2SO4,LC
385 ELSEIF (DRLC.LE.RH .AND. RH.LT.DRNH42S4) THEN
387 CALL CALCI4 ! NA2SO4,(NH4)2SO4
389 ELSEIF (DRNH42S4.LE.RH .AND. RH.LT.DRNA2SO4) THEN
393 ELSEIF (DRNA2SO4.LE.RH) THEN
395 CALL CALCI6 ! NO SOLIDS
399 CALL CALCNHA ! MINOR SPECIES: HNO3, HCl
402 ! *** SULFATE RICH (FREE ACID)
404 ELSEIF (SULRAT.LT.1.0) THEN
406 IF(METSTBL.EQ.1) THEN
408 CALL CALCJ3 ! Only liquid (metastable)
411 IF (RH.LT.DRNH4HS4) THEN
413 CALL CALCJ1 ! NH4HSO4,NAHSO4
415 ELSEIF (DRNH4HS4.LE.RH .AND. RH.LT.DRNAHSO4) THEN
419 ELSEIF (DRNAHSO4.LE.RH) THEN
425 CALL CALCNHA ! MINOR SPECIES: HNO3, HCl
433 ! *** END OF SUBROUTINE ISRP3F *****************************************
435 END SUBROUTINE ISRP3F
436 !=======================================================================
439 ! *** SUBROUTINE CALCA2
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 !=======================================================================
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 *****************************************
476 ! *** INITIAL VALUES FOR BISECTION ************************************
480 IF (ABS(Y1).LE.EPS) RETURN
482 ! *** ROOT TRACKING ; FOR THE RANGE OF HI AND LO **********************
484 DX = (OMEHI-OMELO)/REAL(NDIV)
486 X2 = MAX(X1-DX, OMELO)
488 IF (SIGN(1.d0,Y1)*SIGN(1.d0,Y2).LT.ZERO) GOTO 20 ! (Y1*Y2.LT.ZERO)
492 IF (ABS(Y2).LE.EPS) THEN
495 CALL PUSHERR (0001, 'CALCA2') ! WARNING ERROR: NO SOLUTION
499 ! *** PERFORM BISECTION ***********************************************
504 IF (SIGN(1.d0,Y1)*SIGN(1.d0,Y3) .LE. ZERO) THEN ! (Y1*Y3 .LE. ZERO)
511 IF (ABS(X2-X1) .LE. EPS*X1) GOTO 40
513 CALL PUSHERR (0002, 'CALCA2') ! WARNING ERROR: NO CONVERGENCE
515 ! *** CONVERGED ; RETURN **********************************************
521 ! *** END OF SUBROUTINE CALCA2 ****************************************
523 END SUBROUTINE CALCA2
527 !=======================================================================
530 ! *** FUNCTION FUNCA2
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)
539 IMPLICIT DOUBLE PRECISION (A-H, O-Z)
540 ! INCLUDE 'ISRPIA.INC'
541 DOUBLE PRECISION LAMDA
543 ! *** SETUP PARAMETERS ************************************************
547 PSI = W(2) ! INITIAL AMOUNT OF (NH4)2SO4 IN SOLUTION
549 ! *** SOLVE EQUATIONS ; WITH ITERATIONS FOR ACTIVITY COEF. ************
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)
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
568 ! *** CALCULATE ACTIVITIES OR TERMINATE INTERNAL LOOP *****************
570 IF (FRST.AND.CALAOU .OR. .NOT.FRST.AND.CALAIN) THEN
577 ! *** CALCULATE OBJECTIVE FUNCTION ************************************
579 20 DENOM = (2.0*MOLAL(5)+MOLAL(6))
580 FUNCA2= (MOLAL(3)/DENOM - ONE) + MOLAL(1)/DENOM
583 ! *** END OF FUNCTION FUNCA2 ********************************************
586 !=======================================================================
589 ! *** SUBROUTINE CALCA1
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
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 !=======================================================================
610 IMPLICIT DOUBLE PRECISION (A-H, O-Z)
611 ! INCLUDE 'ISRPIA.INC'
614 GNH3 = MAX (W(3)-2.0*CNH42S4, ZERO)
617 ! *** END OF SUBROUTINE CALCA1 ******************************************
619 END SUBROUTINE CALCA1
623 !=======================================================================
626 ! *** SUBROUTINE CALCB4
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 !=======================================================================
646 IMPLICIT DOUBLE PRECISION (A-H, O-Z)
647 ! INCLUDE 'ISRPIA.INC'
649 ! *** SOLVE EQUATIONS **************************************************
655 ! *** CALCULATE WATER CONTENT ******************************************
657 CALL CALCB1A ! GET DRY SALT CONTENT, AND USE FOR WATER.
664 WATER = MOLALR(13)/M0(13)+MOLALR(9)/M0(9)+MOLALR(4)/M0(4)
666 MOLAL(3) = W(3) ! NH4I
669 AK1 = XK1*((GAMA(8)/GAMA(7))**2.)*(WATER/GAMA(7))
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
692 ! *** END OF SUBROUTINE CALCB4 ******************************************
694 END SUBROUTINE CALCB4
695 !=======================================================================
698 ! *** SUBROUTINE CALCB3
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 !=======================================================================
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'
729 CALL CALCB3A (TLC,TNH42S4) ! LC + (NH4)2SO4
731 SCASE = 'B3 ; SUBCASE 2'
734 CALL CALCB3B (TLC,TNH4HS4) ! LC + NH4HSO4
739 ! *** END OF SUBROUTINE CALCB3 ******************************************
741 END SUBROUTINE CALCB3
744 !=======================================================================
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)
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 ****************
781 Y1 = FUNCB3A (Z1, TLC, TNH42S4)
782 IF (ABS(Y1).LE.EPS) RETURN
785 ! *** ROOT TRACKING ; FOR THE RANGE OF HI AND LO ***********************
787 DZ = (ZHI-ZLO)/REAL(NDIV)
790 Y2 = FUNCB3A (Z2, TLC, TNH42S4)
791 IF (SIGN(1.d0,Y1)*SIGN(1.d0,Y2).LT.ZERO) GOTO 20 ! (Y1*Y2.LT.ZERO)
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
802 ! *** { YLO, YHI } < 0.0 THE SOLUTION IS ALWAYS UNDERSATURATED WITH LC
804 ELSE IF (YLO.LT.ZERO .AND. YHI.LT.ZERO) THEN
809 ! *** { YLO, YHI } > 0.0 THE SOLUTION IS ALWAYS SUPERSATURATED WITH LC
811 ELSE IF (YLO.GT.ZERO .AND. YHI.GT.ZERO) THEN
816 CALL PUSHERR (0001, 'CALCB3A') ! WARNING ERROR: NO SOLUTION
820 ! *** PERFORM BISECTION ***********************************************
824 Y3 = FUNCB3A (Z3, TLC, TNH42S4)
825 IF (SIGN(1.d0,Y1)*SIGN(1.d0,Y3) .LE. ZERO) THEN ! (Y1*Y3 .LE. ZERO)
832 IF (ABS(Z2-Z1) .LE. EPS*Z1) GOTO 40
834 CALL PUSHERR (0002, 'CALCB3A') ! WARNING ERROR: NO CONVERGENCE
836 ! *** CONVERGED ; RETURN ************************************************
839 Y3 = FUNCB3A (ZK, TLC, TNH42S4)
843 ! *** END OF SUBROUTINE CALCB3A ******************************************
845 END SUBROUTINE CALCB3A
849 !=======================================================================
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)
861 IMPLICIT DOUBLE PRECISION (A-H, O-Z)
862 ! INCLUDE 'ISRPIA.INC'
865 ! *** SOLVE EQUATIONS ; WITH ITERATIONS FOR ACTIVITY COEF. ************
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 ***************************************
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
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
899 ! *** END OF FUNCTION FUNCB3A ********************************************
905 !=======================================================================
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)
926 IMPLICIT DOUBLE PRECISION (A-H, O-Z)
927 ! INCLUDE 'ISRPIA.INC'
930 CALAOU = .FALSE. ! Outer loop activity calculation flag
934 ! *** SOLVE EQUATIONS ; WITH ITERATIONS FOR ACTIVITY COEF. ************
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 ***************************************
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
957 ! *** END OF SUBROUTINE CALCB3B ******************************************
959 END SUBROUTINE CALCB3B
960 !=======================================================================
963 ! *** SUBROUTINE CALCB2
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 !=======================================================================
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
998 SCASE = 'B2 ; SUBCASE 2'
999 CALL CALCB2B (Y,X-Y) ! LC ONLY POSSIBLE
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)
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
1047 SCASE = 'B2 ; SUBCASE A1'
1049 SCASE = 'B2 ; SUBCASE A2'
1050 CALL CALCB2A2 (TLC, TNH42S4) ! LIQUID & SOLID PHASE POSSIBLE
1051 SCASE = 'B2 ; SUBCASE A2'
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)
1087 IMPLICIT DOUBLE PRECISION (A-H, O-Z)
1088 ! INCLUDE 'ISRPIA.INC'
1090 ! *** FIND WEIGHT FACTOR **********************************************
1092 IF (WFTYP.EQ.0) THEN
1094 ELSEIF (WFTYP.EQ.1) THEN
1097 WF = (DRLC-RH)/(DRLC-DRMLCAS)
1101 ! *** FIND FIRST SECTION ; DRY ONE ************************************
1103 CLCO = TLC ! FIRST (DRY) SOLUTION
1106 ! *** FIND SECOND SECTION ; DRY & LIQUID ******************************
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
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
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)
1159 IMPLICIT DOUBLE PRECISION (A-H, O-Z)
1160 ! INCLUDE 'ISRPIA.INC'
1162 CALAOU = .TRUE. ! Outer loop activity calculation flag
1164 ZHI = TLC ! High limit: all of it in liquid phase
1166 ! *** INITIAL VALUES FOR BISECTION **************************************
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 ************************
1178 Y2 = FUNCB2B (X2,TNH4HS4,TLC)
1179 IF (SIGN(1.d0,Y1)*SIGN(1.d0,Y2).LT.ZERO) GOTO 20 ! (Y1*Y2.LT.ZERO)
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
1190 ! *** { YLO, YHI } < 0.0 THE SOLUTION IS ALWAYS UNDERSATURATED WITH LC
1192 ELSE IF (YLO.LT.ZERO .AND. YHI.LT.ZERO) THEN
1197 ! *** { YLO, YHI } > 0.0 THE SOLUTION IS ALWAYS SUPERSATURATED WITH LC
1199 ELSE IF (YLO.GT.ZERO .AND. YHI.GT.ZERO) THEN
1204 CALL PUSHERR (0001, 'CALCB2B') ! WARNING ERROR: NO SOLUTION
1208 ! *** PERFORM BISECTION *************************************************
1212 Y3 = FUNCB2B (X3,TNH4HS4,TLC)
1213 IF (SIGN(1.d0,Y1)*SIGN(1.d0,Y3) .LE. ZERO) THEN ! (Y1*Y3 .LE. ZERO)
1220 IF (ABS(X2-X1) .LE. EPS*X1) GOTO 40
1222 CALL PUSHERR (0002, 'CALCB2B') ! WARNING ERROR: NO CONVERGENCE
1224 ! *** CONVERGED ; RETURN ************************************************
1227 Y3 = FUNCB2B (X3,TNH4HS4,TLC)
1231 ! *** END OF SUBROUTINE CALCB2B *****************************************
1233 END SUBROUTINE CALCB2B
1237 !=======================================================================
1239 ! *** ISORROPIA CODE
1240 ! *** FUNCTION FUNCB2B
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)
1249 IMPLICIT DOUBLE PRECISION (A-H, O-Z)
1250 ! INCLUDE 'ISRPIA.INC'
1252 ! *** SOLVE EQUATIONS **************************************************
1257 GRAT2 = XK1*WATER*(GAMA(8)/GAMA(7))**2./GAMA(7)
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
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
1287 ! *** END OF FUNCTION FUNCB2B *******************************************
1289 END FUNCTION FUNCB2B
1292 !=======================================================================
1294 ! *** ISORROPIA CODE
1295 ! *** SUBROUTINE CALCB1
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 !=======================================================================
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'
1326 SCASE = 'B1 ; SUBCASE 2'
1327 CALL CALCB1B ! LIQUID & SOLID PHASE POSSIBLE
1328 SCASE = 'B1 ; SUBCASE 2'
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:
1347 ! 2. THERE IS NO LIQUID PHASE
1348 ! 3. SOLIDS POSSIBLE: LC, { (NH4)2SO4 XOR NH4HSO4 } (ONE OF TWO
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 !=======================================================================
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
1380 CLC = Y ! NH4HSO4 < (NH4)2S04
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:
1401 ! 2. THERE IS BOTH A LIQUID & SOLID PHASE
1402 ! 3. SOLIDS POSSIBLE: LC, { (NH4)2SO4 XOR NH4HSO4 } (ONE OF TWO
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 !=======================================================================
1419 IMPLICIT DOUBLE PRECISION (A-H, O-Z)
1420 ! INCLUDE 'ISRPIA.INC'
1422 ! *** FIND WEIGHT FACTOR **********************************************
1424 IF (WFTYP.EQ.0) THEN
1426 ELSEIF (WFTYP.EQ.1) THEN
1429 WF = (DRNH4HS4-RH)/(DRNH4HS4-DRMLCAB)
1433 ! *** FIND FIRST SECTION ; DRY ONE ************************************
1436 CLCO = CLC ! FIRST (DRY) SOLUTION
1440 ! *** FIND SECOND SECTION ; DRY & LIQUID ******************************
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
1463 ! *** END OF SUBROUTINE CALCB1B *****************************************
1465 END SUBROUTINE CALCB1B
1468 !=======================================================================
1470 ! *** ISORROPIA CODE
1471 ! *** SUBROUTINE CALCC2
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 !=======================================================================
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
1495 ! *** SOLVE EQUATIONS **************************************************
1497 LAMDA = W(3) ! NH4HSO4 INITIALLY IN SOLUTION
1498 PSI = W(2)-W(3) ! H2SO4 IN SOLUTION
1500 PARM = WATER*XK1/GAMA(7)*(GAMA(8)/GAMA(7))**2.
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
1522 ! *** END OF SUBROUTINE CALCC2 *****************************************
1524 END SUBROUTINE CALCC2
1528 !=======================================================================
1530 ! *** ISORROPIA CODE
1531 ! *** SUBROUTINE CALCC1
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 !=======================================================================
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
1556 ! *** INITIAL VALUES FOR BISECTION *************************************
1560 IF (ABS(Y1).LE.EPS) GOTO 50
1563 ! *** ROOT TRACKING ; FOR THE RANGE OF HI AND LO ***********************
1565 DX = (KHI-KLO)/REAL(NDIV)
1569 IF (SIGN(1.d0,Y1)*SIGN(1.d0,Y2) .LT. ZERO) GOTO 20 ! (Y1*Y2 .LT. ZERO)
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
1580 ! *** { YLO, YHI } < 0.0 SOLUTION IS ALWAYS UNDERSATURATED WITH NH4HS04
1582 ELSE IF (YLO.LT.ZERO .AND. YHI.LT.ZERO) THEN
1585 ! *** { YLO, YHI } > 0.0 SOLUTION IS ALWAYS SUPERSATURATED WITH NH4HS04
1587 ELSE IF (YLO.GT.ZERO .AND. YHI.GT.ZERO) THEN
1592 CALL PUSHERR (0001, 'CALCC1') ! WARNING ERROR: NO SOLUTION
1596 ! *** PERFORM BISECTION OF DISSOLVED NH4HSO4 **************************
1601 IF (SIGN(1.d0,Y1)*SIGN(1.d0,Y3) .LE. ZERO) THEN ! (Y1*Y3 .LE. ZERO)
1608 IF (ABS(X2-X1) .LE. EPS*X1) GOTO 40
1610 CALL PUSHERR (0002, 'CALCC1') ! WARNING ERROR: NO CONVERGENCE
1612 ! *** CONVERGED ; RETURN ***********************************************
1619 ! *** END OF SUBROUTINE CALCC1 *****************************************
1621 END SUBROUTINE CALCC1
1625 !=======================================================================
1627 ! *** ISORROPIA CODE
1628 ! *** FUNCTION FUNCC1
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)
1642 IMPLICIT DOUBLE PRECISION (A-H, O-Z)
1643 ! INCLUDE 'ISRPIA.INC'
1644 DOUBLE PRECISION KAPA, LAMDA
1646 ! *** SOLVE EQUATIONS **************************************************
1653 PAR1 = XK1*WATER/GAMA(7)*(GAMA(8)/GAMA(7))**2.0
1654 PAR2 = XK12*(WATER/GAMA(9))**2.0
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
1678 ! *** CALCULATE ZERO FUNCTION *******************************************
1680 !CC30 FUNCC1= (NH4I*HSO4I/PAR2) - ONE
1681 30 FUNCC1= (MOLAL(3)*MOLAL(6)/PAR2) - ONE
1684 ! *** END OF FUNCTION FUNCC1 ********************************************
1688 !=======================================================================
1690 ! *** ISORROPIA CODE
1691 ! *** SUBROUTINE CALCD3
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 !=======================================================================
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 **********************************************
1718 ! *** SETUP PARAMETERS ************************************************
1720 CHI1 = CNH4NO3 ! Save from CALCD1 run
1725 PSI1 = CNH4NO3 ! ASSIGN INITIAL PSI's
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 ************************************
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)
1753 IF (SIGN(1.d0,Y1)*SIGN(1.d0,Y2).LT.ZERO) GOTO 20 ! (Y1*Y2.LT.ZERO)
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
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
1769 ELSE IF (YLO.LT.ZERO .AND. YHI.LT.ZERO) THEN
1770 P4 = TINY ! PSI4LO ! CHI4
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
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
1790 CALL CALCMR ! Initial water
1791 GOTO 60 ! Redo root tracking
1795 ! *** PERFORM BISECTION ***********************************************
1800 IF (SIGN(1.d0,Y1)*SIGN(1.d0,Y3) .LE. ZERO) THEN ! (Y1*Y3 .LE. ZERO)
1807 IF (ABS(X2-X1) .LE. EPS*ABS(X1)) GOTO 40
1809 CALL PUSHERR (0002, 'CALCD3') ! WARNING ERROR: NO CONVERGENCE
1811 ! *** CONVERGED ; RETURN **********************************************
1816 ! *** CALCULATE HSO4 SPECIATION AND RETURN *******************************
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
1827 ! *** END OF SUBROUTINE CALCD3 ******************************************
1829 END SUBROUTINE CALCD3
1833 !=======================================================================
1835 ! *** ISORROPIA CODE
1836 ! *** FUNCTION FUNCD3
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)
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 ************************************************
1858 ! *** SOLVE EQUATIONS ; WITH ITERATIONS FOR ACTIVITY COEF. ************
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)
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
1876 DENM = (BB+ABB) + 2.0*A7/ABB ! Taylor expansion of SQRT
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
1902 ! *** CALCULATE OBJECTIVE FUNCTION ************************************
1905 !CC FUNCD3= NH4I/HI/MAX(GNH3,TINY)/A4 - ONE
1906 FUNCD3= MOLAL(3)/MOLAL(1)/MAX(GNH3,TINY)/A4 - ONE
1909 ! *** END OF FUNCTION FUNCD3 ********************************************
1912 !=======================================================================
1914 ! *** ISORROPIA CODE
1915 ! *** SUBROUTINE CALCD2
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 !=======================================================================
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 **********************************************
1943 ! *** SETUP PARAMETERS ************************************************
1945 CHI1 = CNH4NO3 ! Save from CALCD1 run
1950 PSI1 = CNH4NO3 ! ASSIGN INITIAL PSI's
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 ************************************
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)
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)
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
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
1999 ELSE IF (YLO.LT.ZERO .AND. YHI.LT.ZERO) THEN
2000 P4 = TINY ! PSI4LO ! CHI4
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
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
2020 CALL CALCMR ! Initial water
2021 GOTO 60 ! Redo root tracking
2025 ! *** PERFORM BISECTION ***********************************************
2030 IF (SIGN(1.d0,Y1)*SIGN(1.d0,Y3) .LE. ZERO) THEN ! (Y1*Y3 .LE. ZERO)
2037 IF (ABS(X2-X1) .LE. EPS*ABS(X1)) GOTO 40
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.
2046 ! *** CALCULATE HSO4 SPECIATION AND RETURN *******************************
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
2057 ! *** END OF SUBROUTINE CALCD2 ******************************************
2059 END SUBROUTINE CALCD2
2063 !=======================================================================
2065 ! *** ISORROPIA CODE
2066 ! *** FUNCTION FUNCD2
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)
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
2090 ! *** SOLVE EQUATIONS ; WITH ITERATIONS FOR ACTIVITY COEF. ************
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
2100 CALL POLY3 (PSI14,0.25*PSI14**2.,-A2/4.D0, PSI2, ISLV) ! PSI2
2102 PSI2 = MIN (PSI2, CHI2)
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
2119 DENM = (BB+ABB) + 2.d0*A7/ABB ! Taylor expansion of SQRT
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
2145 ! *** CALCULATE OBJECTIVE FUNCTION ************************************
2148 !CC FUNCD2= NH4I/HI/MAX(GNH3,TINY)/A4 - ONE
2149 FUNCD2= MOLAL(3)/MOLAL(1)/MAX(GNH3,TINY)/A4 - ONE
2152 ! *** END OF FUNCTION FUNCD2 ********************************************
2155 !=======================================================================
2157 ! *** ISORROPIA CODE
2158 ! *** SUBROUTINE CALCD1
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 !=======================================================================
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
2188 SCASE = 'D1 ; SUBCASE 1'
2190 SCASE = 'D1 ; SUBCASE 2' ! LIQUID & SOLID PHASE POSSIBLE
2191 CALL CALCMDRH (RH, DRMASAN, DRNH4NO3, CALCD1A, CALCD2)
2192 SCASE = 'D1 ; SUBCASE 2'
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
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 !=======================================================================
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 *********************************
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)
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
2253 ! *** END OF SUBROUTINE CALCD1A *****************************************
2255 END SUBROUTINE CALCD1A
2256 !=======================================================================
2258 ! *** ISORROPIA CODE
2259 ! *** SUBROUTINE CALCG5
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 !=======================================================================
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 ************************************************
2289 CHI2 = MAX (W(2)-CHI1, ZERO)
2291 CHI4 = MAX (W(3)-2.D0*CHI2, ZERO)
2298 PSI6HI = CHI6-TINY ! MIN(CHI6-TINY, CHI4)
2300 WATER = CHI2/M0(4) + CHI1/M0(2)
2302 ! *** INITIAL VALUES FOR BISECTION ************************************
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)
2316 IF (SIGN(1.d0,Y1)*SIGN(1.d0,Y2).LT.ZERO) GOTO 20 ! (Y1*Y2.LT.ZERO)
2321 ! *** NO SUBDIVISION WITH SOLUTION; IF ABS(Y2)<EPS SOLUTION IS ASSUMED
2323 IF (ABS(Y2) .GT. EPS) Y2 = FUNCG5A (PSI6LO)
2326 ! *** PERFORM BISECTION ***********************************************
2331 IF (SIGN(1.d0,Y1)*SIGN(1.d0,Y3) .LE. ZERO) THEN ! (Y1*Y3 .LE. ZERO)
2338 IF (ABS(X2-X1) .LE. EPS*X1) GOTO 40
2340 CALL PUSHERR (0002, 'CALCG5') ! WARNING ERROR: NO CONVERGENCE
2342 ! *** CONVERGED ; RETURN **********************************************
2347 ! *** CALCULATE HSO4 SPECIATION AND RETURN *******************************
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
2359 ! *** END OF SUBROUTINE CALCG5 *******************************************
2361 END SUBROUTINE CALCG5
2366 !=======================================================================
2368 ! *** ISORROPIA CODE
2369 ! *** SUBROUTINE FUNCG5A
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)
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 ************************************************
2401 ! *** SOLVE EQUATIONS ; WITH ITERATIONS FOR ACTIVITY COEF. ************
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
2412 ! CALCULATE DISSOCIATION QUANTITIES
2414 IF (CHI5.GE.TINY) THEN
2415 PSI5 = PSI6*CHI5/(A6/A5*(CHI6-PSI6) + PSI6)
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))
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
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)
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
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
2469 ! *** END OF FUNCTION FUNCG5A *******************************************
2471 END FUNCTION FUNCG5A
2473 !=======================================================================
2475 ! *** ISORROPIA CODE
2476 ! *** SUBROUTINE CALCG4
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 !=======================================================================
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 ************************************************
2506 CHI2 = MAX (W(2)-CHI1, ZERO)
2508 CHI4 = MAX (W(3)-2.D0*CHI2, ZERO)
2514 PSI6HI = CHI6-TINY ! MIN(CHI6-TINY, CHI4)
2516 WATER = CHI2/M0(4) + CHI1/M0(2)
2518 ! *** INITIAL VALUES FOR BISECTION ************************************
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)
2532 IF (SIGN(1.d0,Y1)*SIGN(1.d0,Y2).LT.ZERO) GOTO 20 ! (Y1*Y2.LT.ZERO)
2537 ! *** NO SUBDIVISION WITH SOLUTION; IF ABS(Y2)<EPS SOLUTION IS ASSUMED
2539 IF (ABS(Y2) .GT. EPS) Y2 = FUNCG4A (PSI6LO)
2542 ! *** PERFORM BISECTION ***********************************************
2547 IF (SIGN(1.d0,Y1)*SIGN(1.d0,Y3) .LE. ZERO) THEN ! (Y1*Y3 .LE. ZERO)
2554 IF (ABS(X2-X1) .LE. EPS*X1) GOTO 40
2556 CALL PUSHERR (0002, 'CALCG4') ! WARNING ERROR: NO CONVERGENCE
2558 ! *** CONVERGED ; RETURN **********************************************
2563 ! *** CALCULATE HSO4 SPECIATION AND RETURN *******************************
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
2575 ! *** END OF SUBROUTINE CALCG4 *******************************************
2577 END SUBROUTINE CALCG4
2582 !=======================================================================
2584 ! *** ISORROPIA CODE
2585 ! *** SUBROUTINE FUNCG4A
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)
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 ************************************************
2619 ! *** SOLVE EQUATIONS ; WITH ITERATIONS FOR ACTIVITY COEF. ************
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)
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))
2647 ! CALCULATE CONCENTRATIONS
2649 NH4I = 2.0*PSI2 + PSI4
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)
2662 PSI1 = MIN (PSI1, CHI1)
2670 ! *** SAVE CONCENTRATIONS IN MOLAL ARRAY ******************************
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)
2689 CNA2SO4 = MAX(CHI1-PSI1,ZERO)
2691 ! *** CALCULATE MOLALR ARRAY, WATER AND ACTIVITIES **********************
2695 ! *** CALCULATE ACTIVITIES OR TERMINATE INTERNAL LOOP *****************
2697 IF (FRST.AND.CALAOU .OR. .NOT.FRST.AND.CALAIN) THEN
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
2711 ! *** END OF FUNCTION FUNCG4A *******************************************
2713 END FUNCTION FUNCG4A
2715 !=======================================================================
2717 ! *** ISORROPIA CODE
2718 ! *** SUBROUTINE CALCG3
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 !=======================================================================
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'
2744 SCASE = 'G3 ; SUBCASE 1'
2745 ELSE ! NO3, CL NON EXISTANT
2746 SCASE = 'G1 ; SUBCASE 1'
2748 SCASE = 'G1 ; SUBCASE 1'
2751 IF (WATER.LE.TINY) THEN
2752 IF (RH.LT.DRMG3) THEN ! ONLY SOLIDS
2758 SCASE = 'G3 ; SUBCASE 2'
2761 SCASE = 'G3 ; SUBCASE 3' ! MDRH REGION (NA2SO4, NH42S4)
2762 CALL CALCMDRH (RH, DRMG3, DRNH42S4, CALCG1A, CALCG4)
2763 SCASE = 'G3 ; SUBCASE 3'
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 !=======================================================================
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 ************************************************
2807 CHI2 = MAX (W(2)-CHI1, ZERO)
2809 CHI4 = MAX (W(3)-2.D0*CHI2, ZERO)
2814 PSI6HI = CHI6-TINY ! MIN(CHI6-TINY, CHI4)
2818 ! *** INITIAL VALUES FOR BISECTION ************************************
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)
2833 IF (SIGN(1.d0,Y1)*SIGN(1.d0,Y2).LT.ZERO) GOTO 20 ! (Y1*Y2.LT.ZERO)
2838 ! *** NO SUBDIVISION WITH SOLUTION; IF ABS(Y2)<EPS SOLUTION IS ASSUMED
2840 IF (ABS(Y2) .GT. EPS) Y2 = FUNCG3A (PSI6LO)
2843 ! *** PERFORM BISECTION ***********************************************
2848 IF (SIGN(1.d0,Y1)*SIGN(1.d0,Y3) .LE. ZERO) THEN ! (Y1*Y3 .LE. ZERO)
2855 IF (ABS(X2-X1) .LE. EPS*X1) GOTO 40
2857 CALL PUSHERR (0002, 'CALCG3A') ! WARNING ERROR: NO CONVERGENCE
2859 ! *** CONVERGED ; RETURN **********************************************
2864 ! *** FINAL CALCULATIONS *************************************************
2868 ! *** Na2SO4 DISSOLUTION
2870 IF (CHI1.GT.TINY .AND. WATER.GT.TINY) THEN ! PSI1
2871 CALL POLY3 (PSI2, ZERO, -A1/4.D0, PSI1, ISLV)
2873 PSI1 = MIN (PSI1, CHI1)
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
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)
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 ************************************************
2938 ! *** SOLVE EQUATIONS ; WITH ITERATIONS FOR ACTIVITY COEF. ************
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)
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))
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)
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)
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
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
3010 ! *** END OF FUNCTION FUNCG3A *******************************************
3012 END FUNCTION FUNCG3A
3014 !=======================================================================
3016 ! *** ISORROPIA CODE
3017 ! *** SUBROUTINE CALCG2
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 !=======================================================================
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'
3043 SCASE = 'G2 ; SUBCASE 1'
3044 ELSE ! NO3 NON EXISTANT, WATER NOT POSSIBLE
3045 SCASE = 'G1 ; SUBCASE 1'
3047 SCASE = 'G1 ; SUBCASE 1'
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
3059 SCASE = 'G2 ; SUBCASE 2'
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'
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'
3076 SCASE = 'G2 ; SUBCASE 2'
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 !=======================================================================
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 ************************************************
3121 CHI2 = MAX (W(2)-CHI1, ZERO)
3123 CHI4 = MAX (W(3)-2.D0*CHI2, ZERO)
3132 ! *** INITIAL VALUES FOR BISECTION ************************************
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)
3146 IF (SIGN(1.d0,Y1)*SIGN(1.d0,Y2).LT.ZERO) GOTO 20 ! (Y1*Y2.LT.ZERO)
3151 ! *** NO SUBDIVISION WITH SOLUTION; IF ABS(Y2)<EPS SOLUTION IS ASSUMED
3153 IF (ABS(Y2) .GT. EPS) WATER = TINY
3156 ! *** PERFORM BISECTION ***********************************************
3161 IF (SIGN(1.d0,Y1)*SIGN(1.d0,Y3) .LE. ZERO) THEN ! (Y1*Y3 .LE. ZERO)
3168 IF (ABS(X2-X1) .LE. EPS*X1) GOTO 40
3170 CALL PUSHERR (0002, 'CALCG2A') ! WARNING ERROR: NO CONVERGENCE
3172 ! *** CONVERGED ; RETURN **********************************************
3175 IF (X3.LE.TINY2) THEN ! PRACTICALLY NO NITRATES, SO DRY SOLUTION
3181 ! *** CALCULATE HSO4 SPECIATION AND RETURN *******************************
3185 ! *** Na2SO4 DISSOLUTION
3187 IF (CHI1.GT.TINY .AND. WATER.GT.TINY) THEN ! PSI1
3188 CALL POLY3 (PSI2, ZERO, -A1/4.D0, PSI1, ISLV)
3190 PSI1 = MIN (PSI1, CHI1)
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
3212 ! *** END OF SUBROUTINE CALCG2A ******************************************
3214 END SUBROUTINE CALCG2A
3219 !=======================================================================
3221 ! *** ISORROPIA CODE
3222 ! *** SUBROUTINE FUNCG2A
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)
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 ************************************************
3256 ! *** SOLVE EQUATIONS ; WITH ITERATIONS FOR ACTIVITY COEF. ************
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)
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)
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)
3299 ! *** NH4Cl(s) calculations
3301 A3 = XK6 /(R*TEMP*R*TEMP)
3302 IF (GNH3*GHCL.GT.A3) THEN
3303 DELT = MIN(GNH3, GHCL)
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
3311 ELSEIF (DELT-PSI32.GT.ZERO .AND. PSI32.GT.ZERO) THEN
3320 ! *** CALCULATE GAS / SOLID SPECIES (LIQUID IN MOLAL ALREADY) *********
3322 GNH3 = MAX(GNH3 - PSI3, TINY)
3323 GHCL = MAX(GHCL - PSI3, TINY)
3326 ! *** CALCULATE MOLALR ARRAY, WATER AND ACTIVITIES **********************
3330 ! *** CALCULATE ACTIVITIES OR TERMINATE INTERNAL LOOP *****************
3332 IF (FRST.AND.CALAOU .OR. .NOT.FRST.AND.CALAIN) THEN
3339 ! *** CALCULATE FUNCTION VALUE FOR OUTER LOOP ***************************
3341 20 IF (CHI4.LE.TINY) THEN
3342 FUNCG2A = MOLAL(1)*MOLAL(4)/GHCL/A6 - ONE
3344 FUNCG2A = MOLAL(3)*MOLAL(4)/GHCL/GNH3/A6/A4 - ONE
3349 ! *** END OF FUNCTION FUNCG2A *******************************************
3351 END FUNCTION FUNCG2A
3353 !=======================================================================
3355 ! *** ISORROPIA CODE
3356 ! *** SUBROUTINE CALCG1
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 !=======================================================================
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'
3388 SCASE = 'G1 ; SUBCASE 2' ! LIQUID & SOLID PHASE POSSIBLE
3389 CALL CALCMDRH (RH, DRMG1, DRNH4NO3, CALCG1A, CALCG2A)
3390 SCASE = 'G1 ; SUBCASE 2'
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
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 !=======================================================================
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 ***********************************
3432 CNH42S4 = W(2) - CNA2SO4
3434 ! *** CALCULATE VOLATILE SPECIES **************************************
3436 ALF = W(3) - 2.0*CNH42S4
3440 RTSQ = R*TEMP*R*TEMP
3444 THETA1 = GAM - BET*(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
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
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
3480 ! SEPERATE SOLUTION OF NH4CL & NH4NO3 EQUILIBRIA
3484 DD1 = (ALF+BET)*(ALF+BET) - 4.0D0*(ALF*BET-A1)
3485 DD2 = (ALF+GAM)*(ALF+GAM) - 4.0D0*(ALF*GAM-A2)
3489 IF (DD1.GE.ZERO) THEN
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
3496 ELSE IF (KAPA2.GE.ZERO .AND. KAPA2.LE.MIN(ALF,BET)) THEN
3503 ! NH4NO3 EQUILIBRIUM
3505 IF (DD2.GE.ZERO) THEN
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
3512 ELSE IF (LAMDA2.GE.ZERO .AND. LAMDA2.LE.MIN(ALF,GAM)) THEN
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
3529 ! *** CALCULATE COMPOSITION OF VOLATILE SPECIES ***********************
3535 GNH3 = MAX(ALF - KAPA - LAMDA, ZERO)
3536 GHNO3 = MAX(GAM - LAMDA, ZERO)
3537 GHCL = MAX(BET - KAPA, ZERO)
3541 ! *** END OF SUBROUTINE CALCG1A *****************************************
3543 END SUBROUTINE CALCG1A
3544 !=======================================================================
3546 ! *** ISORROPIA CODE
3547 ! *** SUBROUTINE CALCH6
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 !=======================================================================
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 ************************************************
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)
3586 PSI6HI = CHI6-TINY ! MIN(CHI6-TINY, CHI4)
3588 ! *** INITIAL VALUES FOR BISECTION ************************************
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)
3600 IF (SIGN(1.d0,Y1)*SIGN(1.d0,Y2).LT.ZERO) GOTO 20 ! (Y1*Y2.LT.ZERO)
3605 ! *** NO SUBDIVISION WITH SOLUTION; IF ABS(Y2)<EPS SOLUTION IS ASSUMED
3607 IF (ABS(Y2) .GT. EPS) Y2 = FUNCH6A (PSI6LO)
3610 ! *** PERFORM BISECTION ***********************************************
3615 IF (SIGN(1.d0,Y1)*SIGN(1.d0,Y3) .LE. ZERO) THEN ! (Y1*Y3 .LE. ZERO)
3622 IF (ABS(X2-X1) .LE. EPS*X1) GOTO 40
3624 CALL PUSHERR (0002, 'CALCH6') ! WARNING ERROR: NO CONVERGENCE
3626 ! *** CONVERGED ; RETURN **********************************************
3631 ! *** CALCULATE HSO4 SPECIATION AND RETURN *******************************
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
3643 ! *** END OF SUBROUTINE CALCH6 ******************************************
3645 END SUBROUTINE CALCH6
3650 !=======================================================================
3652 ! *** ISORROPIA CODE
3653 ! *** SUBROUTINE FUNCH6A
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)
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 ************************************************
3689 ! *** SOLVE EQUATIONS ; WITH ITERATIONS FOR ACTIVITY COEF. ************
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)
3711 PSI4 =0.5d0*(-BB - SQRT(DD))
3712 PSI4 = MIN(PSI4,CHI4)
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)
3730 GNH3 = MAX(CHI4 - PSI4, TINY)
3731 GHNO3 = MAX(CHI5 - PSI5, TINY)
3732 GHCL = MAX(CHI6 - PSI6, TINY)
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
3751 ! *** CALCULATE FUNCTION VALUE FOR OUTER LOOP ***************************
3753 20 FUNCH6A = MOLAL(3)*MOLAL(4)/GHCL/GNH3/A6/A4 - ONE
3757 ! *** END OF FUNCTION FUNCH6A *******************************************
3759 END FUNCTION FUNCH6A
3761 !=======================================================================
3763 ! *** ISORROPIA CODE
3764 ! *** SUBROUTINE CALCH5
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 !=======================================================================
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
3798 ! *** SETUP PARAMETERS ************************************************
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)
3812 PSI6HI = CHI6-TINY ! MIN(CHI6-TINY, CHI4)
3814 ! *** INITIAL VALUES FOR BISECTION ************************************
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)
3826 IF (SIGN(1.d0,Y1)*SIGN(1.d0,Y2).LT.ZERO) GOTO 20 ! (Y1*Y2.LT.ZERO)
3831 ! *** NO SUBDIVISION WITH SOLUTION; IF ABS(Y2)<EPS SOLUTION IS ASSUMED
3833 IF (ABS(Y2) .GT. EPS) Y2 = FUNCH5A (PSI6LO)
3836 ! *** PERFORM BISECTION ***********************************************
3841 IF (SIGN(1.d0,Y1)*SIGN(1.d0,Y3) .LE. ZERO) THEN ! (Y1*Y3 .LE. ZERO)
3848 IF (ABS(X2-X1) .LE. EPS*X1) GOTO 40
3850 CALL PUSHERR (0002, 'CALCH5') ! WARNING ERROR: NO CONVERGENCE
3852 ! *** CONVERGED ; RETURN **********************************************
3857 ! *** CALCULATE HSO4 SPECIATION AND RETURN *******************************
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
3869 ! *** END OF SUBROUTINE CALCH5 ******************************************
3871 END SUBROUTINE CALCH5
3876 !=======================================================================
3878 ! *** ISORROPIA CODE
3879 ! *** SUBROUTINE FUNCH5A
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)
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 ************************************************
3915 ! *** SOLVE EQUATIONS ; WITH ITERATIONS FOR ACTIVITY COEF. ************
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)
3937 PSI4 =0.5d0*(-BB - SQRT(DD))
3938 PSI4 = MIN(PSI4,CHI4)
3943 IF (CHI1.GT.TINY .AND. WATER.GT.TINY) THEN ! NA2SO4 DISSOLUTION
3947 CALL POLY3 (AA, BB, CC, PSI1, ISLV)
3949 PSI1 = MIN (PSI1, CHI1)
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
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)
3968 GNH3 = MAX(CHI4 - PSI4, TINY)
3969 GHNO3 = MAX(CHI5 - PSI5, TINY)
3970 GHCL = MAX(CHI6 - PSI6, TINY)
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
3989 ! *** CALCULATE FUNCTION VALUE FOR OUTER LOOP ***************************
3991 20 FUNCH5A = MOLAL(3)*MOLAL(4)/GHCL/GNH3/A6/A4 - ONE
3995 ! *** END OF FUNCTION FUNCH5A *******************************************
3997 END FUNCTION FUNCH5A
3999 !=======================================================================
4001 ! *** ISORROPIA CODE
4002 ! *** SUBROUTINE CALCH4
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 !=======================================================================
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
4036 ! *** SETUP PARAMETERS ************************************************
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)
4050 PSI6HI = CHI6-TINY ! MIN(CHI6-TINY, CHI4)
4052 ! *** INITIAL VALUES FOR BISECTION ************************************
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)
4064 IF (SIGN(1.d0,Y1)*SIGN(1.d0,Y2).LT.ZERO) GOTO 20 ! (Y1*Y2.LT.ZERO)
4069 ! *** NO SUBDIVISION WITH SOLUTION; IF ABS(Y2)<EPS SOLUTION IS ASSUMED
4071 IF (ABS(Y2) .GT. EPS) Y2 = FUNCH4A (PSI6LO)
4074 ! *** PERFORM BISECTION ***********************************************
4079 IF (SIGN(1.d0,Y1)*SIGN(1.d0,Y3) .LE. ZERO) THEN ! (Y1*Y3 .LE. ZERO)
4086 IF (ABS(X2-X1) .LE. EPS*X1) GOTO 40
4088 CALL PUSHERR (0002, 'CALCH4') ! WARNING ERROR: NO CONVERGENCE
4090 ! *** CONVERGED ; RETURN **********************************************
4095 ! *** CALCULATE HSO4 SPECIATION AND RETURN *******************************
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
4107 ! *** END OF SUBROUTINE CALCH4 ******************************************
4109 END SUBROUTINE CALCH4
4114 !=======================================================================
4116 ! *** ISORROPIA CODE
4117 ! *** SUBROUTINE FUNCH4A
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)
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 ************************************************
4153 ! *** SOLVE EQUATIONS ; WITH ITERATIONS FOR ACTIVITY COEF. ************
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)
4175 PSI4 =0.5d0*(-BB - SQRT(DD))
4176 PSI4 = MIN(PSI4,CHI4)
4181 IF (CHI1.GT.TINY .AND. WATER.GT.TINY) THEN ! NA2SO4 DISSOLUTION
4185 CALL POLY3 (AA, BB, CC, PSI1, ISLV)
4187 PSI1 = MIN (PSI1, CHI1)
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
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)
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)
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)
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
4229 ELSEIF (DELT-PSI32.GT.ZERO .AND. PSI32.GT.ZERO) THEN
4235 ! *** CALCULATE GAS / SOLID SPECIES (LIQUID IN MOLAL ALREADY) *********
4237 GNH3 = MAX(GNH3 - PSI3, TINY)
4238 GHCL = MAX(GHCL - PSI3, TINY)
4241 CALL CALCMR ! Water content
4243 ! *** CALCULATE ACTIVITIES OR TERMINATE INTERNAL LOOP *****************
4245 IF (FRST.AND.CALAOU .OR. .NOT.FRST.AND.CALAIN) THEN
4252 ! *** CALCULATE FUNCTION VALUE FOR OUTER LOOP ***************************
4254 20 FUNCH4A = MOLAL(3)*MOLAL(4)/GHCL/GNH3/A6/A4 - ONE
4258 ! *** END OF FUNCTION FUNCH4A *******************************************
4260 END FUNCTION FUNCH4A
4262 !=======================================================================
4264 ! *** ISORROPIA CODE
4265 ! *** SUBROUTINE CALCH3
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 !=======================================================================
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
4299 ! *** SETUP PARAMETERS ************************************************
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)
4313 PSI6HI = CHI6-TINY ! MIN(CHI6-TINY, CHI4)
4315 ! *** INITIAL VALUES FOR BISECTION ************************************
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)
4327 IF (SIGN(1.d0,Y1)*SIGN(1.d0,Y2).LT.ZERO) GOTO 20 ! (Y1*Y2.LT.ZERO)
4332 ! *** NO SUBDIVISION WITH SOLUTION; IF ABS(Y2)<EPS SOLUTION IS ASSUMED
4334 IF (ABS(Y2) .GT. EPS) Y2 = FUNCH3A (PSI6LO)
4337 ! *** PERFORM BISECTION ***********************************************
4342 IF (SIGN(1.d0,Y1)*SIGN(1.d0,Y3) .LE. ZERO) THEN ! (Y1*Y3 .LE. ZERO)
4349 IF (ABS(X2-X1) .LE. EPS*X1) GOTO 40
4351 CALL PUSHERR (0002, 'CALCH3') ! WARNING ERROR: NO CONVERGENCE
4353 ! *** CONVERGED ; RETURN **********************************************
4358 ! *** CALCULATE HSO4 SPECIATION AND RETURN *******************************
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
4370 ! *** END OF SUBROUTINE CALCH3 ******************************************
4372 END SUBROUTINE CALCH3
4377 !=======================================================================
4379 ! *** ISORROPIA CODE
4380 ! *** SUBROUTINE FUNCH3A
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)
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 ************************************************
4416 ! *** SOLVE EQUATIONS ; WITH ITERATIONS FOR ACTIVITY COEF. ************
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)
4438 PSI4 =0.5d0*(-BB - SQRT(DD))
4439 PSI4 = MIN(PSI4,CHI4)
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)
4450 IF (CHI1.GT.TINY .AND. WATER.GT.TINY) THEN ! NA2SO4 DISSOLUTION
4454 CALL POLY3 (AA, BB, CC, PSI1, ISLV)
4456 PSI1 = MIN (PSI1, CHI1)
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
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)
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)
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)
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
4498 ELSEIF (DELT-PSI32.GT.ZERO .AND. PSI32.GT.ZERO) THEN
4504 ! *** CALCULATE GAS / SOLID SPECIES (LIQUID IN MOLAL ALREADY) *********
4506 GNH3 = MAX(GNH3 - PSI3, TINY)
4507 GHCL = MAX(GHCL - PSI3, TINY)
4510 CALL CALCMR ! Water content
4512 ! *** CALCULATE ACTIVITIES OR TERMINATE INTERNAL LOOP *****************
4514 IF (FRST.AND.CALAOU .OR. .NOT.FRST.AND.CALAIN) THEN
4521 ! *** CALCULATE FUNCTION VALUE FOR OUTER LOOP ***************************
4523 20 FUNCH3A = MOLAL(3)*MOLAL(4)/GHCL/GNH3/A6/A4 - ONE
4527 ! *** END OF FUNCTION FUNCH3A *******************************************
4529 END FUNCTION FUNCH3A
4531 !=======================================================================
4533 ! *** ISORROPIA CODE
4534 ! *** SUBROUTINE CALCH2
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 !=======================================================================
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'
4568 SCASE = 'H2 ; SUBCASE 1'
4569 ELSE ! NO3 NON EXISTANT, WATER NOT POSSIBLE
4570 SCASE = 'H2 ; SUBCASE 1'
4572 SCASE = 'H2 ; SUBCASE 1'
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'
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 !=======================================================================
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 ************************************************
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)
4635 PSI6HI = CHI6-TINY ! MIN(CHI6-TINY, CHI4)
4637 ! *** INITIAL VALUES FOR BISECTION ************************************
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)
4649 IF (SIGN(1.d0,Y1)*SIGN(1.d0,Y2).LT.ZERO) GOTO 20 ! (Y1*Y2.LT.ZERO)
4654 ! *** NO SUBDIVISION WITH SOLUTION; IF ABS(Y2)<EPS SOLUTION IS ASSUMED
4656 IF (Y2 .GT. EPS) Y2 = FUNCH2A (PSI6LO)
4659 ! *** PERFORM BISECTION ***********************************************
4664 IF (SIGN(1.d0,Y1)*SIGN(1.d0,Y3) .LE. ZERO) THEN ! (Y1*Y3 .LE. ZERO)
4671 IF (ABS(X2-X1) .LE. EPS*X1) GOTO 40
4673 CALL PUSHERR (0002, 'CALCH2A') ! WARNING ERROR: NO CONVERGENCE
4675 ! *** CONVERGED ; RETURN **********************************************
4680 ! *** CALCULATE HSO4 SPECIATION AND RETURN *******************************
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
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)
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 ************************************************
4738 ! *** SOLVE EQUATIONS ; WITH ITERATIONS FOR ACTIVITY COEF. ************
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)
4762 PSI4 =0.5d0*(-BB - SQRT(DD))
4763 PSI4 = MIN(PSI4,CHI4)
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)
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)
4780 IF (CHI1.GT.TINY .AND. WATER.GT.TINY) THEN ! NA2SO4 DISSOLUTION
4784 CALL POLY3 (AA, BB, CC, PSI1, ISLV)
4786 PSI1 = MIN (PSI1, CHI1)
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)
4805 GNH3 = MAX(CHI4 - PSI4, TINY)
4806 GHNO3 = MAX(CHI5 - PSI5, TINY)
4807 GHCL = MAX(CHI6 - PSI6, TINY)
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)
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
4826 ELSEIF (DELT-PSI32.GT.ZERO .AND. PSI32.GT.ZERO) THEN
4832 ! *** CALCULATE GAS / SOLID SPECIES (LIQUID IN MOLAL ALREADY) *********
4834 GNH3 = MAX(GNH3 - PSI3, TINY)
4835 GHCL = MAX(GHCL - PSI3, TINY)
4838 CALL CALCMR ! Water content
4840 ! *** CALCULATE ACTIVITIES OR TERMINATE INTERNAL LOOP *****************
4842 IF (FRST.AND.CALAOU .OR. .NOT.FRST.AND.CALAIN) THEN
4849 ! *** CALCULATE FUNCTION VALUE FOR OUTER LOOP ***************************
4851 20 FUNCH2A = MOLAL(3)*MOLAL(4)/GHCL/GNH3/A64 - ONE
4855 ! *** END OF FUNCTION FUNCH2A *******************************************
4857 END FUNCTION FUNCH2A
4860 !=======================================================================
4862 ! *** ISORROPIA CODE
4863 ! *** SUBROUTINE CALCH1
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 !=======================================================================
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'
4895 SCASE = 'H1 ; SUBCASE 2' ! LIQUID & SOLID PHASE POSSIBLE
4896 CALL CALCMDRH (RH, DRMH1, DRNH4NO3, CALCH1A, CALCH2A)
4897 SCASE = 'H1 ; SUBCASE 2'
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 !=======================================================================
4927 IMPLICIT DOUBLE PRECISION (A-H, O-Z)
4928 ! INCLUDE 'ISRPIA.INC'
4929 DOUBLE PRECISION LAMDA, LAMDA1, LAMDA2, KAPA, KAPA1, KAPA2, NAFR, &
4932 ! *** CALCULATE NON VOLATILE SOLIDS ***********************************
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
4952 THETA1 = GAM - BET*(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
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
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
4988 ! SEPERATE SOLUTION OF NH4CL & NH4NO3 EQUILIBRIA
4992 DD1 = (ALF+BET)*(ALF+BET) - 4.0D0*(ALF*BET-A1)
4993 DD2 = (ALF+GAM)*(ALF+GAM) - 4.0D0*(ALF*GAM-A2)
4997 IF (DD1.GE.ZERO) THEN
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
5004 ELSE IF (KAPA2.GE.ZERO .AND. KAPA2.LE.MIN(ALF,BET)) THEN
5011 ! NH4NO3 EQUILIBRIUM
5013 IF (DD2.GE.ZERO) THEN
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
5020 ELSE IF (LAMDA2.GE.ZERO .AND. LAMDA2.LE.MIN(ALF,GAM)) THEN
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
5037 ! *** CALCULATE COMPOSITION OF VOLATILE SPECIES ***********************
5043 GNH3 = ALF - KAPA - LAMDA
5049 ! *** END OF SUBROUTINE CALCH1A *****************************************
5051 END SUBROUTINE CALCH1A
5052 !=======================================================================
5054 ! *** ISORROPIA CODE
5055 ! *** SUBROUTINE CALCI6
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 !=======================================================================
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 **********************************************
5083 ! *** SETUP PARAMETERS ************************************************
5085 CHI1 = CNH4HS4 ! Save from CALCI1 run
5091 PSI1 = CNH4HS4 ! ASSIGN INITIAL PSI's
5097 CALAOU = .TRUE. ! Outer loop activity calculation flag
5101 ! *** SOLVE EQUATIONS ; WITH ITERATIONS FOR ACTIVITY COEF. ************
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
5123 CNA2SO4 = CHI4 - PSI4
5126 CALL CALCMR ! Water content
5128 ! *** CALCULATE ACTIVITIES OR TERMINATE INTERNAL LOOP *****************
5130 IF (FRST.AND.CALAOU .OR. .NOT.FRST.AND.CALAIN) THEN
5139 ! *** END OF SUBROUTINE CALCI6 *****************************************
5141 END SUBROUTINE CALCI6
5143 !=======================================================================
5145 ! *** ISORROPIA CODE
5146 ! *** SUBROUTINE CALCI5
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 !=======================================================================
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 **********************************************
5174 ! *** SETUP PARAMETERS ************************************************
5176 CHI1 = CNH4HS4 ! Save from CALCI1 run
5182 PSI1 = CNH4HS4 ! ASSIGN INITIAL PSI's
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
5199 ! *** INITIAL VALUES FOR BISECTION ************************************
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)
5215 IF (SIGN(1.d0,Y1)*SIGN(1.d0,Y2).LT.ZERO) GOTO 20 ! (Y1*Y2.LT.ZERO)
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
5226 ELSE IF (ABS(Y2) .LT. EPS) THEN ! X2 IS A SOLUTION
5229 CALL PUSHERR (0001, 'CALCI5') ! WARNING ERROR: NO SOLUTION
5233 ! *** PERFORM BISECTION ***********************************************
5238 IF (SIGN(1.d0,Y1)*SIGN(1.d0,Y3) .LE. ZERO) THEN ! (Y1*Y3 .LE. ZERO)
5245 IF (ABS(X2-X1) .LE. EPS*X1) GOTO 40
5247 CALL PUSHERR (0002, 'CALCI5') ! WARNING ERROR: NO CONVERGENCE
5249 ! *** CONVERGED ; RETURN **********************************************
5256 ! *** END OF SUBROUTINE CALCI5 *****************************************
5258 END SUBROUTINE CALCI5
5263 !=======================================================================
5265 ! *** ISORROPIA CODE
5266 ! *** SUBROUTINE FUNCI5A
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)
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
5296 ! *** SOLVE EQUATIONS ; WITH ITERATIONS FOR ACTIVITY COEF. ************
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
5320 CNA2SO4 = CHI4 - PSI4
5323 CALL CALCMR ! Water content
5325 ! *** CALCULATE ACTIVITIES OR TERMINATE INTERNAL LOOP *****************
5327 IF (FRST.AND.CALAOU .OR. .NOT.FRST.AND.CALAIN) THEN
5334 ! *** CALCULATE OBJECTIVE FUNCTION ************************************
5336 20 A4 = XK5 *(WATER/GAMA(2))**3.0
5337 FUNCI5A= MOLAL(5)*MOLAL(2)*MOLAL(2)/A4 - ONE
5340 ! *** END OF FUNCTION FUNCI5A ********************************************
5342 END FUNCTION FUNCI5A
5343 !=======================================================================
5345 ! *** ISORROPIA CODE
5346 ! *** SUBROUTINE CALCI4
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 !=======================================================================
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 **********************************************
5374 ! *** SETUP PARAMETERS ************************************************
5376 CHI1 = CNH4HS4 ! Save from CALCI1 run
5382 PSI1 = CNH4HS4 ! ASSIGN INITIAL PSI's
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
5399 ! *** INITIAL VALUES FOR BISECTION ************************************
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)
5415 IF (SIGN(1.d0,Y1)*SIGN(1.d0,Y2).LT.ZERO) GOTO 20 ! (Y1*Y2.LT.ZERO)
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
5426 ELSE IF (ABS(Y2) .LT. EPS) THEN ! X2 IS A SOLUTION
5429 CALL PUSHERR (0001, 'CALCI4') ! WARNING ERROR: NO SOLUTION
5433 ! *** PERFORM BISECTION ***********************************************
5438 IF (SIGN(1.d0,Y1)*SIGN(1.d0,Y3) .LE. ZERO) THEN ! (Y1*Y3 .LE. ZERO)
5445 IF (ABS(X2-X1) .LE. EPS*X1) GOTO 40
5447 CALL PUSHERR (0002, 'CALCI4') ! WARNING ERROR: NO CONVERGENCE
5449 ! *** CONVERGED ; RETURN **********************************************
5456 ! *** END OF SUBROUTINE CALCI4 *****************************************
5458 END SUBROUTINE CALCI4
5463 !=======================================================================
5465 ! *** ISORROPIA CODE
5466 ! *** SUBROUTINE FUNCI4A
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)
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
5496 ! *** SOLVE EQUATIONS ; WITH ITERATIONS FOR ACTIVITY COEF. ************
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.
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
5524 CNA2SO4 = CHI4 - PSI4
5525 CNH42S4 = CHI5 - PSI5
5527 CALL CALCMR ! Water content
5529 ! *** CALCULATE ACTIVITIES OR TERMINATE INTERNAL LOOP *****************
5531 IF (FRST.AND.CALAOU .OR. .NOT.FRST.AND.CALAIN) THEN
5538 ! *** CALCULATE OBJECTIVE FUNCTION ************************************
5540 20 A4 = XK5 *(WATER/GAMA(2))**3.0
5541 FUNCI4A= MOLAL(5)*MOLAL(2)*MOLAL(2)/A4 - ONE
5544 ! *** END OF FUNCTION FUNCI4A ********************************************
5546 END FUNCTION FUNCI4A
5547 !=======================================================================
5549 ! *** ISORROPIA CODE
5550 ! *** SUBROUTINE CALCI3
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
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 !=======================================================================
5575 IMPLICIT DOUBLE PRECISION (A-H, O-Z)
5576 ! INCLUDE 'ISRPIA.INC'
5577 EXTERNAL CALCI1A, CALCI4
5579 ! *** FIND DRY COMPOSITION **********************************************
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'
5591 IF (WATER.LE.TINY) THEN
5592 IF (RH.LT.DRMI3) THEN ! SOLID SOLUTION
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'
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 !=======================================================================
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
5654 PSI1 = CNH4HS4 ! ASSIGN INITIAL PSI's
5660 CALAOU = .TRUE. ! Outer loop activity calculation flag
5661 PSI2LO = ZERO ! Low limit
5662 PSI2HI = CHI2 ! High limit
5664 ! *** INITIAL VALUES FOR BISECTION ************************************
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)
5678 X2 = MAX(X1-DX, PSI2LO)
5680 IF (SIGN(1.d0,Y1)*SIGN(1.d0,Y2).LT.ZERO) GOTO 20 ! (Y1*Y2.LT.ZERO)
5685 ! *** { YLO, YHI } > 0.0 THE SOLUTION IS ALWAYS SUPERSATURATED WITH LC
5687 IF (Y2.GT.EPS) Y2 = FUNCI3A (ZERO)
5690 ! *** PERFORM BISECTION ***********************************************
5695 IF (SIGN(1.d0,Y1)*SIGN(1.d0,Y3) .LE. ZERO) THEN ! (Y1*Y3 .LE. ZERO)
5702 IF (ABS(X2-X1) .LE. EPS*X1) GOTO 40
5704 CALL PUSHERR (0002, 'CALCI3A') ! WARNING ERROR: NO CONVERGENCE
5706 ! *** CONVERGED ; 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)
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)
5757 ! *** INITIAL VALUES FOR BISECTION ************************************
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)
5772 X2 = MAX(X1-DX, PSI4LO)
5774 IF (SIGN(1.d0,Y1)*SIGN(1.d0,Y2).LT.ZERO) GOTO 20 ! (Y1*Y2.LT.ZERO)
5779 ! *** { YLO, YHI } > 0.0 THE SOLUTION IS ALWAYS SUPERSATURATED WITH NA2SO4
5781 IF (Y2.GT.EPS) Y2 = FUNCI3B (PSI4LO)
5784 ! *** PERFORM BISECTION ***********************************************
5789 IF (SIGN(1.d0,Y1)*SIGN(1.d0,Y3) .LE. ZERO) THEN ! (Y1*Y3 .LE. ZERO)
5796 IF (ABS(X2-X1) .LE. EPS*X1) GOTO 40
5798 CALL PUSHERR (0004, 'FUNCI3A') ! WARNING ERROR: NO CONVERGENCE
5800 ! *** INNER LOOP CONVERGED **********************************************
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
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)
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 ************************************************
5846 ! *** SETUP PARAMETERS ************************************************
5851 ! *** SOLVE EQUATIONS ; WITH ITERATIONS FOR ACTIVITY COEF. ************
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.
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)
5879 CNA2SO4 = MAX(CHI4 - PSI4, ZERO)
5880 CNH42S4 = MAX(CHI5 - PSI5, ZERO)
5882 CALL CALCMR ! Water content
5884 ! *** CALCULATE ACTIVITIES OR TERMINATE INTERNAL LOOP *****************
5886 IF (FRST.AND.CALAOU .OR. .NOT.FRST.AND.CALAIN) THEN
5893 ! *** CALCULATE OBJECTIVE FUNCTION ************************************
5895 20 A4 = XK5 *(WATER/GAMA(2))**3.0
5896 FUNCI3B= MOLAL(5)*MOLAL(2)*MOLAL(2)/A4 - ONE
5899 ! *** END OF FUNCTION FUNCI3B ********************************************
5901 END FUNCTION FUNCI3B
5902 !=======================================================================
5904 ! *** ISORROPIA CODE
5905 ! *** SUBROUTINE CALCI2
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
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 !=======================================================================
5930 IMPLICIT DOUBLE PRECISION (A-H, O-Z)
5931 ! INCLUDE 'ISRPIA.INC'
5932 EXTERNAL CALCI1A, CALCI3A
5934 ! *** FIND DRY COMPOSITION **********************************************
5938 ! *** REGIME DEPENDS UPON THE POSSIBLE SOLIDS & RH **********************
5940 IF (CNH4HS4.GT.TINY) THEN
5941 SCASE = 'I2 ; SUBCASE 1'
5943 SCASE = 'I2 ; SUBCASE 1'
5946 IF (WATER.LE.TINY) THEN
5947 IF (RH.LT.DRMI2) THEN ! SOLID SOLUTION ONLY
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'
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 !=======================================================================
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
6008 PSI1 = CNH4HS4 ! ASSIGN INITIAL PSI's
6014 CALAOU = .TRUE. ! Outer loop activity calculation flag
6015 PSI2LO = ZERO ! Low limit
6016 PSI2HI = CHI2 ! High limit
6018 ! *** INITIAL VALUES FOR BISECTION ************************************
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)
6032 X2 = MAX(X1-DX, PSI2LO)
6034 IF (SIGN(1.d0,Y1)*SIGN(1.d0,Y2).LT.ZERO) GOTO 20 ! (Y1*Y2.LT.ZERO)
6039 ! *** { YLO, YHI } > 0.0 THE SOLUTION IS ALWAYS SUPERSATURATED WITH LC
6041 IF (Y2.GT.EPS) Y2 = FUNCI3A (ZERO)
6044 ! *** PERFORM BISECTION ***********************************************
6049 IF (SIGN(1.d0,Y1)*SIGN(1.d0,Y3) .LE. ZERO) THEN ! (Y1*Y3 .LE. ZERO)
6056 IF (ABS(X2-X1) .LE. EPS*X1) GOTO 40
6058 CALL PUSHERR (0002, 'CALCI2A') ! WARNING ERROR: NO CONVERGENCE
6060 ! *** CONVERGED ; 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)
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 ************************************************
6105 PSI2 = P2 ! Save PSI2 in COMMON BLOCK
6111 ! *** SOLVE EQUATIONS ; WITH ITERATIONS FOR ACTIVITY COEF. ************
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.
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)
6128 IF (CHI4.GT.TINY .AND. WATER.GT.TINY) THEN
6129 AA = PSI2+PSI5+PSI6+PSI3
6131 CC = 0.25D0*(PSI3*PSI3*(PSI2+PSI5+PSI6)-A4)
6132 CALL POLY3 (AA, BB, CC, PSI4, ISLV)
6134 PSI4 = MIN (PSI4, CHI4)
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
6144 CALL POLY3 (AA, BB, CC, PSI3, ISLV)
6146 PSI3 = MIN (PSI3, CHI3)
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
6165 CNAHSO4 = CHI3 - PSI3
6166 CNA2SO4 = CHI4 - PSI4
6167 CNH42S4 = CHI5 - PSI5
6169 CALL CALCMR ! Water content
6171 ! *** CALCULATE ACTIVITIES OR TERMINATE INTERNAL LOOP *****************
6173 IF (FRST.AND.CALAOU .OR. .NOT.FRST.AND.CALAIN) THEN
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
6186 ! *** END OF FUNCTION FUNCI2A *******************************************
6188 END FUNCTION FUNCI2A
6190 !=======================================================================
6192 ! *** ISORROPIA CODE
6193 ! *** SUBROUTINE CALCI1
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 !=======================================================================
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'
6225 SCASE = 'I1 ; SUBCASE 2' ! LIQUID & SOLID PHASE POSSIBLE
6226 CALL CALCMDRH (RH, DRMI1, DRNH4HS4, CALCI1A, CALCI2A)
6227 SCASE = 'I1 ; SUBCASE 2'
6230 ! *** AMMONIA IN GAS PHASE **********************************************
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 !=======================================================================
6261 IMPLICIT DOUBLE PRECISION (A-H, O-Z)
6262 ! INCLUDE 'ISRPIA.INC'
6264 ! *** CALCULATE NON VOLATILE SOLIDS ***********************************
6266 CNA2SO4 = 0.5D0*W(1)
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)
6290 ! *** CALCULATE GAS SPECIES *********************************************
6298 ! *** END OF SUBROUTINE CALCI1A *****************************************
6300 END SUBROUTINE CALCI1A
6301 !=======================================================================
6303 ! *** ISORROPIA CODE
6304 ! *** SUBROUTINE CALCJ3
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 !=======================================================================
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
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
6335 PSI2 = CHI2 ! ALL NH4HSO4 DELIQUESCED
6337 ! *** SOLVE EQUATIONS ; WITH ITERATIONS FOR ACTIVITY COEF. ************
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)
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
6363 CALL CALCMR ! Water content
6365 ! *** CALCULATE ACTIVITIES OR TERMINATE INTERNAL LOOP *****************
6367 IF (FRST.AND.CALAOU .OR. .NOT.FRST.AND.CALAIN) THEN
6376 ! *** END OF SUBROUTINE CALCJ3 ******************************************
6378 END SUBROUTINE CALCJ3
6379 !=======================================================================
6381 ! *** ISORROPIA CODE
6382 ! *** SUBROUTINE CALCJ2
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 !=======================================================================
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, &
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 ************************************
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)
6431 IF (SIGN(1.d0,Y1)*SIGN(1.d0,Y2).LT.ZERO) GOTO 20 ! (Y1*Y2.LT.ZERO)
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
6442 ELSE IF (ABS(Y2) .LT. EPS) THEN ! X2 IS A SOLUTION
6445 CALL PUSHERR (0001, 'CALCJ2') ! WARNING ERROR: NO SOLUTION
6449 ! *** PERFORM BISECTION ***********************************************
6454 IF (SIGN(1.d0,Y1)*SIGN(1.d0,Y3) .LE. ZERO) THEN ! (Y1*Y3 .LE. ZERO)
6461 IF (ABS(X2-X1) .LE. EPS*X1) GOTO 40
6463 CALL PUSHERR (0002, 'CALCJ2') ! WARNING ERROR: NO CONVERGENCE
6465 ! *** CONVERGED ; RETURN **********************************************
6472 ! *** END OF SUBROUTINE CALCJ2 ******************************************
6474 END SUBROUTINE CALCJ2
6479 !=======================================================================
6481 ! *** ISORROPIA CODE
6482 ! *** SUBROUTINE FUNCJ2
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)
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, &
6507 ! *** SETUP PARAMETERS ************************************************
6512 LAMDA = MAX(W(2) - W(3) - W(1), TINY) ! FREE H2SO4
6514 PSI2 = CHI2 ! ALL NH4HSO4 DELIQUESCED
6516 ! *** SOLVE EQUATIONS ; WITH ITERATIONS FOR ACTIVITY COEF. ************
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)
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)
6543 CALL CALCMR ! Water content
6545 ! *** CALCULATE ACTIVITIES OR TERMINATE INTERNAL LOOP *****************
6547 IF (FRST.AND.CALAOU .OR. .NOT.FRST.AND.CALAIN) THEN
6554 ! *** CALCULATE OBJECTIVE FUNCTION ************************************
6556 20 FUNCJ2 = MOLAL(2)*MOLAL(6)/A1 - ONE
6558 ! *** END OF FUNCTION FUNCJ2 *******************************************
6562 !=======================================================================
6564 ! *** ISORROPIA CODE
6565 ! *** SUBROUTINE CALCJ1
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 !=======================================================================
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, &
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 ************************************
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)
6615 IF (SIGN(1.d0,Y1)*SIGN(1.d0,Y2).LT.ZERO) GOTO 20 ! (Y1*Y2.LT.ZERO)
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
6626 ELSE IF (ABS(Y2) .LT. EPS) THEN ! X2 IS A SOLUTION
6629 CALL PUSHERR (0001, 'CALCJ1') ! WARNING ERROR: NO SOLUTION
6633 ! *** PERFORM BISECTION ***********************************************
6638 IF (SIGN(1.d0,Y1)*SIGN(1.d0,Y3) .LE. ZERO) THEN ! (Y1*Y3 .LE. ZERO)
6645 IF (ABS(X2-X1) .LE. EPS*X1) GOTO 40
6647 CALL PUSHERR (0002, 'CALCJ1') ! WARNING ERROR: NO CONVERGENCE
6649 ! *** CONVERGED ; RETURN **********************************************
6656 ! *** END OF SUBROUTINE CALCJ1 ******************************************
6658 END SUBROUTINE CALCJ1
6663 !=======================================================================
6665 ! *** ISORROPIA CODE
6666 ! *** SUBROUTINE FUNCJ1
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)
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, &
6690 ! *** SETUP PARAMETERS ************************************************
6695 LAMDA = MAX(W(2) - W(3) - W(1), TINY) ! FREE H2SO4
6698 ! *** SOLVE EQUATIONS ; WITH ITERATIONS FOR ACTIVITY COEF. ************
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)
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
6720 MOLAL (5) = KAPA ! SO4I
6721 MOLAL (6) = LAMDA + PSI1 + PSI2 - KAPA ! HSO4I
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
6738 ! *** CALCULATE OBJECTIVE FUNCTION ************************************
6740 20 FUNCJ1 = MOLAL(2)*MOLAL(6)/A1 - ONE
6742 ! *** END OF FUNCTION FUNCJ1 *******************************************