!======================================================================= ! ! *** ISORROPIA CODE ! *** SUBROUTINE ISRP1F ! *** THIS SUBROUTINE IS THE DRIVER ROUTINE FOR THE FOREWARD PROBLEM OF ! AN AMMONIUM-SULFATE AEROSOL SYSTEM. ! THE COMPOSITION REGIME IS DETERMINED BY THE SULFATE RATIO AND BY ! THE AMBIENT RELATIVE HUMIDITY. ! ! *** COPYRIGHT 1996-2006, UNIVERSITY OF MIAMI, CARNEGIE MELLON UNIVERSITY, ! *** GEORGIA INSTITUTE OF TECHNOLOGY ! *** WRITTEN BY ATHANASIOS NENES ! *** UPDATED BY CHRISTOS FOUNTOUKIS ! ! REVISION HISTORY: * ! Original code was provided by Dr. ATHANASIOS NENES, Georgia Tech, in 2000 ! Revised by Y. Zhang, AER, Inc. to incorporate v1.5 into MADRID, 2000 ! Revised by Y. Zhang and Xiao-Ming Hu to incorporate it along with MADRID into WRF/Chem, 2005 ! Updated by Xiao-Ming Hu and Y. Zhang, NCSU to v. 1.7, Oct., 2007 !======================================================================= ! SUBROUTINE ISRP1F (WI, RHI, TEMPI) USE ISRPIA IMPLICIT DOUBLE PRECISION (A-H, O-Z) ! INCLUDE 'ISRPIA.INC' DIMENSION WI(NCOMP) ! ! *** INITIALIZE ALL VARIABLES IN COMMON BLOCK ************************** ! CALL INIT1 (WI, RHI, TEMPI) ! ! *** CALCULATE SULFATE RATIO ******************************************* ! SULRAT = W(3)/W(2) ! ! *** FIND CALCULATION REGIME FROM (SULRAT,RH) ************************** ! ! *** SULFATE POOR ! IF (2.0.LE.SULRAT) THEN DC = W(3) - 2.001D0*W(2) ! For numerical stability W(3) = W(3) + MAX(-DC, ZERO) ! IF(METSTBL.EQ.1) THEN SCASE = 'A2' CALL CALCA2 ! Only liquid (metastable) ELSE ! IF (RH.LT.DRNH42S4) THEN SCASE = 'A1' CALL CALCA1 ! NH42SO4 ; case A1 ! ELSEIF (DRNH42S4.LE.RH) THEN SCASE = 'A2' CALL CALCA2 ! Only liquid ; case A2 ENDIF ENDIF ! ! *** SULFATE RICH (NO ACID) ! ELSEIF (1.0.LE.SULRAT .AND. SULRAT.LT.2.0) THEN ! IF(METSTBL.EQ.1) THEN SCASE = 'B4' CALL CALCB4 ! Only liquid (metastable) ELSE ! IF (RH.LT.DRNH4HS4) THEN SCASE = 'B1' CALL CALCB1 ! NH4HSO4,LC,NH42SO4 ; case B1 ! ELSEIF (DRNH4HS4.LE.RH .AND. RH.LT.DRLC) THEN SCASE = 'B2' CALL CALCB2 ! LC,NH42S4 ; case B2 ! ELSEIF (DRLC.LE.RH .AND. RH.LT.DRNH42S4) THEN SCASE = 'B3' CALL CALCB3 ! NH42S4 ; case B3 ! ELSEIF (DRNH42S4.LE.RH) THEN SCASE = 'B4' CALL CALCB4 ! Only liquid ; case B4 ENDIF ENDIF CALL CALCNH3 ! ! *** SULFATE RICH (FREE ACID) ! ELSEIF (SULRAT.LT.1.0) THEN ! IF(METSTBL.EQ.1) THEN SCASE = 'C2' CALL CALCC2 ! Only liquid (metastable) ELSE ! IF (RH.LT.DRNH4HS4) THEN SCASE = 'C1' CALL CALCC1 ! NH4HSO4 ; case C1 ! ELSEIF (DRNH4HS4.LE.RH) THEN SCASE = 'C2' CALL CALCC2 ! Only liquid ; case C2 ! ENDIF ENDIF CALL CALCNH3 ENDIF ! ! *** RETURN POINT ! RETURN ! ! *** END OF SUBROUTINE ISRP1F ***************************************** ! END SUBROUTINE ISRP1F !======================================================================= ! ! *** ISORROPIA CODE ! *** SUBROUTINE ISRP2F ! *** THIS SUBROUTINE IS THE DRIVER ROUTINE FOR THE FOREWARD PROBLEM OF ! AN AMMONIUM-SULFATE-NITRATE AEROSOL SYSTEM. ! THE COMPOSITION REGIME IS DETERMINED BY THE SULFATE RATIO AND BY ! THE AMBIENT RELATIVE HUMIDITY. ! ! *** COPYRIGHT 1996-2006, UNIVERSITY OF MIAMI, CARNEGIE MELLON UNIVERSITY, ! *** GEORGIA INSTITUTE OF TECHNOLOGY ! *** WRITTEN BY ATHANASIOS NENES ! *** UPDATED BY CHRISTOS FOUNTOUKIS ! !======================================================================= ! SUBROUTINE ISRP2F (WI, RHI, TEMPI) USE ISRPIA IMPLICIT DOUBLE PRECISION (A-H, O-Z) ! INCLUDE 'ISRPIA.INC' DIMENSION WI(NCOMP) ! ! *** INITIALIZE ALL VARIABLES IN COMMON BLOCK ************************** ! CALL INIT2 (WI, RHI, TEMPI) ! ! *** CALCULATE SULFATE RATIO ******************************************* ! SULRAT = W(3)/W(2) ! ! *** FIND CALCULATION REGIME FROM (SULRAT,RH) ************************** ! ! *** SULFATE POOR ! IF (2.0.LE.SULRAT) THEN ! IF(METSTBL.EQ.1) THEN SCASE = 'D3' CALL CALCD3 ! Only liquid (metastable) ELSE ! IF (RH.LT.DRNH4NO3) THEN SCASE = 'D1' CALL CALCD1 ! NH42SO4,NH4NO3 ; case D1 ! ELSEIF (DRNH4NO3.LE.RH .AND. RH.LT.DRNH42S4) THEN SCASE = 'D2' CALL CALCD2 ! NH42S4 ; case D2 ! ELSEIF (DRNH42S4.LE.RH) THEN SCASE = 'D3' CALL CALCD3 ! Only liquid ; case D3 ENDIF ENDIF ! ! *** SULFATE RICH (NO ACID) ! FOR SOLVING THIS CASE, NITRIC ACID IS ASSUMED A MINOR SPECIES, ! THAT DOES NOT SIGNIFICANTLY PERTURB THE HSO4-SO4 EQUILIBRIUM. ! SUBROUTINES CALCB? ARE CALLED, AND THEN THE NITRIC ACID IS DISSOLVED ! FROM THE HNO3(G) -> (H+) + (NO3-) EQUILIBRIUM. ! ELSEIF (1.0.LE.SULRAT .AND. SULRAT.LT.2.0) THEN ! IF(METSTBL.EQ.1) THEN SCASE = 'B4' CALL CALCB4 ! Only liquid (metastable) SCASE = 'E4' ELSE ! IF (RH.LT.DRNH4HS4) THEN SCASE = 'B1' CALL CALCB1 ! NH4HSO4,LC,NH42SO4 ; case E1 SCASE = 'E1' ! ELSEIF (DRNH4HS4.LE.RH .AND. RH.LT.DRLC) THEN SCASE = 'B2' CALL CALCB2 ! LC,NH42S4 ; case E2 SCASE = 'E2' ! ELSEIF (DRLC.LE.RH .AND. RH.LT.DRNH42S4) THEN SCASE = 'B3' CALL CALCB3 ! NH42S4 ; case E3 SCASE = 'E3' ! ELSEIF (DRNH42S4.LE.RH) THEN SCASE = 'B4' CALL CALCB4 ! Only liquid ; case E4 SCASE = 'E4' ENDIF ENDIF ! CALL CALCNA ! HNO3(g) DISSOLUTION ! ! *** SULFATE RICH (FREE ACID) ! FOR SOLVING THIS CASE, NITRIC ACID IS ASSUMED A MINOR SPECIES, ! THAT DOES NOT SIGNIFICANTLY PERTURB THE HSO4-SO4 EQUILIBRIUM ! SUBROUTINE CALCC? IS CALLED, AND THEN THE NITRIC ACID IS DISSOLVED ! FROM THE HNO3(G) -> (H+) + (NO3-) EQUILIBRIUM. ! ELSEIF (SULRAT.LT.1.0) THEN ! IF(METSTBL.EQ.1) THEN SCASE = 'C2' CALL CALCC2 ! Only liquid (metastable) SCASE = 'F2' ELSE ! IF (RH.LT.DRNH4HS4) THEN SCASE = 'C1' CALL CALCC1 ! NH4HSO4 ; case F1 SCASE = 'F1' ! ELSEIF (DRNH4HS4.LE.RH) THEN SCASE = 'C2' CALL CALCC2 ! Only liquid ; case F2 SCASE = 'F2' ENDIF ENDIF ! CALL CALCNA ! HNO3(g) DISSOLUTION ENDIF ! ! *** RETURN POINT ! RETURN ! ! *** END OF SUBROUTINE ISRP2F ***************************************** ! END SUBROUTINE ISRP2F !======================================================================= ! ! *** ISORROPIA CODE ! *** SUBROUTINE ISRP3F ! *** THIS SUBROUTINE IS THE DRIVER ROUTINE FOR THE FORWARD PROBLEM OF ! AN AMMONIUM-SULFATE-NITRATE-CHLORIDE-SODIUM AEROSOL SYSTEM. ! THE COMPOSITION REGIME IS DETERMINED BY THE SULFATE & SODIUM ! RATIOS AND BY THE AMBIENT RELATIVE HUMIDITY. ! ! *** COPYRIGHT 1996-2006, UNIVERSITY OF MIAMI, CARNEGIE MELLON UNIVERSITY, ! *** GEORGIA INSTITUTE OF TECHNOLOGY ! *** WRITTEN BY ATHANASIOS NENES ! *** UPDATED BY CHRISTOS FOUNTOUKIS ! !======================================================================= ! SUBROUTINE ISRP3F (WI, RHI, TEMPI) USE ISRPIA IMPLICIT DOUBLE PRECISION (A-H, O-Z) ! INCLUDE 'ISRPIA.INC' DIMENSION WI(NCOMP) ! ! *** ADJUST FOR TOO LITTLE AMMONIUM AND CHLORIDE *********************** ! WI(3) = MAX (WI(3), 1.D-10) ! NH4+ : 1e-4 umoles/m3 WI(5) = MAX (WI(5), 1.D-10) ! Cl- : 1e-4 umoles/m3 ! ! *** ADJUST FOR TOO LITTLE SODIUM, SULFATE AND NITRATE COMBINED ******** ! IF (WI(1)+WI(2)+WI(4) .LE. 1d-10) THEN WI(1) = 1.D-10 ! Na+ : 1e-4 umoles/m3 WI(2) = 1.D-10 ! SO4- : 1e-4 umoles/m3 ENDIF ! ! *** INITIALIZE ALL VARIABLES IN COMMON BLOCK ************************** ! CALL ISOINIT3 (WI, RHI, TEMPI) ! ! *** CHECK IF TOO MUCH SODIUM ; ADJUST AND ISSUE ERROR MESSAGE ********* ! REST = 2.D0*W(2) + W(4) + W(5) IF (W(1).GT.REST) THEN ! NA > 2*SO4+CL+NO3 ? W(1) = (ONE-1D-6)*REST ! Adjust Na amount CALL PUSHERR (0050, 'ISRP3F') ! Warning error: Na adjusted ENDIF ! ! *** CALCULATE SULFATE & SODIUM RATIOS ********************************* ! SULRAT = (W(1)+W(3))/W(2) SODRAT = W(1)/W(2) ! ! *** FIND CALCULATION REGIME FROM (SULRAT,RH) ************************** ! ! *** SULFATE POOR ; SODIUM POOR ! IF (2.0.LE.SULRAT .AND. SODRAT.LT.2.0) THEN ! IF(METSTBL.EQ.1) THEN SCASE = 'G5' CALL CALCG5 ! Only liquid (metastable) ELSE ! IF (RH.LT.DRNH4NO3) THEN SCASE = 'G1' CALL CALCG1 ! NH42SO4,NH4NO3,NH4CL,NA2SO4 ! ELSEIF (DRNH4NO3.LE.RH .AND. RH.LT.DRNH4CL) THEN SCASE = 'G2' CALL CALCG2 ! NH42SO4,NH4CL,NA2SO4 ! ELSEIF (DRNH4CL.LE.RH .AND. RH.LT.DRNH42S4) THEN SCASE = 'G3' CALL CALCG3 ! NH42SO4,NA2SO4 ! ELSEIF (DRNH42S4.LE.RH .AND. RH.LT.DRNA2SO4) THEN SCASE = 'G4' CALL CALCG4 ! NA2SO4 ! ELSEIF (DRNA2SO4.LE.RH) THEN SCASE = 'G5' CALL CALCG5 ! Only liquid ENDIF ENDIF ! ! *** SULFATE POOR ; SODIUM RICH ! ELSE IF (SULRAT.GE.2.0 .AND. SODRAT.GE.2.0) THEN ! IF(METSTBL.EQ.1) THEN SCASE = 'H6' CALL CALCH6 ! Only liquid (metastable) ELSE ! IF (RH.LT.DRNH4NO3) THEN SCASE = 'H1' CALL CALCH1 ! NH4NO3,NH4CL,NA2SO4,NACL,NANO3 ! ELSEIF (DRNH4NO3.LE.RH .AND. RH.LT.DRNANO3) THEN SCASE = 'H2' CALL CALCH2 ! NH4CL,NA2SO4,NACL,NANO3 ! ELSEIF (DRNANO3.LE.RH .AND. RH.LT.DRNACL) THEN SCASE = 'H3' CALL CALCH3 ! NH4CL,NA2SO4,NACL ! ELSEIF (DRNACL.LE.RH .AND. RH.LT.DRNH4Cl) THEN SCASE = 'H4' CALL CALCH4 ! NH4CL,NA2SO4 ! ELSEIF (DRNH4Cl.LE.RH .AND. RH.LT.DRNA2SO4) THEN SCASE = 'H5' CALL CALCH5 ! NA2SO4 ! ELSEIF (DRNA2SO4.LE.RH) THEN SCASE = 'H6' CALL CALCH6 ! NO SOLID ENDIF ENDIF ! ! *** SULFATE RICH (NO ACID) ! ELSEIF (1.0.LE.SULRAT .AND. SULRAT.LT.2.0) THEN ! IF(METSTBL.EQ.1) THEN SCASE = 'I6' CALL CALCI6 ! Only liquid (metastable) ELSE ! IF (RH.LT.DRNH4HS4) THEN SCASE = 'I1' CALL CALCI1 ! NA2SO4,(NH4)2SO4,NAHSO4,NH4HSO4,LC ! ELSEIF (DRNH4HS4.LE.RH .AND. RH.LT.DRNAHSO4) THEN SCASE = 'I2' CALL CALCI2 ! NA2SO4,(NH4)2SO4,NAHSO4,LC ! ELSEIF (DRNAHSO4.LE.RH .AND. RH.LT.DRLC) THEN SCASE = 'I3' CALL CALCI3 ! NA2SO4,(NH4)2SO4,LC ! ELSEIF (DRLC.LE.RH .AND. RH.LT.DRNH42S4) THEN SCASE = 'I4' CALL CALCI4 ! NA2SO4,(NH4)2SO4 ! ELSEIF (DRNH42S4.LE.RH .AND. RH.LT.DRNA2SO4) THEN SCASE = 'I5' CALL CALCI5 ! NA2SO4 ! ELSEIF (DRNA2SO4.LE.RH) THEN SCASE = 'I6' CALL CALCI6 ! NO SOLIDS ENDIF ENDIF ! CALL CALCNHA ! MINOR SPECIES: HNO3, HCl CALL CALCNH3 ! NH3 ! ! *** SULFATE RICH (FREE ACID) ! ELSEIF (SULRAT.LT.1.0) THEN ! IF(METSTBL.EQ.1) THEN SCASE = 'J3' CALL CALCJ3 ! Only liquid (metastable) ELSE ! IF (RH.LT.DRNH4HS4) THEN SCASE = 'J1' CALL CALCJ1 ! NH4HSO4,NAHSO4 ! ELSEIF (DRNH4HS4.LE.RH .AND. RH.LT.DRNAHSO4) THEN SCASE = 'J2' CALL CALCJ2 ! NAHSO4 ! ELSEIF (DRNAHSO4.LE.RH) THEN SCASE = 'J3' CALL CALCJ3 ENDIF ENDIF ! CALL CALCNHA ! MINOR SPECIES: HNO3, HCl CALL CALCNH3 ! NH3 ENDIF ! ! *** RETURN POINT ! RETURN ! ! *** END OF SUBROUTINE ISRP3F ***************************************** ! END SUBROUTINE ISRP3F !======================================================================= ! ! *** ISORROPIA CODE ! *** SUBROUTINE CALCA2 ! *** CASE A2 ! ! THE MAIN CHARACTERISTICS OF THIS REGIME ARE: ! 1. SULFATE POOR (SULRAT > 2.0) ! 2. LIQUID AEROSOL PHASE ONLY POSSIBLE ! ! FOR CALCULATIONS, A BISECTION IS PERFORMED TOWARDS X, THE ! AMOUNT OF HYDROGEN IONS (H+) FOUND IN THE LIQUID PHASE. ! FOR EACH ESTIMATION OF H+, FUNCTION FUNCB2A CALCULATES THE ! CONCENTRATION OF IONS FROM THE NH3(GAS) - NH4+(LIQ) EQUILIBRIUM. ! ELECTRONEUTRALITY IS USED AS THE OBJECTIVE FUNCTION. ! ! *** COPYRIGHT 1996-2006, UNIVERSITY OF MIAMI, CARNEGIE MELLON UNIVERSITY, ! *** GEORGIA INSTITUTE OF TECHNOLOGY ! *** WRITTEN BY ATHANASIOS NENES ! *** UPDATED BY CHRISTOS FOUNTOUKIS ! !======================================================================= ! SUBROUTINE CALCA2 USE ISRPIA IMPLICIT DOUBLE PRECISION (A-H, O-Z) ! INCLUDE 'ISRPIA.INC' ! ! *** SETUP PARAMETERS ************************************************ ! CALAOU =.TRUE. ! Outer loop activity calculation flag OMELO = TINY ! Low limit: SOLUTION IS VERY BASIC OMEHI = 2.0D0*W(2) ! High limit: FROM NH4+ -> NH3(g) + H+(aq) ! ! *** CALCULATE WATER CONTENT ***************************************** ! MOLAL(5) = W(2) MOLAL(6) = ZERO CALL CALCMR ! ! *** INITIAL VALUES FOR BISECTION ************************************ ! X1 = OMEHI Y1 = FUNCA2 (X1) IF (ABS(Y1).LE.EPS) RETURN ! ! *** ROOT TRACKING ; FOR THE RANGE OF HI AND LO ********************** ! DX = (OMEHI-OMELO)/REAL(NDIV) DO 10 I=1,NDIV X2 = MAX(X1-DX, OMELO) Y2 = FUNCA2 (X2) IF (SIGN(1.d0,Y1)*SIGN(1.d0,Y2).LT.ZERO) GOTO 20 ! (Y1*Y2.LT.ZERO) X1 = X2 Y1 = Y2 10 CONTINUE IF (ABS(Y2).LE.EPS) THEN RETURN ELSE CALL PUSHERR (0001, 'CALCA2') ! WARNING ERROR: NO SOLUTION RETURN ENDIF ! ! *** PERFORM BISECTION *********************************************** ! 20 DO 30 I=1,MAXIT X3 = 0.5*(X1+X2) Y3 = FUNCA2 (X3) IF (SIGN(1.d0,Y1)*SIGN(1.d0,Y3) .LE. ZERO) THEN ! (Y1*Y3 .LE. ZERO) Y2 = Y3 X2 = X3 ELSE Y1 = Y3 X1 = X3 ENDIF IF (ABS(X2-X1) .LE. EPS*X1) GOTO 40 30 CONTINUE CALL PUSHERR (0002, 'CALCA2') ! WARNING ERROR: NO CONVERGENCE ! ! *** CONVERGED ; RETURN ********************************************** ! 40 X3 = 0.5*(X1+X2) Y3 = FUNCA2 (X3) RETURN ! ! *** END OF SUBROUTINE CALCA2 **************************************** ! END SUBROUTINE CALCA2 !======================================================================= ! ! *** ISORROPIA CODE ! *** FUNCTION FUNCA2 ! *** CASE A2 ! FUNCTION THAT SOLVES THE SYSTEM OF EQUATIONS FOR CASE A2 ; ! AND RETURNS THE VALUE OF THE ZEROED FUNCTION IN FUNCA2. ! !======================================================================= ! DOUBLE PRECISION FUNCTION FUNCA2 (OMEGI) USE ISRPIA IMPLICIT DOUBLE PRECISION (A-H, O-Z) ! INCLUDE 'ISRPIA.INC' DOUBLE PRECISION LAMDA ! ! *** SETUP PARAMETERS ************************************************ ! FRST = .TRUE. CALAIN = .TRUE. PSI = W(2) ! INITIAL AMOUNT OF (NH4)2SO4 IN SOLUTION ! ! *** SOLVE EQUATIONS ; WITH ITERATIONS FOR ACTIVITY COEF. ************ ! DO 10 I=1,NSWEEP A1 = XK1*WATER/GAMA(7)*(GAMA(8)/GAMA(7))**2. A2 = XK2*R*TEMP/XKW*(GAMA(8)/GAMA(9))**2. A3 = XKW*RH*WATER*WATER ! LAMDA = PSI/(A1/OMEGI+ONE) ZETA = A3/OMEGI ! ! *** SPECIATION & WATER CONTENT *************************************** ! MOLAL (1) = OMEGI ! HI MOLAL (3) = W(3)/(ONE/A2/OMEGI + ONE) ! NH4I MOLAL (5) = MAX(PSI-LAMDA,TINY) ! SO4I MOLAL (6) = LAMDA ! HSO4I GNH3 = MAX (W(3)-MOLAL(3), TINY) ! NH3GI COH = ZETA ! OHI ! ! *** CALCULATE ACTIVITIES OR TERMINATE INTERNAL LOOP ***************** ! IF (FRST.AND.CALAOU .OR. .NOT.FRST.AND.CALAIN) THEN CALL CALCACT ELSE GOTO 20 ENDIF 10 CONTINUE ! ! *** CALCULATE OBJECTIVE FUNCTION ************************************ ! 20 DENOM = (2.0*MOLAL(5)+MOLAL(6)) FUNCA2= (MOLAL(3)/DENOM - ONE) + MOLAL(1)/DENOM RETURN ! ! *** END OF FUNCTION FUNCA2 ******************************************** ! END FUNCTION FUNCA2 !======================================================================= ! ! *** ISORROPIA CODE ! *** SUBROUTINE CALCA1 ! *** CASE A1 ! ! THE MAIN CHARACTERISTICS OF THIS REGIME ARE: ! 1. SULFATE POOR (SULRAT > 2.0) ! 2. SOLID AEROSOL ONLY ! 3. SOLIDS POSSIBLE : (NH4)2SO4 ! ! A SIMPLE MATERIAL BALANCE IS PERFORMED, AND THE SOLID (NH4)2SO4 ! IS CALCULATED FROM THE SULFATES. THE EXCESS AMMONIA REMAINS IN ! THE GAS PHASE. ! ! *** COPYRIGHT 1996-2006, UNIVERSITY OF MIAMI, CARNEGIE MELLON UNIVERSITY, ! *** GEORGIA INSTITUTE OF TECHNOLOGY ! *** WRITTEN BY ATHANASIOS NENES ! *** UPDATED BY CHRISTOS FOUNTOUKIS ! !======================================================================= ! SUBROUTINE CALCA1 USE ISRPIA IMPLICIT DOUBLE PRECISION (A-H, O-Z) ! INCLUDE 'ISRPIA.INC' ! CNH42S4 = W(2) GNH3 = MAX (W(3)-2.0*CNH42S4, ZERO) RETURN ! ! *** END OF SUBROUTINE CALCA1 ****************************************** ! END SUBROUTINE CALCA1 !======================================================================= ! ! *** ISORROPIA CODE ! *** SUBROUTINE CALCB4 ! *** CASE B4 ! ! THE MAIN CHARACTERISTICS OF THIS REGIME ARE: ! 1. SULFATE RICH, NO FREE ACID (1.0 <= SULRAT < 2.0) ! 2. LIQUID AEROSOL PHASE ONLY POSSIBLE ! ! FOR CALCULATIONS, A BISECTION IS PERFORMED WITH RESPECT TO H+. ! THE OBJECTIVE FUNCTION IS THE DIFFERENCE BETWEEN THE ESTIMATED H+ ! AND THAT CALCULATED FROM ELECTRONEUTRALITY. ! ! *** COPYRIGHT 1996-2006, UNIVERSITY OF MIAMI, CARNEGIE MELLON UNIVERSITY, ! *** GEORGIA INSTITUTE OF TECHNOLOGY ! *** WRITTEN BY ATHANASIOS NENES ! *** UPDATED BY CHRISTOS FOUNTOUKIS ! !======================================================================= ! SUBROUTINE CALCB4 USE ISRPIA IMPLICIT DOUBLE PRECISION (A-H, O-Z) ! INCLUDE 'ISRPIA.INC' ! ! *** SOLVE EQUATIONS ************************************************** ! FRST = .TRUE. CALAIN = .TRUE. CALAOU = .TRUE. ! ! *** CALCULATE WATER CONTENT ****************************************** ! CALL CALCB1A ! GET DRY SALT CONTENT, AND USE FOR WATER. MOLALR(13) = CLC MOLALR(9) = CNH4HS4 MOLALR(4) = CNH42S4 CLC = ZERO CNH4HS4 = ZERO CNH42S4 = ZERO WATER = MOLALR(13)/M0(13)+MOLALR(9)/M0(9)+MOLALR(4)/M0(4) ! MOLAL(3) = W(3) ! NH4I ! DO 20 I=1,NSWEEP AK1 = XK1*((GAMA(8)/GAMA(7))**2.)*(WATER/GAMA(7)) BET = W(2) GAM = MOLAL(3) ! BB = BET + AK1 - GAM CC =-AK1*BET DD = BB*BB - 4.D0*CC ! ! *** SPECIATION & WATER CONTENT *************************************** ! MOLAL (5) = MAX(TINY,MIN(0.5*(-BB + SQRT(DD)), W(2))) ! SO4I MOLAL (6) = MAX(TINY,MIN(W(2)-MOLAL(5),W(2))) ! HSO4I MOLAL (1) = MAX(TINY,MIN(AK1*MOLAL(6)/MOLAL(5),W(2))) ! HI CALL CALCMR ! Water content ! ! *** CALCULATE ACTIVITIES OR TERMINATE INTERNAL LOOP ***************** ! IF (.NOT.CALAIN) GOTO 30 CALL CALCACT 20 CONTINUE ! 30 RETURN ! ! *** END OF SUBROUTINE CALCB4 ****************************************** ! END SUBROUTINE CALCB4 !======================================================================= ! ! *** ISORROPIA CODE ! *** SUBROUTINE CALCB3 ! *** CASE B3 ! ! THE MAIN CHARACTERISTICS OF THIS REGIME ARE: ! 1. SULFATE RICH, NO FREE ACID (1.0 <= SULRAT < 2.0) ! 2. BOTH LIQUID & SOLID PHASE IS POSSIBLE ! 3. SOLIDS POSSIBLE: (NH4)2SO4 ! ! *** COPYRIGHT 1996-2006, UNIVERSITY OF MIAMI, CARNEGIE MELLON UNIVERSITY, ! *** GEORGIA INSTITUTE OF TECHNOLOGY ! *** WRITTEN BY ATHANASIOS NENES ! *** UPDATED BY CHRISTOS FOUNTOUKIS ! !======================================================================= ! SUBROUTINE CALCB3 USE ISRPIA IMPLICIT DOUBLE PRECISION (A-H, O-Z) ! INCLUDE 'ISRPIA.INC' ! ! *** CALCULATE EQUIVALENT AMOUNT OF HSO4 AND SO4 *********************** ! X = MAX(2*W(2)-W(3), ZERO) ! Equivalent NH4HSO4 Y = MAX(W(3) -W(2), ZERO) ! Equivalent NH42SO4 ! ! *** CALCULATE SPECIES ACCORDING TO RELATIVE ABUNDANCE OF HSO4 ********* ! IF (X.LT.Y) THEN ! LC is the MIN (x,y) SCASE = 'B3 ; SUBCASE 1' TLC = X TNH42S4 = Y-X CALL CALCB3A (TLC,TNH42S4) ! LC + (NH4)2SO4 ELSE SCASE = 'B3 ; SUBCASE 2' TLC = Y TNH4HS4 = X-Y CALL CALCB3B (TLC,TNH4HS4) ! LC + NH4HSO4 ENDIF ! RETURN ! ! *** END OF SUBROUTINE CALCB3 ****************************************** ! END SUBROUTINE CALCB3 !======================================================================= ! ! *** ISORROPIA CODE ! *** SUBROUTINE CALCB3A ! *** CASE B3 ; SUBCASE 1 ! ! THE MAIN CHARACTERISTICS OF THIS REGIME ARE: ! 1. SULFATE RICH (1.0 < SULRAT < 2.0) ! 2. BOTH LIQUID & SOLID PHASE IS POSSIBLE ! 3. SOLIDS POSSIBLE: (NH4)2SO4 ! ! FOR CALCULATIONS, A BISECTION IS PERFORMED TOWARDS ZETA, THE ! AMOUNT OF SOLID (NH4)2SO4 DISSOLVED IN THE LIQUID PHASE. ! FOR EACH ESTIMATION OF ZETA, FUNCTION FUNCB3A CALCULATES THE ! AMOUNT OF H+ PRODUCED (BASED ON THE SO4 RELEASED INTO THE ! SOLUTION). THE SOLUBILITY PRODUCT OF (NH4)2SO4 IS USED AS THE ! OBJECTIVE FUNCTION. ! ! *** COPYRIGHT 1996-2006, UNIVERSITY OF MIAMI, CARNEGIE MELLON UNIVERSITY, ! *** GEORGIA INSTITUTE OF TECHNOLOGY ! *** WRITTEN BY ATHANASIOS NENES ! *** UPDATED BY CHRISTOS FOUNTOUKIS ! !======================================================================= ! SUBROUTINE CALCB3A (TLC, TNH42S4) USE ISRPIA IMPLICIT DOUBLE PRECISION (A-H, O-Z) ! INCLUDE 'ISRPIA.INC' ! CALAOU = .TRUE. ! Outer loop activity calculation flag ZLO = ZERO ! MIN DISSOLVED (NH4)2SO4 ZHI = TNH42S4 ! MAX DISSOLVED (NH4)2SO4 ! ! *** INITIAL VALUES FOR BISECTION (DISSOLVED (NH4)2SO4 **************** ! Z1 = ZLO Y1 = FUNCB3A (Z1, TLC, TNH42S4) IF (ABS(Y1).LE.EPS) RETURN YLO= Y1 ! ! *** ROOT TRACKING ; FOR THE RANGE OF HI AND LO *********************** ! DZ = (ZHI-ZLO)/REAL(NDIV) DO 10 I=1,NDIV Z2 = Z1+DZ Y2 = FUNCB3A (Z2, TLC, TNH42S4) IF (SIGN(1.d0,Y1)*SIGN(1.d0,Y2).LT.ZERO) GOTO 20 ! (Y1*Y2.LT.ZERO) Z1 = Z2 Y1 = Y2 10 CONTINUE ! ! *** NO SUBDIVISION WITH SOLUTION FOUND ! YHI= Y1 ! Save Y-value at HI position IF (ABS(Y2) .LT. EPS) THEN ! X2 IS A SOLUTION RETURN ! ! *** { YLO, YHI } < 0.0 THE SOLUTION IS ALWAYS UNDERSATURATED WITH LC ! ELSE IF (YLO.LT.ZERO .AND. YHI.LT.ZERO) THEN Z1 = ZHI Z2 = ZHI GOTO 40 ! ! *** { YLO, YHI } > 0.0 THE SOLUTION IS ALWAYS SUPERSATURATED WITH LC ! ELSE IF (YLO.GT.ZERO .AND. YHI.GT.ZERO) THEN Z1 = ZLO Z2 = ZLO GOTO 40 ELSE CALL PUSHERR (0001, 'CALCB3A') ! WARNING ERROR: NO SOLUTION RETURN ENDIF ! ! *** PERFORM BISECTION *********************************************** ! 20 DO 30 I=1,MAXIT Z3 = 0.5*(Z1+Z2) Y3 = FUNCB3A (Z3, TLC, TNH42S4) IF (SIGN(1.d0,Y1)*SIGN(1.d0,Y3) .LE. ZERO) THEN ! (Y1*Y3 .LE. ZERO) Y2 = Y3 Z2 = Z3 ELSE Y1 = Y3 Z1 = Z3 ENDIF IF (ABS(Z2-Z1) .LE. EPS*Z1) GOTO 40 30 CONTINUE CALL PUSHERR (0002, 'CALCB3A') ! WARNING ERROR: NO CONVERGENCE ! ! *** CONVERGED ; RETURN ************************************************ ! 40 ZK = 0.5*(Z1+Z2) Y3 = FUNCB3A (ZK, TLC, TNH42S4) ! RETURN ! ! *** END OF SUBROUTINE CALCB3A ****************************************** ! END SUBROUTINE CALCB3A !======================================================================= ! ! *** ISORROPIA CODE ! *** FUNCTION FUNCB3A ! *** CASE B3 ; SUBCASE 1 ! FUNCTION THAT SOLVES THE SYSTEM OF EQUATIONS FOR CASE B3 ! AND RETURNS THE VALUE OF THE ZEROED FUNCTION IN FUNCA3. ! !======================================================================= ! DOUBLE PRECISION FUNCTION FUNCB3A (ZK, Y, X) USE ISRPIA IMPLICIT DOUBLE PRECISION (A-H, O-Z) ! INCLUDE 'ISRPIA.INC' DOUBLE PRECISION KK ! ! *** SOLVE EQUATIONS ; WITH ITERATIONS FOR ACTIVITY COEF. ************ ! FRST = .TRUE. CALAIN = .TRUE. DO 20 I=1,NSWEEP GRAT1 = XK1*WATER/GAMA(7)*(GAMA(8)/GAMA(7))**2. DD = SQRT( (ZK+GRAT1+Y)**2. + 4.0*Y*GRAT1) KK = 0.5*(-(ZK+GRAT1+Y) + DD ) ! ! *** SPECIATION & WATER CONTENT *************************************** ! MOLAL (1) = KK ! HI MOLAL (5) = KK+ZK+Y ! SO4I MOLAL (6) = MAX (Y-KK, TINY) ! HSO4I MOLAL (3) = 3.0*Y+2*ZK ! NH4I CNH42S4 = X-ZK ! Solid (NH4)2SO4 CALL CALCMR ! Water content ! ! *** CALCULATE ACTIVITIES OR TERMINATE INTERNAL LOOP ***************** ! IF (FRST.AND.CALAOU .OR. .NOT.FRST.AND.CALAIN) THEN CALL CALCACT ELSE GOTO 30 ENDIF 20 CONTINUE ! ! *** CALCULATE OBJECTIVE FUNCTION ************************************ ! !CC30 FUNCB3A= ( SO4I*NH4I**2.0 )/( XK7*(WATER/GAMA(4))**3.0 ) 30 FUNCB3A= MOLAL(5)*MOLAL(3)**2.0 FUNCB3A= FUNCB3A/(XK7*(WATER/GAMA(4))**3.0) - ONE RETURN ! ! *** END OF FUNCTION FUNCB3A ******************************************** ! END FUNCTION FUNCB3A !======================================================================= ! ! *** ISORROPIA CODE ! *** SUBROUTINE CALCB3B ! *** CASE B3 ; SUBCASE 2 ! ! THE MAIN CHARACTERISTICS OF THIS REGIME ARE: ! 1. SULFATE RICH (1.0 < SULRAT < 2.0) ! 2. LIQUID PHASE ONLY IS POSSIBLE ! ! SPECIATION CALCULATIONS IS BASED ON THE HSO4 <--> SO4 EQUILIBRIUM. ! ! *** COPYRIGHT 1996-2006, UNIVERSITY OF MIAMI, CARNEGIE MELLON UNIVERSITY, ! *** GEORGIA INSTITUTE OF TECHNOLOGY ! *** WRITTEN BY ATHANASIOS NENES ! *** UPDATED BY CHRISTOS FOUNTOUKIS ! !======================================================================= ! SUBROUTINE CALCB3B (Y, X) USE ISRPIA IMPLICIT DOUBLE PRECISION (A-H, O-Z) ! INCLUDE 'ISRPIA.INC' DOUBLE PRECISION KK ! CALAOU = .FALSE. ! Outer loop activity calculation flag FRST = .FALSE. CALAIN = .TRUE. ! ! *** SOLVE EQUATIONS ; WITH ITERATIONS FOR ACTIVITY COEF. ************ ! DO 20 I=1,NSWEEP GRAT1 = XK1*WATER/GAMA(7)*(GAMA(8)/GAMA(7))**2. DD = SQRT( (GRAT1+Y)**2. + 4.0*(X+Y)*GRAT1) KK = 0.5*(-(GRAT1+Y) + DD ) ! ! *** SPECIATION & WATER CONTENT *************************************** ! MOLAL (1) = KK ! HI MOLAL (5) = Y+KK ! SO4I MOLAL (6) = MAX (X+Y-KK, TINY) ! HSO4I MOLAL (3) = 3.0*Y+X ! NH4I CALL CALCMR ! Water content ! ! *** CALCULATE ACTIVITIES OR TERMINATE INTERNAL LOOP ***************** ! IF (.NOT.CALAIN) GOTO 30 CALL CALCACT 20 CONTINUE ! 30 RETURN ! ! *** END OF SUBROUTINE CALCB3B ****************************************** ! END SUBROUTINE CALCB3B !======================================================================= ! ! *** ISORROPIA CODE ! *** SUBROUTINE CALCB2 ! *** CASE B2 ! ! THE MAIN CHARACTERISTICS OF THIS REGIME ARE: ! 1. SULFATE RICH, NO FREE ACID (1.0 <= SULRAT < 2.0) ! 2. THERE IS BOTH A LIQUID & SOLID PHASE ! 3. SOLIDS POSSIBLE : LC, (NH4)2SO4 ! ! THERE ARE TWO POSSIBLE REGIMES HERE, DEPENDING ON THE SULFATE RATIO: ! 1. WHEN BOTH LC AND (NH4)2SO4 ARE POSSIBLE (SUBROUTINE CALCB2A) ! 2. WHEN ONLY LC IS POSSIBLE (SUBROUTINE CALCB2B). ! ! *** COPYRIGHT 1996-2006, UNIVERSITY OF MIAMI, CARNEGIE MELLON UNIVERSITY, ! *** GEORGIA INSTITUTE OF TECHNOLOGY ! *** WRITTEN BY ATHANASIOS NENES ! *** UPDATED BY CHRISTOS FOUNTOUKIS ! !======================================================================= ! SUBROUTINE CALCB2 USE ISRPIA IMPLICIT DOUBLE PRECISION (A-H, O-Z) ! INCLUDE 'ISRPIA.INC' ! ! *** CALCULATE EQUIVALENT AMOUNT OF HSO4 AND SO4 *********************** ! X = MAX(2*W(2)-W(3), TINY) ! Equivalent NH4HSO4 Y = MAX(W(3) -W(2), TINY) ! Equivalent NH42SO4 ! ! *** CALCULATE SPECIES ACCORDING TO RELATIVE ABUNDANCE OF HSO4 ********* ! IF (X.LE.Y) THEN ! LC is the MIN (x,y) SCASE = 'B2 ; SUBCASE 1' CALL CALCB2A (X,Y-X) ! LC + (NH4)2SO4 POSSIBLE ELSE SCASE = 'B2 ; SUBCASE 2' CALL CALCB2B (Y,X-Y) ! LC ONLY POSSIBLE ENDIF ! RETURN ! ! *** END OF SUBROUTINE CALCB2 ****************************************** ! END SUBROUTINE CALCB2 !======================================================================= ! ! *** ISORROPIA CODE ! *** SUBROUTINE CALCB2 ! *** CASE B2 ; SUBCASE A. ! ! THE MAIN CHARACTERISTICS OF THIS REGIME ARE: ! 1. SULFATE RICH (1.0 < SULRAT < 2.0) ! 2. SOLID PHASE ONLY POSSIBLE ! 3. SOLIDS POSSIBLE: LC, (NH4)2SO4 ! ! THERE ARE TWO POSSIBLE REGIMES HERE, DEPENDING ON RELATIVE HUMIDITY: ! 1. WHEN RH >= MDRH ; LIQUID PHASE POSSIBLE (MDRH REGION) ! 2. WHEN RH < MDRH ; ONLY SOLID PHASE POSSIBLE ! ! FOR SOLID CALCULATIONS, A MATERIAL BALANCE BASED ON THE STOICHIMETRIC ! PROPORTION OF AMMONIUM AND SULFATE IS DONE TO CALCULATE THE AMOUNT ! OF LC AND (NH4)2SO4 IN THE SOLID PHASE. ! ! *** COPYRIGHT 1996-2006, UNIVERSITY OF MIAMI, CARNEGIE MELLON UNIVERSITY, ! *** GEORGIA INSTITUTE OF TECHNOLOGY ! *** WRITTEN BY ATHANASIOS NENES ! *** UPDATED BY CHRISTOS FOUNTOUKIS ! !======================================================================= ! SUBROUTINE CALCB2A (TLC, TNH42S4) USE ISRPIA IMPLICIT DOUBLE PRECISION (A-H, O-Z) ! INCLUDE 'ISRPIA.INC' ! ! *** REGIME DEPENDS UPON THE AMBIENT RELATIVE HUMIDITY ***************** ! IF (RH.LT.DRMLCAS) THEN SCASE = 'B2 ; SUBCASE A1' ! SOLIDS POSSIBLE ONLY CLC = TLC CNH42S4 = TNH42S4 SCASE = 'B2 ; SUBCASE A1' ELSE SCASE = 'B2 ; SUBCASE A2' CALL CALCB2A2 (TLC, TNH42S4) ! LIQUID & SOLID PHASE POSSIBLE SCASE = 'B2 ; SUBCASE A2' ENDIF ! RETURN ! ! *** END OF SUBROUTINE CALCB2A ***************************************** ! END SUBROUTINE CALCB2A !======================================================================= ! ! *** ISORROPIA CODE ! *** SUBROUTINE CALCB2A2 ! *** CASE B2 ; SUBCASE A2. ! ! THE MAIN CHARACTERISTICS OF THIS REGIME ARE: ! 1. SULFATE RICH (1.0 < SULRAT < 2.0) ! 2. SOLID PHASE ONLY POSSIBLE ! 3. SOLIDS POSSIBLE: LC, (NH4)2SO4 ! ! THIS IS THE CASE WHERE THE RELATIVE HUMIDITY IS IN THE MUTUAL ! DRH REGION. THE SOLUTION IS ASSUMED TO BE THE SUM OF TWO WEIGHTED ! SOLUTIONS ; THE SOLID PHASE ONLY (SUBROUTINE CALCB2A1) AND THE ! THE SOLID WITH LIQUID PHASE (SUBROUTINE CALCB3). ! ! *** COPYRIGHT 1996-2006, UNIVERSITY OF MIAMI, CARNEGIE MELLON UNIVERSITY, ! *** GEORGIA INSTITUTE OF TECHNOLOGY ! *** WRITTEN BY ATHANASIOS NENES ! *** UPDATED BY CHRISTOS FOUNTOUKIS ! !======================================================================= ! SUBROUTINE CALCB2A2 (TLC, TNH42S4) USE ISRPIA IMPLICIT DOUBLE PRECISION (A-H, O-Z) ! INCLUDE 'ISRPIA.INC' ! ! *** FIND WEIGHT FACTOR ********************************************** ! IF (WFTYP.EQ.0) THEN WF = ZERO ELSEIF (WFTYP.EQ.1) THEN WF = 0.5D0 ELSE WF = (DRLC-RH)/(DRLC-DRMLCAS) ENDIF ONEMWF = ONE - WF ! ! *** FIND FIRST SECTION ; DRY ONE ************************************ ! CLCO = TLC ! FIRST (DRY) SOLUTION CNH42SO = TNH42S4 ! ! *** FIND SECOND SECTION ; DRY & LIQUID ****************************** ! CLC = ZERO CNH42S4 = ZERO CALL CALCB3 ! SECOND (LIQUID) SOLUTION ! ! *** FIND SOLUTION AT MDRH BY WEIGHTING DRY & LIQUID SOLUTIONS. ! MOLAL(1)= ONEMWF*MOLAL(1) ! H+ MOLAL(3)= ONEMWF*(2.D0*(CNH42SO-CNH42S4) + 3.D0*(CLCO-CLC)) ! NH4+ MOLAL(5)= ONEMWF*(CNH42SO-CNH42S4 + CLCO-CLC) ! SO4-- MOLAL(6)= ONEMWF*(CLCO-CLC) ! HSO4- ! WATER = ONEMWF*WATER ! CLC = WF*CLCO + ONEMWF*CLC CNH42S4 = WF*CNH42SO + ONEMWF*CNH42S4 ! RETURN ! ! *** END OF SUBROUTINE CALCB2A2 **************************************** ! END SUBROUTINE CALCB2A2 !======================================================================= ! ! *** ISORROPIA CODE ! *** SUBROUTINE CALCB2 ! *** CASE B2 ; SUBCASE B ! ! THE MAIN CHARACTERISTICS OF THIS REGIME ARE: ! 1. SULFATE RICH (1.0 < SULRAT < 2.0) ! 2. BOTH LIQUID & SOLID PHASE IS POSSIBLE ! 3. SOLIDS POSSIBLE: LC ! ! FOR CALCULATIONS, A BISECTION IS PERFORMED TOWARDS ZETA, THE ! AMOUNT OF SOLID LC DISSOLVED IN THE LIQUID PHASE. ! FOR EACH ESTIMATION OF ZETA, FUNCTION FUNCB2A CALCULATES THE ! AMOUNT OF H+ PRODUCED (BASED ON THE HSO4, SO4 RELEASED INTO THE ! SOLUTION). THE SOLUBILITY PRODUCT OF LC IS USED AS THE OBJECTIVE ! FUNCTION. ! ! *** COPYRIGHT 1996-2006, UNIVERSITY OF MIAMI, CARNEGIE MELLON UNIVERSITY, ! *** GEORGIA INSTITUTE OF TECHNOLOGY ! *** WRITTEN BY ATHANASIOS NENES ! *** UPDATED BY CHRISTOS FOUNTOUKIS ! !======================================================================= ! SUBROUTINE CALCB2B (TLC,TNH4HS4) USE ISRPIA IMPLICIT DOUBLE PRECISION (A-H, O-Z) ! INCLUDE 'ISRPIA.INC' ! CALAOU = .TRUE. ! Outer loop activity calculation flag ZLO = ZERO ZHI = TLC ! High limit: all of it in liquid phase ! ! *** INITIAL VALUES FOR BISECTION ************************************** ! X1 = ZHI Y1 = FUNCB2B (X1,TNH4HS4,TLC) IF (ABS(Y1).LE.EPS) RETURN YHI= Y1 ! Save Y-value at Hi position ! ! *** ROOT TRACKING ; FOR THE RANGE OF HI AND LO ************************ ! DX = (ZHI-ZLO)/NDIV DO 10 I=1,NDIV X2 = X1-DX Y2 = FUNCB2B (X2,TNH4HS4,TLC) IF (SIGN(1.d0,Y1)*SIGN(1.d0,Y2).LT.ZERO) GOTO 20 ! (Y1*Y2.LT.ZERO) X1 = X2 Y1 = Y2 10 CONTINUE ! ! *** NO SUBDIVISION WITH SOLUTION FOUND ! YLO= Y1 ! Save Y-value at LO position IF (ABS(Y2) .LT. EPS) THEN ! X2 IS A SOLUTION RETURN ! ! *** { YLO, YHI } < 0.0 THE SOLUTION IS ALWAYS UNDERSATURATED WITH LC ! ELSE IF (YLO.LT.ZERO .AND. YHI.LT.ZERO) THEN X1 = ZHI X2 = ZHI GOTO 40 ! ! *** { YLO, YHI } > 0.0 THE SOLUTION IS ALWAYS SUPERSATURATED WITH LC ! ELSE IF (YLO.GT.ZERO .AND. YHI.GT.ZERO) THEN X1 = ZLO X2 = ZLO GOTO 40 ELSE CALL PUSHERR (0001, 'CALCB2B') ! WARNING ERROR: NO SOLUTION RETURN ENDIF ! ! *** PERFORM BISECTION ************************************************* ! 20 DO 30 I=1,MAXIT X3 = 0.5*(X1+X2) Y3 = FUNCB2B (X3,TNH4HS4,TLC) IF (SIGN(1.d0,Y1)*SIGN(1.d0,Y3) .LE. ZERO) THEN ! (Y1*Y3 .LE. ZERO) Y2 = Y3 X2 = X3 ELSE Y1 = Y3 X1 = X3 ENDIF IF (ABS(X2-X1) .LE. EPS*X1) GOTO 40 30 CONTINUE CALL PUSHERR (0002, 'CALCB2B') ! WARNING ERROR: NO CONVERGENCE ! ! *** CONVERGED ; RETURN ************************************************ ! 40 X3 = 0.5*(X1+X2) Y3 = FUNCB2B (X3,TNH4HS4,TLC) ! RETURN ! ! *** END OF SUBROUTINE CALCB2B ***************************************** ! END SUBROUTINE CALCB2B !======================================================================= ! ! *** ISORROPIA CODE ! *** FUNCTION FUNCB2B ! *** CASE B2 ; ! FUNCTION THAT SOLVES THE SYSTEM OF EQUATIONS FOR CASE B2 ; SUBCASE 2 ! AND RETURNS THE VALUE OF THE ZEROED FUNCTION IN FUNCB2B. ! !======================================================================= ! DOUBLE PRECISION FUNCTION FUNCB2B (X,TNH4HS4,TLC) USE ISRPIA IMPLICIT DOUBLE PRECISION (A-H, O-Z) ! INCLUDE 'ISRPIA.INC' ! ! *** SOLVE EQUATIONS ************************************************** ! FRST = .TRUE. CALAIN = .TRUE. DO 20 I=1,NSWEEP GRAT2 = XK1*WATER*(GAMA(8)/GAMA(7))**2./GAMA(7) PARM = X+GRAT2 DELTA = PARM*PARM + 4.0*(X+TNH4HS4)*GRAT2 ! Diakrinousa OMEGA = 0.5*(-PARM + SQRT(DELTA)) ! Thetiki riza (ie:H+>0) ! ! *** SPECIATION & WATER CONTENT *************************************** ! MOLAL (1) = OMEGA ! HI MOLAL (3) = 3.0*X+TNH4HS4 ! NH4I MOLAL (5) = X+OMEGA ! SO4I MOLAL (6) = MAX (X+TNH4HS4-OMEGA, TINY) ! HSO4I CLC = MAX(TLC-X,ZERO) ! Solid LC CALL CALCMR ! Water content ! ! *** CALCULATE ACTIVITIES OR TERMINATE INTERNAL LOOP ****************** ! IF (FRST.AND.CALAOU .OR. .NOT.FRST.AND.CALAIN) THEN CALL CALCACT ELSE GOTO 30 ENDIF 20 CONTINUE ! ! *** CALCULATE OBJECTIVE FUNCTION ************************************** ! !CC30 FUNCB2B= ( NH4I**3.*SO4I*HSO4I )/( XK13*(WATER/GAMA(13))**5. ) 30 FUNCB2B= (MOLAL(3)**3.)*MOLAL(5)*MOLAL(6) FUNCB2B= FUNCB2B/(XK13*(WATER/GAMA(13))**5.) - ONE RETURN ! ! *** END OF FUNCTION FUNCB2B ******************************************* ! END FUNCTION FUNCB2B !======================================================================= ! ! *** ISORROPIA CODE ! *** SUBROUTINE CALCB1 ! *** CASE B1 ! ! THE MAIN CHARACTERISTICS OF THIS REGIME ARE: ! 1. SULFATE RICH, NO FREE ACID (1.0 <= SULRAT < 2.0) ! 2. THERE IS BOTH A LIQUID & SOLID PHASE ! 3. SOLIDS POSSIBLE : LC, (NH4)2SO4, NH4HSO4 ! ! THERE ARE TWO POSSIBLE REGIMES HERE, DEPENDING ON RELATIVE HUMIDITY: ! 1. WHEN RH >= MDRH ; LIQUID PHASE POSSIBLE (MDRH REGION) ! 2. WHEN RH < MDRH ; ONLY SOLID PHASE POSSIBLE (SUBROUTINE CALCB1A) ! ! *** COPYRIGHT 1996-2006, UNIVERSITY OF MIAMI, CARNEGIE MELLON UNIVERSITY, ! *** GEORGIA INSTITUTE OF TECHNOLOGY ! *** WRITTEN BY ATHANASIOS NENES ! *** UPDATED BY CHRISTOS FOUNTOUKIS ! !======================================================================= ! SUBROUTINE CALCB1 USE ISRPIA IMPLICIT DOUBLE PRECISION (A-H, O-Z) ! INCLUDE 'ISRPIA.INC' ! ! *** REGIME DEPENDS UPON THE AMBIENT RELATIVE HUMIDITY ***************** ! IF (RH.LT.DRMLCAB) THEN SCASE = 'B1 ; SUBCASE 1' CALL CALCB1A ! SOLID PHASE ONLY POSSIBLE SCASE = 'B1 ; SUBCASE 1' ELSE SCASE = 'B1 ; SUBCASE 2' CALL CALCB1B ! LIQUID & SOLID PHASE POSSIBLE SCASE = 'B1 ; SUBCASE 2' ENDIF ! RETURN ! ! *** END OF SUBROUTINE CALCB1 ****************************************** ! END SUBROUTINE CALCB1 !======================================================================= ! ! *** ISORROPIA CODE ! *** SUBROUTINE CALCB1A ! *** CASE B1 ; SUBCASE 1 ! ! THE MAIN CHARACTERISTICS OF THIS REGIME ARE: ! 1. SULFATE RICH ! 2. THERE IS NO LIQUID PHASE ! 3. SOLIDS POSSIBLE: LC, { (NH4)2SO4 XOR NH4HSO4 } (ONE OF TWO ! BUT NOT BOTH) ! ! A SIMPLE MATERIAL BALANCE IS PERFORMED, AND THE AMOUNT OF LC ! IS CALCULATED FROM THE (NH4)2SO4 AND NH4HSO4 WHICH IS LEAST ! ABUNDANT (STOICHIMETRICALLY). THE REMAINING EXCESS OF SALT ! IS MIXED WITH THE LC. ! ! *** COPYRIGHT 1996-2006, UNIVERSITY OF MIAMI, CARNEGIE MELLON UNIVERSITY, ! *** GEORGIA INSTITUTE OF TECHNOLOGY ! *** WRITTEN BY ATHANASIOS NENES ! *** UPDATED BY CHRISTOS FOUNTOUKIS ! !======================================================================= ! SUBROUTINE CALCB1A USE ISRPIA IMPLICIT DOUBLE PRECISION (A-H, O-Z) ! INCLUDE 'ISRPIA.INC' ! ! *** SETUP PARAMETERS ************************************************ ! X = 2*W(2)-W(3) ! Equivalent NH4HSO4 Y = W(3)-W(2) ! Equivalent (NH4)2SO4 ! ! *** CALCULATE COMPOSITION ******************************************* ! IF (X.LE.Y) THEN ! LC is the MIN (x,y) CLC = X ! NH4HSO4 >= (NH4)2S04 CNH4HS4 = ZERO CNH42S4 = Y-X ELSE CLC = Y ! NH4HSO4 < (NH4)2S04 CNH4HS4 = X-Y CNH42S4 = ZERO ENDIF RETURN ! ! *** END OF SUBROUTINE CALCB1 ****************************************** ! END SUBROUTINE CALCB1A !======================================================================= ! ! *** ISORROPIA CODE ! *** SUBROUTINE CALCB1B ! *** CASE B1 ; SUBCASE 2 ! ! THE MAIN CHARACTERISTICS OF THIS REGIME ARE: ! 1. SULFATE RICH ! 2. THERE IS BOTH A LIQUID & SOLID PHASE ! 3. SOLIDS POSSIBLE: LC, { (NH4)2SO4 XOR NH4HSO4 } (ONE OF TWO ! BUT NOT BOTH) ! ! THIS IS THE CASE WHERE THE RELATIVE HUMIDITY IS IN THE MUTUAL ! DRH REGION. THE SOLUTION IS ASSUMED TO BE THE SUM OF TWO WEIGHTED ! SOLUTIONS ; THE SOLID PHASE ONLY (SUBROUTINE CALCB1A) AND THE ! THE SOLID WITH LIQUID PHASE (SUBROUTINE CALCB2). ! ! *** COPYRIGHT 1996-2006, UNIVERSITY OF MIAMI, CARNEGIE MELLON UNIVERSITY, ! *** GEORGIA INSTITUTE OF TECHNOLOGY ! *** WRITTEN BY ATHANASIOS NENES ! *** UPDATED BY CHRISTOS FOUNTOUKIS ! !======================================================================= ! SUBROUTINE CALCB1B USE ISRPIA IMPLICIT DOUBLE PRECISION (A-H, O-Z) ! INCLUDE 'ISRPIA.INC' ! ! *** FIND WEIGHT FACTOR ********************************************** ! IF (WFTYP.EQ.0) THEN WF = ZERO ELSEIF (WFTYP.EQ.1) THEN WF = 0.5D0 ELSE WF = (DRNH4HS4-RH)/(DRNH4HS4-DRMLCAB) ENDIF ONEMWF = ONE - WF ! ! *** FIND FIRST SECTION ; DRY ONE ************************************ ! CALL CALCB1A CLCO = CLC ! FIRST (DRY) SOLUTION CNH42SO = CNH42S4 CNH4HSO = CNH4HS4 ! ! *** FIND SECOND SECTION ; DRY & LIQUID ****************************** ! CLC = ZERO CNH42S4 = ZERO CNH4HS4 = ZERO CALL CALCB2 ! SECOND (LIQUID) SOLUTION ! ! *** FIND SOLUTION AT MDRH BY WEIGHTING DRY & LIQUID SOLUTIONS. ! MOLAL(1)= ONEMWF*MOLAL(1) ! H+ MOLAL(3)= ONEMWF*(2.D0*(CNH42SO-CNH42S4) + (CNH4HSO-CNH4HS4) & + 3.D0*(CLCO-CLC)) ! NH4+ MOLAL(5)= ONEMWF*(CNH42SO-CNH42S4 + CLCO-CLC) ! SO4-- MOLAL(6)= ONEMWF*(CNH4HSO-CNH4HS4 + CLCO-CLC) ! HSO4- ! WATER = ONEMWF*WATER ! CLC = WF*CLCO + ONEMWF*CLC CNH42S4 = WF*CNH42SO + ONEMWF*CNH42S4 CNH4HS4 = WF*CNH4HSO + ONEMWF*CNH4HS4 ! RETURN ! ! *** END OF SUBROUTINE CALCB1B ***************************************** ! END SUBROUTINE CALCB1B !======================================================================= ! ! *** ISORROPIA CODE ! *** SUBROUTINE CALCC2 ! *** CASE C2 ! ! THE MAIN CHARACTERISTICS OF THIS REGIME ARE: ! 1. SULFATE RICH, FREE ACID (SULRAT < 1.0) ! 2. THERE IS ONLY A LIQUID PHASE ! ! *** COPYRIGHT 1996-2006, UNIVERSITY OF MIAMI, CARNEGIE MELLON UNIVERSITY, ! *** GEORGIA INSTITUTE OF TECHNOLOGY ! *** WRITTEN BY ATHANASIOS NENES ! *** UPDATED BY CHRISTOS FOUNTOUKIS ! !======================================================================= ! SUBROUTINE CALCC2 USE ISRPIA IMPLICIT DOUBLE PRECISION (A-H, O-Z) ! INCLUDE 'ISRPIA.INC' DOUBLE PRECISION LAMDA, KAPA ! CALAOU =.TRUE. ! Outer loop activity calculation flag FRST =.TRUE. CALAIN =.TRUE. ! ! *** SOLVE EQUATIONS ************************************************** ! LAMDA = W(3) ! NH4HSO4 INITIALLY IN SOLUTION PSI = W(2)-W(3) ! H2SO4 IN SOLUTION DO 20 I=1,NSWEEP PARM = WATER*XK1/GAMA(7)*(GAMA(8)/GAMA(7))**2. BB = PSI+PARM CC =-PARM*(LAMDA+PSI) KAPA = 0.5*(-BB+SQRT(BB*BB-4.0*CC)) ! ! *** SPECIATION & WATER CONTENT *************************************** ! MOLAL(1) = PSI+KAPA ! HI MOLAL(3) = LAMDA ! NH4I MOLAL(5) = KAPA ! SO4I MOLAL(6) = MAX(LAMDA+PSI-KAPA, TINY) ! HSO4I CH2SO4 = MAX(MOLAL(5)+MOLAL(6)-MOLAL(3), ZERO) ! Free H2SO4 CALL CALCMR ! Water content ! ! *** CALCULATE ACTIVITIES OR TERMINATE INTERNAL LOOP ***************** ! IF (.NOT.CALAIN) GOTO 30 CALL CALCACT 20 CONTINUE ! 30 RETURN ! ! *** END OF SUBROUTINE CALCC2 ***************************************** ! END SUBROUTINE CALCC2 !======================================================================= ! ! *** ISORROPIA CODE ! *** SUBROUTINE CALCC1 ! *** CASE C1 ! ! THE MAIN CHARACTERISTICS OF THIS REGIME ARE: ! 1. SULFATE RICH, FREE ACID (SULRAT < 1.0) ! 2. THERE IS BOTH A LIQUID & SOLID PHASE ! 3. SOLIDS POSSIBLE: NH4HSO4 ! ! *** COPYRIGHT 1996-2006, UNIVERSITY OF MIAMI, CARNEGIE MELLON UNIVERSITY, ! *** GEORGIA INSTITUTE OF TECHNOLOGY ! *** WRITTEN BY ATHANASIOS NENES ! *** UPDATED BY CHRISTOS FOUNTOUKIS ! !======================================================================= ! SUBROUTINE CALCC1 USE ISRPIA IMPLICIT DOUBLE PRECISION (A-H, O-Z) ! INCLUDE 'ISRPIA.INC' DOUBLE PRECISION KLO, KHI ! CALAOU = .TRUE. ! Outer loop activity calculation flag KLO = TINY KHI = W(3) ! ! *** INITIAL VALUES FOR BISECTION ************************************* ! X1 = KLO Y1 = FUNCC1 (X1) IF (ABS(Y1).LE.EPS) GOTO 50 YLO= Y1 ! ! *** ROOT TRACKING ; FOR THE RANGE OF HI AND LO *********************** ! DX = (KHI-KLO)/REAL(NDIV) DO 10 I=1,NDIV X2 = X1+DX Y2 = FUNCC1 (X2) IF (SIGN(1.d0,Y1)*SIGN(1.d0,Y2) .LT. ZERO) GOTO 20 ! (Y1*Y2 .LT. ZERO) X1 = X2 Y1 = Y2 10 CONTINUE ! ! *** NO SUBDIVISION WITH SOLUTION FOUND ! YHI= Y2 ! Save Y-value at HI position IF (ABS(Y2) .LT. EPS) THEN ! X2 IS A SOLUTION GOTO 50 ! ! *** { YLO, YHI } < 0.0 SOLUTION IS ALWAYS UNDERSATURATED WITH NH4HS04 ! ELSE IF (YLO.LT.ZERO .AND. YHI.LT.ZERO) THEN GOTO 50 ! ! *** { YLO, YHI } > 0.0 SOLUTION IS ALWAYS SUPERSATURATED WITH NH4HS04 ! ELSE IF (YLO.GT.ZERO .AND. YHI.GT.ZERO) THEN X1 = KLO X2 = KLO GOTO 40 ELSE CALL PUSHERR (0001, 'CALCC1') ! WARNING ERROR: NO SOLUTION GOTO 50 ENDIF ! ! *** PERFORM BISECTION OF DISSOLVED NH4HSO4 ************************** ! 20 DO 30 I=1,MAXIT X3 = 0.5*(X1+X2) Y3 = FUNCC1 (X3) IF (SIGN(1.d0,Y1)*SIGN(1.d0,Y3) .LE. ZERO) THEN ! (Y1*Y3 .LE. ZERO) Y2 = Y3 X2 = X3 ELSE Y1 = Y3 X1 = X3 ENDIF IF (ABS(X2-X1) .LE. EPS*X1) GOTO 40 30 CONTINUE CALL PUSHERR (0002, 'CALCC1') ! WARNING ERROR: NO CONVERGENCE ! ! *** CONVERGED ; RETURN *********************************************** ! 40 X3 = 0.5*(X1+X2) Y3 = FUNCC1 (X3) ! 50 RETURN ! ! *** END OF SUBROUTINE CALCC1 ***************************************** ! END SUBROUTINE CALCC1 !======================================================================= ! ! *** ISORROPIA CODE ! *** FUNCTION FUNCC1 ! *** CASE C1 ; ! FUNCTION THAT SOLVES THE SYSTEM OF EQUATIONS FOR CASE C1 ! AND RETURNS THE VALUE OF THE ZEROED FUNCTION IN FUNCC1. ! ! *** COPYRIGHT 1996-2006, UNIVERSITY OF MIAMI, CARNEGIE MELLON UNIVERSITY, ! *** GEORGIA INSTITUTE OF TECHNOLOGY ! *** WRITTEN BY ATHANASIOS NENES ! *** UPDATED BY CHRISTOS FOUNTOUKIS ! !======================================================================= ! DOUBLE PRECISION FUNCTION FUNCC1 (KAPA) USE ISRPIA IMPLICIT DOUBLE PRECISION (A-H, O-Z) ! INCLUDE 'ISRPIA.INC' DOUBLE PRECISION KAPA, LAMDA ! ! *** SOLVE EQUATIONS ************************************************** ! FRST = .TRUE. CALAIN = .TRUE. ! PSI = W(2)-W(3) DO 20 I=1,NSWEEP PAR1 = XK1*WATER/GAMA(7)*(GAMA(8)/GAMA(7))**2.0 PAR2 = XK12*(WATER/GAMA(9))**2.0 BB = PSI + PAR1 CC =-PAR1*(PSI+KAPA) LAMDA = 0.5*(-BB+SQRT(BB*BB-4*CC)) ! ! *** SAVE CONCENTRATIONS IN MOLAL ARRAY ******************************* ! MOLAL(1) = PSI+LAMDA ! HI MOLAL(3) = KAPA ! NH4I MOLAL(5) = LAMDA ! SO4I MOLAL(6) = MAX (ZERO, PSI+KAPA-LAMDA) ! HSO4I CNH4HS4 = MAX(W(3)-MOLAL(3), ZERO) ! Solid NH4HSO4 CH2SO4 = MAX(PSI, ZERO) ! Free H2SO4 CALL CALCMR ! Water content ! ! *** CALCULATE ACTIVITIES OR TERMINATE INTERNAL LOOP ***************** ! IF (FRST.AND.CALAOU .OR. .NOT.FRST.AND.CALAIN) THEN CALL CALCACT ELSE GOTO 30 ENDIF 20 CONTINUE ! ! *** CALCULATE ZERO FUNCTION ******************************************* ! !CC30 FUNCC1= (NH4I*HSO4I/PAR2) - ONE 30 FUNCC1= (MOLAL(3)*MOLAL(6)/PAR2) - ONE RETURN ! ! *** END OF FUNCTION FUNCC1 ******************************************** ! END FUNCTION FUNCC1 !======================================================================= ! ! *** ISORROPIA CODE ! *** SUBROUTINE CALCD3 ! *** CASE D3 ! ! THE MAIN CHARACTERISTICS OF THIS REGIME ARE: ! 1. SULFATE POOR (SULRAT > 2.0) ! 2. THERE IS OLNY A LIQUID PHASE ! ! *** COPYRIGHT 1996-2006, UNIVERSITY OF MIAMI, CARNEGIE MELLON UNIVERSITY, ! *** GEORGIA INSTITUTE OF TECHNOLOGY ! *** WRITTEN BY ATHANASIOS NENES ! *** UPDATED BY CHRISTOS FOUNTOUKIS ! !======================================================================= ! SUBROUTINE CALCD3 USE ISRPIA USE SOLUT IMPLICIT DOUBLE PRECISION (A-H, O-Z) ! INCLUDE 'ISRPIA.INC' ! COMMON /SOLUT/ CHI1, CHI2, CHI3, CHI4, CHI5, CHI6, CHI7, CHI8, & ! PSI1, PSI2, PSI3, PSI4, PSI5, PSI6, PSI7, PSI8, & ! A1, A2, A3, A4, A5, A6, A7, A8 ! ! *** FIND DRY COMPOSITION ********************************************** ! CALL CALCD1A ! ! *** SETUP PARAMETERS ************************************************ ! CHI1 = CNH4NO3 ! Save from CALCD1 run CHI2 = CNH42S4 CHI3 = GHNO3 CHI4 = GNH3 ! PSI1 = CNH4NO3 ! ASSIGN INITIAL PSI's PSI2 = CHI2 PSI3 = ZERO PSI4 = ZERO ! MOLAL(5) = ZERO MOLAL(6) = ZERO MOLAL(3) = PSI1 MOLAL(7) = PSI1 CALL CALCMR ! Initial water ! CALAOU = .TRUE. ! Outer loop activity calculation flag PSI4LO = TINY ! Low limit PSI4HI = CHI4 ! High limit ! ! *** INITIAL VALUES FOR BISECTION ************************************ ! 60 X1 = PSI4LO Y1 = FUNCD3 (X1) IF (ABS(Y1).LE.EPS) RETURN YLO= Y1 ! Save Y-value at HI position ! ! *** ROOT TRACKING ; FOR THE RANGE OF HI AND LO ********************** ! DX = (PSI4HI-PSI4LO)/REAL(NDIV) DO 10 I=1,NDIV X2 = X1+DX Y2 = FUNCD3 (X2) IF (SIGN(1.d0,Y1)*SIGN(1.d0,Y2).LT.ZERO) GOTO 20 ! (Y1*Y2.LT.ZERO) X1 = X2 Y1 = Y2 10 CONTINUE ! ! *** NO SUBDIVISION WITH SOLUTION FOUND ! YHI= Y1 ! Save Y-value at Hi position IF (ABS(Y2) .LT. EPS) THEN ! X2 IS A SOLUTION RETURN ! ! *** { YLO, YHI } < 0.0 THE SOLUTION IS ALWAYS UNDERSATURATED WITH NH3 ! Physically I dont know when this might happen, but I have put this ! branch in for completeness. I assume there is no solution; all NO3 goes to the ! gas phase. ! ELSE IF (YLO.LT.ZERO .AND. YHI.LT.ZERO) THEN P4 = TINY ! PSI4LO ! CHI4 YY = FUNCD3(P4) GOTO 50 ! ! *** { YLO, YHI } > 0.0 THE SOLUTION IS ALWAYS SUPERSATURATED WITH NH3 ! This happens when Sul.Rat. = 2.0, so some NH4+ from sulfate evaporates ! and goes to the gas phase ; so I redefine the LO and HI limits of PSI4 ! and proceed again with root tracking. ! ELSE IF (YLO.GT.ZERO .AND. YHI.GT.ZERO) THEN PSI4HI = PSI4LO PSI4LO = PSI4LO - 0.1*(PSI1+PSI2) ! No solution; some NH3 evaporates IF (PSI4LO.LT.-(PSI1+PSI2)) THEN CALL PUSHERR (0001, 'CALCD3') ! WARNING ERROR: NO SOLUTION RETURN ELSE MOLAL(5) = ZERO MOLAL(6) = ZERO MOLAL(3) = PSI1 MOLAL(7) = PSI1 CALL CALCMR ! Initial water GOTO 60 ! Redo root tracking ENDIF ENDIF ! ! *** PERFORM BISECTION *********************************************** ! 20 DO 30 I=1,MAXIT X3 = 0.5*(X1+X2) Y3 = FUNCD3 (X3) IF (SIGN(1.d0,Y1)*SIGN(1.d0,Y3) .LE. ZERO) THEN ! (Y1*Y3 .LE. ZERO) Y2 = Y3 X2 = X3 ELSE Y1 = Y3 X1 = X3 ENDIF IF (ABS(X2-X1) .LE. EPS*ABS(X1)) GOTO 40 30 CONTINUE CALL PUSHERR (0002, 'CALCD3') ! WARNING ERROR: NO CONVERGENCE ! ! *** CONVERGED ; RETURN ********************************************** ! 40 X3 = 0.5*(X1+X2) Y3 = FUNCD3 (X3) ! ! *** CALCULATE HSO4 SPECIATION AND RETURN ******************************* ! 50 CONTINUE IF (MOLAL(1).GT.TINY) THEN CALL CALCHS4 (MOLAL(1), MOLAL(5), ZERO, DELTA) MOLAL(1) = MOLAL(1) - DELTA ! H+ EFFECT MOLAL(5) = MOLAL(5) - DELTA ! SO4 EFFECT MOLAL(6) = DELTA ! HSO4 EFFECT ENDIF RETURN ! ! *** END OF SUBROUTINE CALCD3 ****************************************** ! END SUBROUTINE CALCD3 !======================================================================= ! ! *** ISORROPIA CODE ! *** FUNCTION FUNCD3 ! *** CASE D3 ! FUNCTION THAT SOLVES THE SYSTEM OF EQUATIONS FOR CASE D3 ; ! AND RETURNS THE VALUE OF THE ZEROED FUNCTION IN FUNCD3. ! !======================================================================= ! DOUBLE PRECISION FUNCTION FUNCD3 (P4) USE ISRPIA USE SOLUT IMPLICIT DOUBLE PRECISION (A-H, O-Z) ! INCLUDE 'ISRPIA.INC' ! COMMON /SOLUT/ CHI1, CHI2, CHI3, CHI4, CHI5, CHI6, CHI7, CHI8, & ! PSI1, PSI2, PSI3, PSI4, PSI5, PSI6, PSI7, PSI8, & ! A1, A2, A3, A4, A5, A6, A7, A8 ! ! *** SETUP PARAMETERS ************************************************ ! FRST = .TRUE. CALAIN = .TRUE. PSI4 = P4 ! ! *** SOLVE EQUATIONS ; WITH ITERATIONS FOR ACTIVITY COEF. ************ ! DO 10 I=1,NSWEEP A2 = XK7*(WATER/GAMA(4))**3.0 A3 = XK4*R*TEMP*(WATER/GAMA(10))**2.0 A4 = (XK2/XKW)*R*TEMP*(GAMA(10)/GAMA(5))**2.0 A7 = XKW *RH*WATER*WATER ! PSI3 = A3*A4*CHI3*(CHI4-PSI4) - PSI1*(2.D0*PSI2+PSI1+PSI4) PSI3 = PSI3/(A3*A4*(CHI4-PSI4) + 2.D0*PSI2+PSI1+PSI4) PSI3 = MIN(MAX(PSI3, ZERO), CHI3) ! BB = PSI4 - PSI3 !CCOLD AHI = 0.5*(-BB + SQRT(BB*BB + 4.d0*A7)) ! This is correct also !CC AHI =2.0*A7/(BB+SQRT(BB*BB + 4.d0*A7)) ! Avoid overflow when HI->0 DENM = BB+SQRT(BB*BB + 4.d0*A7) IF (DENM.LE.TINY) THEN ! Avoid overflow when HI->0 ABB = ABS(BB) DENM = (BB+ABB) + 2.0*A7/ABB ! Taylor expansion of SQRT ENDIF AHI = 2.0*A7/DENM ! ! *** SPECIATION & WATER CONTENT *************************************** ! MOLAL (1) = AHI ! HI MOLAL (3) = PSI1 + PSI4 + 2.D0*PSI2 ! NH4I MOLAL (5) = PSI2 ! SO4I MOLAL (6) = ZERO ! HSO4I MOLAL (7) = PSI3 + PSI1 ! NO3I CNH42S4 = CHI2 - PSI2 ! Solid (NH4)2SO4 CNH4NO3 = ZERO ! Solid NH4NO3 GHNO3 = CHI3 - PSI3 ! Gas HNO3 GNH3 = CHI4 - PSI4 ! Gas NH3 CALL CALCMR ! Water content ! ! *** CALCULATE ACTIVITIES OR TERMINATE INTERNAL LOOP ***************** ! IF (FRST.AND.CALAOU .OR. .NOT.FRST.AND.CALAIN) THEN CALL CALCACT ELSE GOTO 20 ENDIF 10 CONTINUE ! ! *** CALCULATE OBJECTIVE FUNCTION ************************************ ! 20 CONTINUE !CC FUNCD3= NH4I/HI/MAX(GNH3,TINY)/A4 - ONE FUNCD3= MOLAL(3)/MOLAL(1)/MAX(GNH3,TINY)/A4 - ONE RETURN ! ! *** END OF FUNCTION FUNCD3 ******************************************** ! END FUNCTION FUNCD3 !======================================================================= ! ! *** ISORROPIA CODE ! *** SUBROUTINE CALCD2 ! *** CASE D2 ! ! THE MAIN CHARACTERISTICS OF THIS REGIME ARE: ! 1. SULFATE POOR (SULRAT > 2.0) ! 2. THERE IS BOTH A LIQUID & SOLID PHASE ! 3. SOLIDS POSSIBLE : (NH4)2SO4 ! ! *** COPYRIGHT 1996-2006, UNIVERSITY OF MIAMI, CARNEGIE MELLON UNIVERSITY, ! *** GEORGIA INSTITUTE OF TECHNOLOGY ! *** WRITTEN BY ATHANASIOS NENES ! *** UPDATED BY CHRISTOS FOUNTOUKIS ! !======================================================================= ! SUBROUTINE CALCD2 USE ISRPIA USE SOLUT IMPLICIT DOUBLE PRECISION (A-H, O-Z) ! INCLUDE 'ISRPIA.INC' ! COMMON /SOLUT/ CHI1, CHI2, CHI3, CHI4, CHI5, CHI6, CHI7, CHI8, & ! PSI1, PSI2, PSI3, PSI4, PSI5, PSI6, PSI7, PSI8, & ! A1, A2, A3, A4, A5, A6, A7, A8 ! ! *** FIND DRY COMPOSITION ********************************************** ! CALL CALCD1A ! ! *** SETUP PARAMETERS ************************************************ ! CHI1 = CNH4NO3 ! Save from CALCD1 run CHI2 = CNH42S4 CHI3 = GHNO3 CHI4 = GNH3 ! PSI1 = CNH4NO3 ! ASSIGN INITIAL PSI's PSI2 = CNH42S4 PSI3 = ZERO PSI4 = ZERO ! MOLAL(5) = ZERO MOLAL(6) = ZERO MOLAL(3) = PSI1 MOLAL(7) = PSI1 CALL CALCMR ! Initial water ! CALAOU = .TRUE. ! Outer loop activity calculation flag PSI4LO = TINY ! Low limit PSI4HI = CHI4 ! High limit ! ! *** INITIAL VALUES FOR BISECTION ************************************ ! 60 X1 = PSI4LO Y1 = FUNCD2 (X1) IF (ABS(Y1).LE.EPS) RETURN YLO= Y1 ! Save Y-value at HI position ! ! *** ROOT TRACKING ; FOR THE RANGE OF HI AND LO ********************** ! DX = (PSI4HI-PSI4LO)/REAL(NDIV) DO 10 I=1,NDIV X2 = X1+DX Y2 = FUNCD2 (X2) IF (SIGN(1.d0,Y1)*SIGN(1.d0,Y2).LT.ZERO) THEN ! ! This is done, in case if Y(PSI4LO)>0, but Y(PSI4LO+DX) < 0 (i.e.undersat) ! IF (Y1 .LE. Y2) GOTO 20 ! (Y1*Y2.LT.ZERO) ENDIF X1 = X2 Y1 = Y2 10 CONTINUE ! ! *** NO SUBDIVISION WITH SOLUTION FOUND ! YHI= Y1 ! Save Y-value at Hi position IF (ABS(Y2) .LT. EPS) THEN ! X2 IS A SOLUTION RETURN ! ! *** { YLO, YHI } < 0.0 THE SOLUTION IS ALWAYS UNDERSATURATED WITH NH3 ! Physically I dont know when this might happen, but I have put this ! branch in for completeness. I assume there is no solution; all NO3 goes to the ! gas phase. ! ELSE IF (YLO.LT.ZERO .AND. YHI.LT.ZERO) THEN P4 = TINY ! PSI4LO ! CHI4 YY = FUNCD2(P4) GOTO 50 ! ! *** { YLO, YHI } > 0.0 THE SOLUTION IS ALWAYS SUPERSATURATED WITH NH3 ! This happens when Sul.Rat. = 2.0, so some NH4+ from sulfate evaporates ! and goes to the gas phase ; so I redefine the LO and HI limits of PSI4 ! and proceed again with root tracking. ! ELSE IF (YLO.GT.ZERO .AND. YHI.GT.ZERO) THEN PSI4HI = PSI4LO PSI4LO = PSI4LO - 0.1*(PSI1+PSI2) ! No solution; some NH3 evaporates IF (PSI4LO.LT.-(PSI1+PSI2)) THEN CALL PUSHERR (0001, 'CALCD2') ! WARNING ERROR: NO SOLUTION RETURN ELSE MOLAL(5) = ZERO MOLAL(6) = ZERO MOLAL(3) = PSI1 MOLAL(7) = PSI1 CALL CALCMR ! Initial water GOTO 60 ! Redo root tracking ENDIF ENDIF ! ! *** PERFORM BISECTION *********************************************** ! 20 DO 30 I=1,MAXIT X3 = 0.5*(X1+X2) Y3 = FUNCD2 (X3) IF (SIGN(1.d0,Y1)*SIGN(1.d0,Y3) .LE. ZERO) THEN ! (Y1*Y3 .LE. ZERO) Y2 = Y3 X2 = X3 ELSE Y1 = Y3 X1 = X3 ENDIF IF (ABS(X2-X1) .LE. EPS*ABS(X1)) GOTO 40 30 CONTINUE CALL PUSHERR (0002, 'CALCD2') ! WARNING ERROR: NO CONVERGENCE ! ! *** CONVERGED ; RETURN ********************************************** ! 40 X3 = MIN(X1,X2) ! 0.5*(X1+X2) ! Get "low" side, it's acidic soln. Y3 = FUNCD2 (X3) ! ! *** CALCULATE HSO4 SPECIATION AND RETURN ******************************* ! 50 CONTINUE IF (MOLAL(1).GT.TINY) THEN CALL CALCHS4 (MOLAL(1), MOLAL(5), ZERO, DELTA) MOLAL(1) = MOLAL(1) - DELTA ! H+ EFFECT MOLAL(5) = MOLAL(5) - DELTA ! SO4 EFFECT MOLAL(6) = DELTA ! HSO4 EFFECT ENDIF RETURN ! ! *** END OF SUBROUTINE CALCD2 ****************************************** ! END SUBROUTINE CALCD2 !======================================================================= ! ! *** ISORROPIA CODE ! *** FUNCTION FUNCD2 ! *** CASE D2 ! FUNCTION THAT SOLVES THE SYSTEM OF EQUATIONS FOR CASE D2 ; ! AND RETURNS THE VALUE OF THE ZEROED FUNCTION IN FUNCD2. ! !======================================================================= ! DOUBLE PRECISION FUNCTION FUNCD2 (P4) USE ISRPIA USE SOLUT IMPLICIT DOUBLE PRECISION (A-H, O-Z) ! INCLUDE 'ISRPIA.INC' ! COMMON /SOLUT/ CHI1, CHI2, CHI3, CHI4, CHI5, CHI6, CHI7, CHI8, & ! PSI1, PSI2, PSI3, PSI4, PSI5, PSI6, PSI7, PSI8, & ! A1, A2, A3, A4, A5, A6, A7, A8 ! ! *** SETUP PARAMETERS ************************************************ ! CALL RSTGAM ! Reset activity coefficients to 0.1 FRST = .TRUE. CALAIN = .TRUE. PSI4 = P4 PSI2 = CHI2 ! ! *** SOLVE EQUATIONS ; WITH ITERATIONS FOR ACTIVITY COEF. ************ ! DO 10 I=1,NSWEEP A2 = XK7*(WATER/GAMA(4))**3.0 A3 = XK4*R*TEMP*(WATER/GAMA(10))**2.0 A4 = (XK2/XKW)*R*TEMP*(GAMA(10)/GAMA(5))**2.0 A7 = XKW *RH*WATER*WATER ! IF (CHI2.GT.TINY .AND. WATER.GT.TINY) THEN PSI14 = PSI1+PSI4 CALL POLY3 (PSI14,0.25*PSI14**2.,-A2/4.D0, PSI2, ISLV) ! PSI2 IF (ISLV.EQ.0) THEN PSI2 = MIN (PSI2, CHI2) ELSE PSI2 = ZERO ENDIF ENDIF ! PSI3 = A3*A4*CHI3*(CHI4-PSI4) - PSI1*(2.D0*PSI2+PSI1+PSI4) PSI3 = PSI3/(A3*A4*(CHI4-PSI4) + 2.D0*PSI2+PSI1+PSI4) !cc PSI3 = MIN(MAX(PSI3, ZERO), CHI3) ! BB = PSI4-PSI3 ! (BB > 0, acidic solution, <0 alkaline) ! ! Do not change computation scheme for H+, all others did not work well. ! DENM = BB+SQRT(BB*BB + 4.d0*A7) IF (DENM.LE.TINY) THEN ! Avoid overflow when HI->0 ABB = ABS(BB) DENM = (BB+ABB) + 2.d0*A7/ABB ! Taylor expansion of SQRT ENDIF AHI = 2.d0*A7/DENM ! ! *** SPECIATION & WATER CONTENT *************************************** ! MOLAL (1) = AHI ! HI MOLAL (3) = PSI1 + PSI4 + 2.D0*PSI2 ! NH4 MOLAL (5) = PSI2 ! SO4 MOLAL (6) = ZERO ! HSO4 MOLAL (7) = PSI3 + PSI1 ! NO3 CNH42S4 = CHI2 - PSI2 ! Solid (NH4)2SO4 CNH4NO3 = ZERO ! Solid NH4NO3 GHNO3 = CHI3 - PSI3 ! Gas HNO3 GNH3 = CHI4 - PSI4 ! Gas NH3 CALL CALCMR ! Water content ! ! *** CALCULATE ACTIVITIES OR TERMINATE INTERNAL LOOP ***************** ! IF (FRST.AND.CALAOU .OR. .NOT.FRST.AND.CALAIN) THEN CALL CALCACT ELSE GOTO 20 ENDIF 10 CONTINUE ! ! *** CALCULATE OBJECTIVE FUNCTION ************************************ ! 20 CONTINUE !CC FUNCD2= NH4I/HI/MAX(GNH3,TINY)/A4 - ONE FUNCD2= MOLAL(3)/MOLAL(1)/MAX(GNH3,TINY)/A4 - ONE RETURN ! ! *** END OF FUNCTION FUNCD2 ******************************************** ! END FUNCTION FUNCD2 !======================================================================= ! ! *** ISORROPIA CODE ! *** SUBROUTINE CALCD1 ! *** CASE D1 ! ! THE MAIN CHARACTERISTICS OF THIS REGIME ARE: ! 1. SULFATE POOR (SULRAT > 2.0) ! 2. SOLID AEROSOL ONLY ! 3. SOLIDS POSSIBLE : (NH4)2SO4, NH4NO3 ! ! THERE ARE TWO REGIMES DEFINED BY RELATIVE HUMIDITY: ! 1. RH < MDRH ; ONLY SOLID PHASE POSSIBLE (SUBROUTINE CALCD1A) ! 2. RH >= MDRH ; LIQUID PHASE POSSIBLE (MDRH REGION) ! ! *** COPYRIGHT 1996-2006, UNIVERSITY OF MIAMI, CARNEGIE MELLON UNIVERSITY, ! *** GEORGIA INSTITUTE OF TECHNOLOGY ! *** WRITTEN BY ATHANASIOS NENES ! *** UPDATED BY CHRISTOS FOUNTOUKIS ! !======================================================================= ! SUBROUTINE CALCD1 USE ISRPIA IMPLICIT DOUBLE PRECISION (A-H, O-Z) ! INCLUDE 'ISRPIA.INC' EXTERNAL CALCD1A, CALCD2 ! ! *** REGIME DEPENDS UPON THE AMBIENT RELATIVE HUMIDITY ***************** ! IF (RH.LT.DRMASAN) THEN SCASE = 'D1 ; SUBCASE 1' ! SOLID PHASE ONLY POSSIBLE CALL CALCD1A SCASE = 'D1 ; SUBCASE 1' ELSE SCASE = 'D1 ; SUBCASE 2' ! LIQUID & SOLID PHASE POSSIBLE CALL CALCMDRH (RH, DRMASAN, DRNH4NO3, CALCD1A, CALCD2) SCASE = 'D1 ; SUBCASE 2' ENDIF ! RETURN ! ! *** END OF SUBROUTINE CALCD1 ****************************************** ! END SUBROUTINE CALCD1 !======================================================================= ! ! *** ISORROPIA CODE ! *** SUBROUTINE CALCD1A ! *** CASE D1 ; SUBCASE 1 ! ! THE MAIN CHARACTERISTICS OF THIS REGIME ARE: ! 1. SULFATE POOR (SULRAT > 2.0) ! 2. SOLID AEROSOL ONLY ! 3. SOLIDS POSSIBLE : (NH4)2SO4, NH4NO3 ! ! THE SOLID (NH4)2SO4 IS CALCULATED FROM THE SULFATES, WHILE NH4NO3 ! IS CALCULATED FROM NH3-HNO3 EQUILIBRIUM. 'ZE' IS THE AMOUNT OF ! NH4NO3 THAT VOLATIZES WHEN ALL POSSILBE NH4NO3 IS INITIALLY IN ! THE SOLID PHASE. ! ! *** COPYRIGHT 1996-2006, UNIVERSITY OF MIAMI, CARNEGIE MELLON UNIVERSITY, ! *** GEORGIA INSTITUTE OF TECHNOLOGY ! *** WRITTEN BY ATHANASIOS NENES ! *** UPDATED BY CHRISTOS FOUNTOUKIS ! !======================================================================= ! SUBROUTINE CALCD1A USE ISRPIA IMPLICIT DOUBLE PRECISION (A-H, O-Z) ! INCLUDE 'ISRPIA.INC' ! ! *** SETUP PARAMETERS ************************************************ ! PARM = XK10/(R*TEMP)/(R*TEMP) ! ! *** CALCULATE NH4NO3 THAT VOLATIZES ********************************* ! CNH42S4 = W(2) X = MAX(ZERO, MIN(W(3)-2.0*CNH42S4, W(4))) ! MAX NH4NO3 PS = MAX(W(3) - X - 2.0*CNH42S4, ZERO) OM = MAX(W(4) - X, ZERO) ! OMPS = OM+PS DIAK = SQRT(OMPS*OMPS + 4.0*PARM) ! DIAKRINOUSA ZE = MIN(X, 0.5*(-OMPS + DIAK)) ! THETIKI RIZA ! ! *** SPECIATION ******************************************************* ! CNH4NO3 = X - ZE ! Solid NH4NO3 GNH3 = PS + ZE ! Gas NH3 GHNO3 = OM + ZE ! Gas HNO3 ! RETURN ! ! *** END OF SUBROUTINE CALCD1A ***************************************** ! END SUBROUTINE CALCD1A !======================================================================= ! ! *** ISORROPIA CODE ! *** SUBROUTINE CALCG5 ! *** CASE G5 ! ! THE MAIN CHARACTERISTICS OF THIS REGIME ARE: ! 1. SULFATE POOR (SULRAT > 2.0) ; SODIUM POOR (SODRAT < 2.0) ! 2. THERE IS BOTH A LIQUID & SOLID PHASE ! 3. SOLIDS POSSIBLE : (NH4)2SO4, NH4CL, NA2SO4 ! ! *** COPYRIGHT 1996-2006, UNIVERSITY OF MIAMI, CARNEGIE MELLON UNIVERSITY, ! *** GEORGIA INSTITUTE OF TECHNOLOGY ! *** WRITTEN BY ATHANASIOS NENES ! *** UPDATED BY CHRISTOS FOUNTOUKIS ! !======================================================================= ! SUBROUTINE CALCG5 USE ISRPIA USE CASEG IMPLICIT DOUBLE PRECISION (A-H, O-Z) ! INCLUDE 'ISRPIA.INC' ! ! DOUBLE PRECISION LAMDA ! COMMON /CASEG/ CHI1, CHI2, CHI3, CHI4, CHI5, CHI6, LAMDA, & ! PSI1, PSI2, PSI3, PSI4, PSI5, PSI6, PSI7, & ! A1, A2, A3, A4, A5, A6, A7 ! ! *** SETUP PARAMETERS ************************************************ ! CALAOU = .TRUE. CHI1 = 0.5*W(1) CHI2 = MAX (W(2)-CHI1, ZERO) CHI3 = ZERO CHI4 = MAX (W(3)-2.D0*CHI2, ZERO) CHI5 = W(4) CHI6 = W(5) ! PSI1 = CHI1 PSI2 = CHI2 PSI6LO = TINY PSI6HI = CHI6-TINY ! MIN(CHI6-TINY, CHI4) ! WATER = CHI2/M0(4) + CHI1/M0(2) ! ! *** INITIAL VALUES FOR BISECTION ************************************ ! X1 = PSI6LO Y1 = FUNCG5A (X1) IF (CHI6.LE.TINY) GOTO 50 !cc IF (ABS(Y1).LE.EPS .OR. CHI6.LE.TINY) GOTO 50 !cc IF (WATER .LE. TINY) RETURN ! No water ! ! *** ROOT TRACKING ; FOR THE RANGE OF HI AND LO ********************** ! DX = (PSI6HI-PSI6LO)/REAL(NDIV) DO 10 I=1,NDIV X2 = X1+DX Y2 = FUNCG5A (X2) IF (SIGN(1.d0,Y1)*SIGN(1.d0,Y2).LT.ZERO) GOTO 20 ! (Y1*Y2.LT.ZERO) X1 = X2 Y1 = Y2 10 CONTINUE ! ! *** NO SUBDIVISION WITH SOLUTION; IF ABS(Y2) 2.0) ; SODIUM POOR (SODRAT < 2.0) ! 2. THERE IS BOTH A LIQUID & SOLID PHASE ! 3. SOLIDS POSSIBLE : (NH4)2SO4, NH4CL, NA2SO4 ! ! *** COPYRIGHT 1996-2006, UNIVERSITY OF MIAMI, CARNEGIE MELLON UNIVERSITY, ! *** GEORGIA INSTITUTE OF TECHNOLOGY ! *** WRITTEN BY ATHANASIOS NENES ! *** UPDATED BY CHRISTOS FOUNTOUKIS ! !======================================================================= ! DOUBLE PRECISION FUNCTION FUNCG5A (X) USE ISRPIA USE CASEG IMPLICIT DOUBLE PRECISION (A-H, O-Z) ! INCLUDE 'ISRPIA.INC' ! ! DOUBLE PRECISION LAMDA ! COMMON /CASEG/ CHI1, CHI2, CHI3, CHI4, CHI5, CHI6, LAMDA, & ! PSI1, PSI2, PSI3, PSI4, PSI5, PSI6, PSI7, & ! A1, A2, A3, A4, A5, A6, A7 ! ! *** SETUP PARAMETERS ************************************************ ! PSI6 = X FRST = .TRUE. CALAIN = .TRUE. ! ! *** SOLVE EQUATIONS ; WITH ITERATIONS FOR ACTIVITY COEF. ************ ! DO 10 I=1,NSWEEP ! A1 = XK5 *(WATER/GAMA(2))**3.0 A2 = XK7 *(WATER/GAMA(4))**3.0 A4 = (XK2/XKW)*R*TEMP*(GAMA(10)/GAMA(5))**2.0 A5 = XK4 *R*TEMP*(WATER/GAMA(10))**2.0 A6 = XK3 *R*TEMP*(WATER/GAMA(11))**2.0 AKK = A4*A6 ! ! CALCULATE DISSOCIATION QUANTITIES ! IF (CHI5.GE.TINY) THEN PSI5 = PSI6*CHI5/(A6/A5*(CHI6-PSI6) + PSI6) ELSE PSI5 = TINY ENDIF ! !CC IF(CHI4.GT.TINY) THEN IF(W(2).GT.TINY) THEN ! Accounts for NH3 evaporation BB =-(CHI4 + PSI6 + PSI5 + 1.d0/A4) CC = CHI4*(PSI5+PSI6) - 2.d0*PSI2/A4 DD = MAX(BB*BB-4.d0*CC,ZERO) ! Patch proposed by Uma Shankar, 19/11/01 PSI4 =0.5d0*(-BB - SQRT(DD)) ELSE PSI4 = TINY ENDIF ! ! *** CALCULATE SPECIATION ******************************************** ! MOLAL (2) = 2.0D0*PSI1 ! NAI MOLAL (3) = 2.0*PSI2 + PSI4 ! NH4I MOLAL (4) = PSI6 ! CLI MOLAL (5) = PSI2 + PSI1 ! SO4I MOLAL (6) = ZERO MOLAL (7) = PSI5 ! NO3I ! SMIN = 2.d0*MOLAL(5)+MOLAL(7)+MOLAL(4)-MOLAL(2)-MOLAL(3) CALL CALCPH (SMIN, HI, OHI) MOLAL (1) = HI ! GNH3 = MAX(CHI4 - PSI4, TINY) ! Gas NH3 GHNO3 = MAX(CHI5 - PSI5, TINY) ! Gas HNO3 GHCL = MAX(CHI6 - PSI6, TINY) ! Gas HCl ! CNH42S4 = ZERO ! Solid (NH4)2SO4 CNH4NO3 = ZERO ! Solid NH4NO3 CNH4CL = ZERO ! Solid NH4Cl ! CALL CALCMR ! Water content ! ! *** CALCULATE ACTIVITIES OR TERMINATE INTERNAL LOOP ***************** ! IF (FRST.AND.CALAOU .OR. .NOT.FRST.AND.CALAIN) THEN CALL CALCACT ELSE GOTO 20 ENDIF 10 CONTINUE ! ! *** CALCULATE FUNCTION VALUE FOR OUTER LOOP *************************** ! 20 FUNCG5A = MOLAL(1)*MOLAL(4)/GHCL/A6 - ONE !CC FUNCG5A = MOLAL(3)*MOLAL(4)/GHCL/GNH3/A6/A4 - ONE ! RETURN ! ! *** END OF FUNCTION FUNCG5A ******************************************* ! END FUNCTION FUNCG5A !======================================================================= ! ! *** ISORROPIA CODE ! *** SUBROUTINE CALCG4 ! *** CASE G4 ! ! THE MAIN CHARACTERISTICS OF THIS REGIME ARE: ! 1. SULFATE POOR (SULRAT > 2.0) ; SODIUM POOR (SODRAT < 2.0) ! 2. THERE IS BOTH A LIQUID & SOLID PHASE ! 3. SOLIDS POSSIBLE : (NH4)2SO4, NH4CL, NA2SO4 ! ! *** COPYRIGHT 1996-2006, UNIVERSITY OF MIAMI, CARNEGIE MELLON UNIVERSITY, ! *** GEORGIA INSTITUTE OF TECHNOLOGY ! *** WRITTEN BY ATHANASIOS NENES ! *** UPDATED BY CHRISTOS FOUNTOUKIS ! !======================================================================= ! SUBROUTINE CALCG4 USE ISRPIA USE CASEG IMPLICIT DOUBLE PRECISION (A-H, O-Z) ! INCLUDE 'ISRPIA.INC' ! ! DOUBLE PRECISION LAMDA ! COMMON /CASEG/ CHI1, CHI2, CHI3, CHI4, CHI5, CHI6, LAMDA, & ! PSI1, PSI2, PSI3, PSI4, PSI5, PSI6, PSI7, & ! A1, A2, A3, A4, A5, A6, A7 ! ! *** SETUP PARAMETERS ************************************************ ! CALAOU = .TRUE. CHI1 = 0.5*W(1) CHI2 = MAX (W(2)-CHI1, ZERO) CHI3 = ZERO CHI4 = MAX (W(3)-2.D0*CHI2, ZERO) CHI5 = W(4) CHI6 = W(5) ! PSI2 = CHI2 PSI6LO = TINY PSI6HI = CHI6-TINY ! MIN(CHI6-TINY, CHI4) ! WATER = CHI2/M0(4) + CHI1/M0(2) ! ! *** INITIAL VALUES FOR BISECTION ************************************ ! X1 = PSI6LO Y1 = FUNCG4A (X1) IF (CHI6.LE.TINY) GOTO 50 !CC IF (ABS(Y1).LE.EPS .OR. CHI6.LE.TINY .OR. WATER .LE. TINY) GOTO 50 !CC IF (WATER .LE. TINY) RETURN ! No water ! ! *** ROOT TRACKING ; FOR THE RANGE OF HI AND LO ********************** ! DX = (PSI6HI-PSI6LO)/REAL(NDIV) DO 10 I=1,NDIV X2 = X1+DX Y2 = FUNCG4A (X2) IF (SIGN(1.d0,Y1)*SIGN(1.d0,Y2).LT.ZERO) GOTO 20 ! (Y1*Y2.LT.ZERO) X1 = X2 Y1 = Y2 10 CONTINUE ! ! *** NO SUBDIVISION WITH SOLUTION; IF ABS(Y2) 2.0) ; SODIUM POOR (SODRAT < 2.0) ! 2. THERE IS BOTH A LIQUID & SOLID PHASE ! 3. SOLIDS POSSIBLE : (NH4)2SO4, NH4CL, NA2SO4 ! ! *** COPYRIGHT 1996-2006, UNIVERSITY OF MIAMI, CARNEGIE MELLON UNIVERSITY, ! *** GEORGIA INSTITUTE OF TECHNOLOGY ! *** WRITTEN BY ATHANASIOS NENES ! *** UPDATED BY CHRISTOS FOUNTOUKIS ! !======================================================================= ! DOUBLE PRECISION FUNCTION FUNCG4A (X) USE ISRPIA USE CASEG IMPLICIT DOUBLE PRECISION (A-H, O-Z) ! INCLUDE 'ISRPIA.INC' ! ! DOUBLE PRECISION LAMDA, NAI, NH4I, NO3I DOUBLE PRECISION NAI, NH4I, NO3I ! COMMON /CASEG/ CHI1, CHI2, CHI3, CHI4, CHI5, CHI6, LAMDA, & ! PSI1, PSI2, PSI3, PSI4, PSI5, PSI6, PSI7, & ! A1, A2, A3, A4, A5, A6, A7 ! ! *** SETUP PARAMETERS ************************************************ ! PSI6 = X PSI1 = CHI1 FRST = .TRUE. CALAIN = .TRUE. ! ! *** SOLVE EQUATIONS ; WITH ITERATIONS FOR ACTIVITY COEF. ************ ! DO 10 I=1,NSWEEP ! A1 = XK5 *(WATER/GAMA(2))**3.0 A2 = XK7 *(WATER/GAMA(4))**3.0 A4 = (XK2/XKW)*R*TEMP*(GAMA(10)/GAMA(5))**2.0 A5 = XK4 *R*TEMP*(WATER/GAMA(10))**2.0 A6 = XK3 *R*TEMP*(WATER/GAMA(11))**2.0 ! ! CALCULATE DISSOCIATION QUANTITIES ! IF (CHI5.GE.TINY) THEN PSI5 = PSI6*CHI5/(A6/A5*(CHI6-PSI6) + PSI6) ELSE PSI5 = TINY ENDIF ! !CC IF(CHI4.GT.TINY) THEN IF(W(2).GT.TINY) THEN ! Accounts for NH3 evaporation BB =-(CHI4 + PSI6 + PSI5 + 1.d0/A4) CC = CHI4*(PSI5+PSI6) - 2.d0*PSI2/A4 DD = MAX(BB*BB-4.d0*CC,ZERO) ! Patch proposed by Uma shankar, 19/11/2001 PSI4 =0.5d0*(-BB - SQRT(DD)) ELSE PSI4 = TINY ENDIF ! ! CALCULATE CONCENTRATIONS ! NH4I = 2.0*PSI2 + PSI4 CLI = PSI6 SO4I = PSI2 + PSI1 NO3I = PSI5 NAI = 2.0D0*PSI1 ! CALL CALCPH(2.d0*SO4I+NO3I+CLI-NAI-NH4I, HI, OHI) ! ! *** Na2SO4 DISSOLUTION ! IF (CHI1.GT.TINY .AND. WATER.GT.TINY) THEN ! PSI1 CALL POLY3 (PSI2, ZERO, -A1/4.D0, PSI1, ISLV) IF (ISLV.EQ.0) THEN PSI1 = MIN (PSI1, CHI1) ELSE PSI1 = ZERO ENDIF ELSE PSI1 = ZERO ENDIF ! ! *** SAVE CONCENTRATIONS IN MOLAL ARRAY ****************************** ! MOLAL (1) = HI MOLAL (2) = NAI MOLAL (3) = NH4I MOLAL (4) = CLI MOLAL (5) = SO4I MOLAL (6) = ZERO MOLAL (7) = NO3I ! ! *** CALCULATE GAS / SOLID SPECIES (LIQUID IN MOLAL ALREADY) ********* ! GNH3 = MAX(CHI4 - PSI4, TINY) GHNO3 = MAX(CHI5 - PSI5, TINY) GHCL = MAX(CHI6 - PSI6, TINY) ! CNH42S4 = ZERO CNH4NO3 = ZERO CNH4CL = ZERO CNA2SO4 = MAX(CHI1-PSI1,ZERO) ! ! *** CALCULATE MOLALR ARRAY, WATER AND ACTIVITIES ********************** ! CALL CALCMR ! ! *** CALCULATE ACTIVITIES OR TERMINATE INTERNAL LOOP ***************** ! IF (FRST.AND.CALAOU .OR. .NOT.FRST.AND.CALAIN) THEN CALL CALCACT ELSE GOTO 20 ENDIF 10 CONTINUE ! ! *** CALCULATE FUNCTION VALUE FOR OUTER LOOP *************************** ! 20 FUNCG4A = MOLAL(1)*MOLAL(4)/GHCL/A6 - ONE !CC FUNCG4A = MOLAL(3)*MOLAL(4)/GHCL/GNH3/A6/A4 - ONE ! RETURN ! ! *** END OF FUNCTION FUNCG4A ******************************************* ! END FUNCTION FUNCG4A !======================================================================= ! ! *** ISORROPIA CODE ! *** SUBROUTINE CALCG3 ! *** CASE G3 ! ! THE MAIN CHARACTERISTICS OF THIS REGIME ARE: ! 1. SULFATE POOR (SULRAT > 2.0) ; SODIUM POOR (SODRAT < 2.0) ! 2. LIQUID & SOLID PHASE ARE BOTH POSSIBLE ! 3. SOLIDS POSSIBLE : (NH4)2SO4, NH4CL, NA2SO4 ! ! *** COPYRIGHT 1996-2006, UNIVERSITY OF MIAMI, CARNEGIE MELLON UNIVERSITY, ! *** GEORGIA INSTITUTE OF TECHNOLOGY ! *** WRITTEN BY ATHANASIOS NENES ! *** UPDATED BY CHRISTOS FOUNTOUKIS ! !======================================================================= ! SUBROUTINE CALCG3 USE ISRPIA IMPLICIT DOUBLE PRECISION (A-H, O-Z) ! INCLUDE 'ISRPIA.INC' EXTERNAL CALCG1A, CALCG4 ! ! *** REGIME DEPENDS ON THE EXISTANCE OF WATER AND OF THE RH ************ ! IF (W(4).GT.TINY .AND. W(5).GT.TINY) THEN ! NO3,CL EXIST, WATER POSSIBLE SCASE = 'G3 ; SUBCASE 1' CALL CALCG3A SCASE = 'G3 ; SUBCASE 1' ELSE ! NO3, CL NON EXISTANT SCASE = 'G1 ; SUBCASE 1' CALL CALCG1A SCASE = 'G1 ; SUBCASE 1' ENDIF ! IF (WATER.LE.TINY) THEN IF (RH.LT.DRMG3) THEN ! ONLY SOLIDS WATER = TINY DO 10 I=1,NIONS MOLAL(I) = ZERO 10 CONTINUE CALL CALCG1A SCASE = 'G3 ; SUBCASE 2' RETURN ELSE SCASE = 'G3 ; SUBCASE 3' ! MDRH REGION (NA2SO4, NH42S4) CALL CALCMDRH (RH, DRMG3, DRNH42S4, CALCG1A, CALCG4) SCASE = 'G3 ; SUBCASE 3' ENDIF ENDIF ! RETURN ! ! *** END OF SUBROUTINE CALCG3 ****************************************** ! END SUBROUTINE CALCG3 !======================================================================= ! ! *** ISORROPIA CODE ! *** SUBROUTINE CALCG3A ! *** CASE G3 ; SUBCASE 1 ! ! THE MAIN CHARACTERISTICS OF THIS REGIME ARE: ! 1. SULFATE POOR (SULRAT > 2.0) ; SODIUM POOR (SODRAT < 2.0) ! 2. THERE IS BOTH A LIQUID & SOLID PHASE ! 3. SOLIDS POSSIBLE : (NH4)2SO4, NH4CL, NA2SO4 ! ! *** COPYRIGHT 1996-2006, UNIVERSITY OF MIAMI, CARNEGIE MELLON UNIVERSITY, ! *** GEORGIA INSTITUTE OF TECHNOLOGY ! *** WRITTEN BY ATHANASIOS NENES ! *** UPDATED BY CHRISTOS FOUNTOUKIS ! !======================================================================= ! SUBROUTINE CALCG3A USE ISRPIA USE CASEG IMPLICIT DOUBLE PRECISION (A-H, O-Z) ! INCLUDE 'ISRPIA.INC' ! ! DOUBLE PRECISION LAMDA ! COMMON /CASEG/ CHI1, CHI2, CHI3, CHI4, CHI5, CHI6, LAMDA, & ! PSI1, PSI2, PSI3, PSI4, PSI5, PSI6, PSI7, & ! A1, A2, A3, A4, A5, A6, A7 ! ! *** SETUP PARAMETERS ************************************************ ! CALAOU = .TRUE. CHI1 = 0.5*W(1) CHI2 = MAX (W(2)-CHI1, ZERO) CHI3 = ZERO CHI4 = MAX (W(3)-2.D0*CHI2, ZERO) CHI5 = W(4) CHI6 = W(5) ! PSI6LO = TINY PSI6HI = CHI6-TINY ! MIN(CHI6-TINY, CHI4) ! WATER = TINY ! ! *** INITIAL VALUES FOR BISECTION ************************************ ! X1 = PSI6LO Y1 = FUNCG3A (X1) IF (CHI6.LE.TINY) GOTO 50 !CC IF (ABS(Y1).LE.EPS .OR. CHI6.LE.TINY .OR. WATER .LE. TINY) GOTO 50 !CC IF (WATER .LE. TINY) RETURN ! No water ! ! *** ROOT TRACKING ; FOR THE RANGE OF HI AND LO ********************** ! DX = (PSI6HI-PSI6LO)/REAL(NDIV) DO 10 I=1,NDIV X2 = X1+DX Y2 = FUNCG3A (X2) ! IF (SIGN(1.d0,Y1)*SIGN(1.d0,Y2).LT.ZERO) GOTO 20 ! (Y1*Y2.LT.ZERO) X1 = X2 Y1 = Y2 10 CONTINUE ! ! *** NO SUBDIVISION WITH SOLUTION; IF ABS(Y2) 2.0) ; SODIUM POOR (SODRAT < 2.0) ! 2. THERE IS BOTH A LIQUID & SOLID PHASE ! 3. SOLIDS POSSIBLE : (NH4)2SO4, NH4CL, NA2SO4 ! ! *** COPYRIGHT 1996-2006, UNIVERSITY OF MIAMI, CARNEGIE MELLON UNIVERSITY, ! *** GEORGIA INSTITUTE OF TECHNOLOGY ! *** WRITTEN BY ATHANASIOS NENES ! *** UPDATED BY CHRISTOS FOUNTOUKIS ! !======================================================================= ! DOUBLE PRECISION FUNCTION FUNCG3A (X) USE ISRPIA USE CASEG IMPLICIT DOUBLE PRECISION (A-H, O-Z) ! INCLUDE 'ISRPIA.INC' ! ! DOUBLE PRECISION LAMDA ! COMMON /CASEG/ CHI1, CHI2, CHI3, CHI4, CHI5, CHI6, LAMDA, & ! PSI1, PSI2, PSI3, PSI4, PSI5, PSI6, PSI7, & ! A1, A2, A3, A4, A5, A6, A7 ! ! *** SETUP PARAMETERS ************************************************ ! PSI6 = X PSI2 = CHI2 FRST = .TRUE. CALAIN = .TRUE. ! ! *** SOLVE EQUATIONS ; WITH ITERATIONS FOR ACTIVITY COEF. ************ ! DO 10 I=1,NSWEEP ! A1 = XK5 *(WATER/GAMA(2))**3.0 A2 = XK7 *(WATER/GAMA(4))**3.0 A4 = (XK2/XKW)*R*TEMP*(GAMA(10)/GAMA(5))**2.0 A5 = XK4 *R*TEMP*(WATER/GAMA(10))**2.0 A6 = XK3 *R*TEMP*(WATER/GAMA(11))**2.0 ! ! CALCULATE DISSOCIATION QUANTITIES ! IF (CHI5.GE.TINY) THEN PSI5 = PSI6*CHI5/(A6/A5*(CHI6-PSI6) + PSI6) ELSE PSI5 = TINY ENDIF ! !CC IF(CHI4.GT.TINY) THEN IF(W(2).GT.TINY) THEN ! Accounts for NH3 evaporation BB =-(CHI4 + PSI6 + PSI5 + 1.d0/A4) CC = CHI4*(PSI5+PSI6) - 2.d0*PSI2/A4 DD = MAX(BB*BB-4.d0*CC,ZERO) ! Patch proposed by Uma Shankar, 19/11/01 PSI4 =0.5d0*(-BB - SQRT(DD)) ELSE PSI4 = TINY ENDIF ! IF (CHI2.GT.TINY .AND. WATER.GT.TINY) THEN CALL POLY3 (PSI4, PSI4*PSI4/4.D0, -A2/4.D0, PSI20, ISLV) IF (ISLV.EQ.0) PSI2 = MIN (PSI20, CHI2) ENDIF ! ! *** CALCULATE GAS / SOLID SPECIES (LIQUID IN MOLAL ALREADY) ********* ! MOLAL (2) = ZERO ! Na MOLAL (3) = 2.0*PSI2 + PSI4 ! NH4I MOLAL (4) = PSI6 ! CLI MOLAL (5) = PSI2 ! SO4I MOLAL (6) = ZERO ! HSO4 MOLAL (7) = PSI5 ! NO3I ! SMIN = 2.d0*MOLAL(5)+MOLAL(7)+MOLAL(4)-MOLAL(2)-MOLAL(3) CALL CALCPH (SMIN, HI, OHI) MOLAL (1) = HI ! GNH3 = MAX(CHI4 - PSI4, TINY) ! Gas NH3 GHNO3 = MAX(CHI5 - PSI5, TINY) ! Gas HNO3 GHCL = MAX(CHI6 - PSI6, TINY) ! Gas HCl ! CNH42S4 = CHI2 - PSI2 ! Solid (NH4)2SO4 CNH4NO3 = ZERO ! Solid NH4NO3 CNH4CL = ZERO ! Solid NH4Cl ! CALL CALCMR ! Water content ! ! *** CALCULATE ACTIVITIES OR TERMINATE INTERNAL LOOP ***************** ! IF (FRST.AND.CALAOU .OR. .NOT.FRST.AND.CALAIN) THEN CALL CALCACT ELSE GOTO 20 ENDIF 10 CONTINUE ! ! *** CALCULATE FUNCTION VALUE FOR OUTER LOOP *************************** ! 20 FUNCG3A = MOLAL(1)*MOLAL(4)/GHCL/A6 - ONE !CC FUNCG3A = MOLAL(3)*MOLAL(4)/GHCL/GNH3/A6/A4 - ONE ! RETURN ! ! *** END OF FUNCTION FUNCG3A ******************************************* ! END FUNCTION FUNCG3A !======================================================================= ! ! *** ISORROPIA CODE ! *** SUBROUTINE CALCG2 ! *** CASE G2 ! ! THE MAIN CHARACTERISTICS OF THIS REGIME ARE: ! 1. SULFATE POOR (SULRAT > 2.0) ; SODIUM POOR (SODRAT < 2.0) ! 2. LIQUID & SOLID PHASE ARE BOTH POSSIBLE ! 3. SOLIDS POSSIBLE : (NH4)2SO4, NH4CL, NA2SO4 ! ! *** COPYRIGHT 1996-2006, UNIVERSITY OF MIAMI, CARNEGIE MELLON UNIVERSITY, ! *** GEORGIA INSTITUTE OF TECHNOLOGY ! *** WRITTEN BY ATHANASIOS NENES ! *** UPDATED BY CHRISTOS FOUNTOUKIS ! !======================================================================= ! SUBROUTINE CALCG2 USE ISRPIA IMPLICIT DOUBLE PRECISION (A-H, O-Z) ! INCLUDE 'ISRPIA.INC' EXTERNAL CALCG1A, CALCG3A, CALCG4 ! ! *** REGIME DEPENDS ON THE EXISTANCE OF NITRATES *********************** ! IF (W(4).GT.TINY) THEN ! NO3 EXISTS, WATER POSSIBLE SCASE = 'G2 ; SUBCASE 1' CALL CALCG2A SCASE = 'G2 ; SUBCASE 1' ELSE ! NO3 NON EXISTANT, WATER NOT POSSIBLE SCASE = 'G1 ; SUBCASE 1' CALL CALCG1A SCASE = 'G1 ; SUBCASE 1' ENDIF ! ! *** REGIME DEPENDS ON THE EXISTANCE OF WATER AND OF THE RH ************ ! IF (WATER.LE.TINY) THEN IF (RH.LT.DRMG2) THEN ! ONLY SOLIDS WATER = TINY DO 10 I=1,NIONS MOLAL(I) = ZERO 10 CONTINUE CALL CALCG1A SCASE = 'G2 ; SUBCASE 2' ELSE IF (W(5).GT. TINY) THEN SCASE = 'G2 ; SUBCASE 3' ! MDRH (NH4CL, NA2SO4, NH42S4) CALL CALCMDRH (RH, DRMG2, DRNH4CL, CALCG1A, CALCG3A) SCASE = 'G2 ; SUBCASE 3' ENDIF IF (WATER.LE.TINY .AND. RH.GE.DRMG3) THEN SCASE = 'G2 ; SUBCASE 4' ! MDRH (NA2SO4, NH42S4) CALL CALCMDRH (RH, DRMG3, DRNH42S4, CALCG1A, CALCG4) SCASE = 'G2 ; SUBCASE 4' ELSE WATER = TINY DO 20 I=1,NIONS MOLAL(I) = ZERO 20 CONTINUE CALL CALCG1A SCASE = 'G2 ; SUBCASE 2' ENDIF ENDIF ENDIF ! RETURN ! ! *** END OF SUBROUTINE CALCG2 ****************************************** ! END SUBROUTINE CALCG2 !======================================================================= ! ! *** ISORROPIA CODE ! *** SUBROUTINE CALCG2A ! *** CASE G2 ; SUBCASE 1 ! ! THE MAIN CHARACTERISTICS OF THIS REGIME ARE: ! 1. SULFATE POOR (SULRAT > 2.0) ; SODIUM POOR (SODRAT < 2.0) ! 2. THERE IS BOTH A LIQUID & SOLID PHASE ! 3. SOLIDS POSSIBLE : (NH4)2SO4, NH4CL, NA2SO4 ! ! *** COPYRIGHT 1996-2006, UNIVERSITY OF MIAMI, CARNEGIE MELLON UNIVERSITY, ! *** GEORGIA INSTITUTE OF TECHNOLOGY ! *** WRITTEN BY ATHANASIOS NENES ! *** UPDATED BY CHRISTOS FOUNTOUKIS ! !======================================================================= ! SUBROUTINE CALCG2A USE ISRPIA USE CASEG IMPLICIT DOUBLE PRECISION (A-H, O-Z) ! INCLUDE 'ISRPIA.INC' ! ! DOUBLE PRECISION LAMDA ! COMMON /CASEG/ CHI1, CHI2, CHI3, CHI4, CHI5, CHI6, LAMDA, & ! PSI1, PSI2, PSI3, PSI4, PSI5, PSI6, PSI7, & ! A1, A2, A3, A4, A5, A6, A7 ! ! *** SETUP PARAMETERS ************************************************ ! CALAOU = .TRUE. CHI1 = 0.5*W(1) CHI2 = MAX (W(2)-CHI1, ZERO) CHI3 = ZERO CHI4 = MAX (W(3)-2.D0*CHI2, ZERO) CHI5 = W(4) CHI6 = W(5) ! PSI6LO = TINY PSI6HI = CHI6-TINY ! WATER = TINY ! ! *** INITIAL VALUES FOR BISECTION ************************************ ! X1 = PSI6LO Y1 = FUNCG2A (X1) IF (CHI6.LE.TINY) GOTO 50 !CC IF (ABS(Y1).LE.EPS .OR. CHI6.LE.TINY) GOTO 50 !CC IF (WATER .LE. TINY) GOTO 50 ! No water ! ! *** ROOT TRACKING ; FOR THE RANGE OF HI AND LO ********************** ! DX = (PSI6HI-PSI6LO)/REAL(NDIV) DO 10 I=1,NDIV X2 = X1+DX Y2 = FUNCG2A (X2) IF (SIGN(1.d0,Y1)*SIGN(1.d0,Y2).LT.ZERO) GOTO 20 ! (Y1*Y2.LT.ZERO) X1 = X2 Y1 = Y2 10 CONTINUE ! ! *** NO SUBDIVISION WITH SOLUTION; IF ABS(Y2) 2.0) ; SODIUM POOR (SODRAT < 2.0) ! 2. THERE IS BOTH A LIQUID & SOLID PHASE ! 3. SOLIDS POSSIBLE : (NH4)2SO4, NH4CL, NA2SO4 ! ! *** COPYRIGHT 1996-2006, UNIVERSITY OF MIAMI, CARNEGIE MELLON UNIVERSITY, ! *** GEORGIA INSTITUTE OF TECHNOLOGY ! *** WRITTEN BY ATHANASIOS NENES ! *** UPDATED BY CHRISTOS FOUNTOUKIS ! !======================================================================= ! DOUBLE PRECISION FUNCTION FUNCG2A (X) USE ISRPIA USE CASEG IMPLICIT DOUBLE PRECISION (A-H, O-Z) ! INCLUDE 'ISRPIA.INC' ! ! DOUBLE PRECISION LAMDA ! COMMON /CASEG/ CHI1, CHI2, CHI3, CHI4, CHI5, CHI6, LAMDA, & ! PSI1, PSI2, PSI3, PSI4, PSI5, PSI6, PSI7, & ! A1, A2, A3, A4, A5, A6, A7 ! ! *** SETUP PARAMETERS ************************************************ ! PSI6 = X PSI2 = CHI2 PSI3 = ZERO FRST = .TRUE. CALAIN = .TRUE. ! ! *** SOLVE EQUATIONS ; WITH ITERATIONS FOR ACTIVITY COEF. ************ ! DO 10 I=1,NSWEEP ! A1 = XK5 *(WATER/GAMA(2))**3.0 A2 = XK7 *(WATER/GAMA(4))**3.0 A4 = (XK2/XKW)*R*TEMP*(GAMA(10)/GAMA(5))**2.0 A5 = XK4 *R*TEMP*(WATER/GAMA(10))**2.0 A6 = XK3 *R*TEMP*(WATER/GAMA(11))**2.0 ! DENO = MAX(CHI6-PSI6-PSI3, ZERO) PSI5 = CHI5/((A6/A5)*(DENO/PSI6) + ONE) ! PSI4 = MIN(PSI5+PSI6,CHI4) ! IF (CHI2.GT.TINY .AND. WATER.GT.TINY) THEN CALL POLY3 (PSI4, PSI4*PSI4/4.D0, -A2/4.D0, PSI20, ISLV) IF (ISLV.EQ.0) PSI2 = MIN (PSI20, CHI2) ENDIF ! ! *** SAVE CONCENTRATIONS IN MOLAL ARRAY ****************************** ! MOLAL (2) = ZERO ! NA MOLAL (3) = 2.0*PSI2 + PSI4 ! NH4I MOLAL (4) = PSI6 ! CLI MOLAL (5) = PSI2 ! SO4I MOLAL (6) = ZERO ! HSO4 MOLAL (7) = PSI5 ! NO3I ! !CC MOLAL (1) = MAX(CHI5 - PSI5, TINY)*A5/PSI5 ! HI SMIN = 2.d0*MOLAL(5)+MOLAL(7)+MOLAL(4)-MOLAL(2)-MOLAL(3) CALL CALCPH (SMIN, HI, OHI) MOLAL (1) = HI ! ! *** CALCULATE GAS / SOLID SPECIES (LIQUID IN MOLAL ALREADY) ********* ! GNH3 = MAX(CHI4 - PSI4, TINY) GHNO3 = MAX(CHI5 - PSI5, TINY) GHCL = MAX(CHI6 - PSI6, TINY) ! CNH42S4 = MAX(CHI2 - PSI2, ZERO) CNH4NO3 = ZERO ! ! *** NH4Cl(s) calculations ! A3 = XK6 /(R*TEMP*R*TEMP) IF (GNH3*GHCL.GT.A3) THEN DELT = MIN(GNH3, GHCL) BB = -(GNH3+GHCL) CC = GNH3*GHCL-A3 DD = BB*BB - 4.D0*CC PSI31 = 0.5D0*(-BB + SQRT(DD)) PSI32 = 0.5D0*(-BB - SQRT(DD)) IF (DELT-PSI31.GT.ZERO .AND. PSI31.GT.ZERO) THEN PSI3 = PSI31 ELSEIF (DELT-PSI32.GT.ZERO .AND. PSI32.GT.ZERO) THEN PSI3 = PSI32 ELSE PSI3 = ZERO ENDIF ELSE PSI3 = ZERO ENDIF ! ! *** CALCULATE GAS / SOLID SPECIES (LIQUID IN MOLAL ALREADY) ********* ! GNH3 = MAX(GNH3 - PSI3, TINY) GHCL = MAX(GHCL - PSI3, TINY) CNH4CL = PSI3 ! ! *** CALCULATE MOLALR ARRAY, WATER AND ACTIVITIES ********************** ! CALL CALCMR ! ! *** CALCULATE ACTIVITIES OR TERMINATE INTERNAL LOOP ***************** ! IF (FRST.AND.CALAOU .OR. .NOT.FRST.AND.CALAIN) THEN CALL CALCACT ELSE GOTO 20 ENDIF 10 CONTINUE ! ! *** CALCULATE FUNCTION VALUE FOR OUTER LOOP *************************** ! 20 IF (CHI4.LE.TINY) THEN FUNCG2A = MOLAL(1)*MOLAL(4)/GHCL/A6 - ONE ELSE FUNCG2A = MOLAL(3)*MOLAL(4)/GHCL/GNH3/A6/A4 - ONE ENDIF ! RETURN ! ! *** END OF FUNCTION FUNCG2A ******************************************* ! END FUNCTION FUNCG2A !======================================================================= ! ! *** ISORROPIA CODE ! *** SUBROUTINE CALCG1 ! *** CASE G1 ! ! THE MAIN CHARACTERISTICS OF THIS REGIME ARE: ! 1. SULFATE POOR (SULRAT > 2.0) ; SODIUM POOR (SODRAT < 2.0) ! 2. SOLID AEROSOL ONLY ! 3. SOLIDS POSSIBLE : (NH4)2SO4, NH4NO3, NH4CL, NA2SO4 ! ! THERE ARE TWO POSSIBLE REGIMES HERE, DEPENDING ON RELATIVE HUMIDITY: ! 1. WHEN RH >= MDRH ; LIQUID PHASE POSSIBLE (MDRH REGION) ! 2. WHEN RH < MDRH ; ONLY SOLID PHASE POSSIBLE (SUBROUTINE CALCG1A) ! ! *** COPYRIGHT 1996-2006, UNIVERSITY OF MIAMI, CARNEGIE MELLON UNIVERSITY, ! *** GEORGIA INSTITUTE OF TECHNOLOGY ! *** WRITTEN BY ATHANASIOS NENES ! *** UPDATED BY CHRISTOS FOUNTOUKIS ! !======================================================================= ! SUBROUTINE CALCG1 USE ISRPIA IMPLICIT DOUBLE PRECISION (A-H, O-Z) ! INCLUDE 'ISRPIA.INC' EXTERNAL CALCG1A, CALCG2A ! ! *** REGIME DEPENDS UPON THE AMBIENT RELATIVE HUMIDITY ***************** ! IF (RH.LT.DRMG1) THEN SCASE = 'G1 ; SUBCASE 1' CALL CALCG1A ! SOLID PHASE ONLY POSSIBLE SCASE = 'G1 ; SUBCASE 1' ELSE SCASE = 'G1 ; SUBCASE 2' ! LIQUID & SOLID PHASE POSSIBLE CALL CALCMDRH (RH, DRMG1, DRNH4NO3, CALCG1A, CALCG2A) SCASE = 'G1 ; SUBCASE 2' ENDIF ! RETURN ! ! *** END OF SUBROUTINE CALCG1 ****************************************** ! END SUBROUTINE CALCG1 !======================================================================= ! ! *** ISORROPIA CODE ! *** SUBROUTINE CALCG1A ! *** CASE G1 ; SUBCASE 1 ! ! THE MAIN CHARACTERISTICS OF THIS REGIME ARE: ! 1. SULFATE POOR (SULRAT > 2.0) ! 2. SOLID AEROSOL ONLY ! 3. SOLIDS POSSIBLE : (NH4)2SO4, NH4NO3 ! ! SOLID (NH4)2SO4 IS CALCULATED FROM THE SULFATES, WHILE NH4NO3 ! IS CALCULATED FROM NH3-HNO3 EQUILIBRIUM. 'ZE' IS THE AMOUNT OF ! NH4NO3 THAT VOLATIZES WHEN ALL POSSILBE NH4NO3 IS INITIALLY IN ! THE SOLID PHASE. ! ! *** COPYRIGHT 1996-2006, UNIVERSITY OF MIAMI, CARNEGIE MELLON UNIVERSITY, ! *** GEORGIA INSTITUTE OF TECHNOLOGY ! *** WRITTEN BY ATHANASIOS NENES ! *** UPDATED BY CHRISTOS FOUNTOUKIS ! !======================================================================= ! SUBROUTINE CALCG1A USE ISRPIA IMPLICIT DOUBLE PRECISION (A-H, O-Z) ! INCLUDE 'ISRPIA.INC' DOUBLE PRECISION LAMDA, LAMDA1, LAMDA2, KAPA, KAPA1, KAPA2 ! ! *** CALCULATE NON VOLATILE SOLIDS *********************************** ! CNA2SO4 = 0.5*W(1) CNH42S4 = W(2) - CNA2SO4 ! ! *** CALCULATE VOLATILE SPECIES ************************************** ! ALF = W(3) - 2.0*CNH42S4 BET = W(5) GAM = W(4) ! RTSQ = R*TEMP*R*TEMP A1 = XK6/RTSQ A2 = XK10/RTSQ ! THETA1 = GAM - BET*(A2/A1) THETA2 = A2/A1 ! ! QUADRATIC EQUATION SOLUTION ! BB = (THETA1-ALF-BET*(ONE+THETA2))/(ONE+THETA2) CC = (ALF*BET-A1-BET*THETA1)/(ONE+THETA2) DD = BB*BB - 4.0D0*CC IF (DD.LT.ZERO) GOTO 100 ! Solve each reaction seperately ! ! TWO ROOTS FOR KAPA, CHECK AND SEE IF ANY VALID ! SQDD = SQRT(DD) KAPA1 = 0.5D0*(-BB+SQDD) KAPA2 = 0.5D0*(-BB-SQDD) LAMDA1 = THETA1 + THETA2*KAPA1 LAMDA2 = THETA1 + THETA2*KAPA2 ! IF (KAPA1.GE.ZERO .AND. LAMDA1.GE.ZERO) THEN IF (ALF-KAPA1-LAMDA1.GE.ZERO .AND. & BET-KAPA1.GE.ZERO .AND. GAM-LAMDA1.GE.ZERO) THEN KAPA = KAPA1 LAMDA= LAMDA1 GOTO 200 ENDIF ENDIF ! IF (KAPA2.GE.ZERO .AND. LAMDA2.GE.ZERO) THEN IF (ALF-KAPA2-LAMDA2.GE.ZERO .AND. & BET-KAPA2.GE.ZERO .AND. GAM-LAMDA2.GE.ZERO) THEN KAPA = KAPA2 LAMDA= LAMDA2 GOTO 200 ENDIF ENDIF ! ! SEPERATE SOLUTION OF NH4CL & NH4NO3 EQUILIBRIA ! 100 KAPA = ZERO LAMDA = ZERO DD1 = (ALF+BET)*(ALF+BET) - 4.0D0*(ALF*BET-A1) DD2 = (ALF+GAM)*(ALF+GAM) - 4.0D0*(ALF*GAM-A2) ! ! NH4CL EQUILIBRIUM ! IF (DD1.GE.ZERO) THEN SQDD1 = SQRT(DD1) KAPA1 = 0.5D0*(ALF+BET + SQDD1) KAPA2 = 0.5D0*(ALF+BET - SQDD1) ! IF (KAPA1.GE.ZERO .AND. KAPA1.LE.MIN(ALF,BET)) THEN KAPA = KAPA1 ELSE IF (KAPA2.GE.ZERO .AND. KAPA2.LE.MIN(ALF,BET)) THEN KAPA = KAPA2 ELSE KAPA = ZERO ENDIF ENDIF ! ! NH4NO3 EQUILIBRIUM ! IF (DD2.GE.ZERO) THEN SQDD2 = SQRT(DD2) LAMDA1= 0.5D0*(ALF+GAM + SQDD2) LAMDA2= 0.5D0*(ALF+GAM - SQDD2) ! IF (LAMDA1.GE.ZERO .AND. LAMDA1.LE.MIN(ALF,GAM)) THEN LAMDA = LAMDA1 ELSE IF (LAMDA2.GE.ZERO .AND. LAMDA2.LE.MIN(ALF,GAM)) THEN LAMDA = LAMDA2 ELSE LAMDA = ZERO ENDIF ENDIF ! ! IF BOTH KAPA, LAMDA ARE > 0, THEN APPLY EXISTANCE CRITERION ! IF (KAPA.GT.ZERO .AND. LAMDA.GT.ZERO) THEN IF (BET .LT. LAMDA/THETA1) THEN KAPA = ZERO ELSE LAMDA= ZERO ENDIF ENDIF ! ! *** CALCULATE COMPOSITION OF VOLATILE SPECIES *********************** ! 200 CONTINUE CNH4NO3 = LAMDA CNH4CL = KAPA ! GNH3 = MAX(ALF - KAPA - LAMDA, ZERO) GHNO3 = MAX(GAM - LAMDA, ZERO) GHCL = MAX(BET - KAPA, ZERO) ! RETURN ! ! *** END OF SUBROUTINE CALCG1A ***************************************** ! END SUBROUTINE CALCG1A !======================================================================= ! ! *** ISORROPIA CODE ! *** SUBROUTINE CALCH6 ! *** CASE H6 ! ! THE MAIN CHARACTERISTICS OF THIS REGIME ARE: ! 1. SULFATE POOR (SULRAT > 2.0) ; SODIUM RICH (SODRAT >= 2.0) ! 2. THERE IS BOTH A LIQUID & SOLID PHASE ! 3. SOLIDS POSSIBLE : (NH4)2SO4, NH4CL, NA2SO4 ! ! *** COPYRIGHT 1996-2006, UNIVERSITY OF MIAMI, CARNEGIE MELLON UNIVERSITY, ! *** GEORGIA INSTITUTE OF TECHNOLOGY ! *** WRITTEN BY ATHANASIOS NENES ! *** UPDATED BY CHRISTOS FOUNTOUKIS ! !======================================================================= ! SUBROUTINE CALCH6 USE ISRPIA USE SOLUT IMPLICIT DOUBLE PRECISION (A-H, O-Z) ! INCLUDE 'ISRPIA.INC' ! ! COMMON /SOLUT/ CHI1, CHI2, CHI3, CHI4, CHI5, CHI6, CHI7, CHI8, & ! PSI1, PSI2, PSI3, PSI4, PSI5, PSI6, PSI7, PSI8, & ! A1, A2, A3, A4, A5, A6, A7, A8 ! ! *** SETUP PARAMETERS ************************************************ ! CALAOU = .TRUE. CHI1 = W(2) ! CNA2SO4 CHI2 = ZERO ! CNH42S4 CHI3 = ZERO ! CNH4CL FRNA = MAX (W(1)-2.D0*CHI1, ZERO) CHI8 = MIN (FRNA, W(4)) ! CNANO3 CHI4 = W(3) ! NH3(g) CHI5 = MAX (W(4)-CHI8, ZERO) ! HNO3(g) CHI7 = MIN (MAX(FRNA-CHI8, ZERO), W(5)) ! CNACL CHI6 = MAX (W(5)-CHI7, ZERO) ! HCL(g) ! PSI6LO = TINY PSI6HI = CHI6-TINY ! MIN(CHI6-TINY, CHI4) ! ! *** INITIAL VALUES FOR BISECTION ************************************ ! X1 = PSI6LO Y1 = FUNCH6A (X1) IF (ABS(Y1).LE.EPS .OR. CHI6.LE.TINY) GOTO 50 ! ! *** ROOT TRACKING ; FOR THE RANGE OF HI AND LO ********************** ! DX = (PSI6HI-PSI6LO)/REAL(NDIV) DO 10 I=1,NDIV X2 = X1+DX Y2 = FUNCH6A (X2) IF (SIGN(1.d0,Y1)*SIGN(1.d0,Y2).LT.ZERO) GOTO 20 ! (Y1*Y2.LT.ZERO) X1 = X2 Y1 = Y2 10 CONTINUE ! ! *** NO SUBDIVISION WITH SOLUTION; IF ABS(Y2) 2.0) ; SODIUM RICH (SODRAT >= 2.0) ! 2. THERE IS BOTH A LIQUID & SOLID PHASE ! 3. SOLIDS POSSIBLE : (NH4)2SO4, NH4CL, NA2SO4 ! ! *** COPYRIGHT 1996-2006, UNIVERSITY OF MIAMI, CARNEGIE MELLON UNIVERSITY, ! *** GEORGIA INSTITUTE OF TECHNOLOGY ! *** WRITTEN BY ATHANASIOS NENES ! *** UPDATED BY CHRISTOS FOUNTOUKIS ! !======================================================================= ! DOUBLE PRECISION FUNCTION FUNCH6A (X) USE ISRPIA USE SOLUT IMPLICIT DOUBLE PRECISION (A-H, O-Z) ! INCLUDE 'ISRPIA.INC' ! ! COMMON /SOLUT/ CHI1, CHI2, CHI3, CHI4, CHI5, CHI6, CHI7, CHI8, & ! PSI1, PSI2, PSI3, PSI4, PSI5, PSI6, PSI7, PSI8, & ! A1, A2, A3, A4, A5, A6, A7, A8 ! ! *** SETUP PARAMETERS ************************************************ ! PSI6 = X PSI1 = CHI1 PSI2 = ZERO PSI3 = ZERO PSI7 = CHI7 PSI8 = CHI8 FRST = .TRUE. CALAIN = .TRUE. ! ! *** SOLVE EQUATIONS ; WITH ITERATIONS FOR ACTIVITY COEF. ************ ! DO 10 I=1,NSWEEP ! A1 = XK5 *(WATER/GAMA(2))**3.0 A4 = (XK2/XKW)*R*TEMP*(GAMA(10)/GAMA(5))**2.0 A5 = XK4 *R*TEMP*(WATER/GAMA(10))**2.0 A6 = XK3 *R*TEMP*(WATER/GAMA(11))**2.0 A7 = XK8 *(WATER/GAMA(1))**2.0 A8 = XK9 *(WATER/GAMA(3))**2.0 A9 = XK1*WATER/GAMA(7)*(GAMA(8)/GAMA(7))**2. ! ! CALCULATE DISSOCIATION QUANTITIES ! PSI5 = CHI5*(PSI6+PSI7) - A6/A5*PSI8*(CHI6-PSI6-PSI3) PSI5 = PSI5/(A6/A5*(CHI6-PSI6-PSI3) + PSI6 + PSI7) PSI5 = MAX(PSI5, TINY) ! IF (W(3).GT.TINY .AND. WATER.GT.TINY) THEN ! First try 3rd order soln BB =-(CHI4 + PSI6 + PSI5 + 1.d0/A4) CC = CHI4*(PSI5+PSI6) DD = BB*BB-4.d0*CC PSI4 =0.5d0*(-BB - SQRT(DD)) PSI4 = MIN(PSI4,CHI4) ELSE PSI4 = TINY ENDIF ! ! *** CALCULATE SPECIATION ******************************************** ! MOLAL (2) = PSI8 + PSI7 + 2.D0*PSI1 ! NAI MOLAL (3) = PSI4 ! NH4I MOLAL (4) = PSI6 + PSI7 ! CLI MOLAL (5) = PSI2 + PSI1 ! SO4I MOLAL (6) = ZERO ! HSO4I MOLAL (7) = PSI5 + PSI8 ! NO3I ! SMIN = 2.d0*MOLAL(5)+MOLAL(7)+MOLAL(4)-MOLAL(2)-MOLAL(3) CALL CALCPH (SMIN, HI, OHI) MOLAL (1) = HI ! GNH3 = MAX(CHI4 - PSI4, TINY) GHNO3 = MAX(CHI5 - PSI5, TINY) GHCL = MAX(CHI6 - PSI6, TINY) ! CNH42S4 = ZERO CNH4NO3 = ZERO CNACL = MAX(CHI7 - PSI7, ZERO) CNANO3 = MAX(CHI8 - PSI8, ZERO) CNA2SO4 = MAX(CHI1 - PSI1, ZERO) ! CALL CALCMR ! Water content ! ! *** CALCULATE ACTIVITIES OR TERMINATE INTERNAL LOOP ***************** ! IF (FRST.AND.CALAOU .OR. .NOT.FRST.AND.CALAIN) THEN CALL CALCACT ELSE GOTO 20 ENDIF 10 CONTINUE ! ! *** CALCULATE FUNCTION VALUE FOR OUTER LOOP *************************** ! 20 FUNCH6A = MOLAL(3)*MOLAL(4)/GHCL/GNH3/A6/A4 - ONE ! RETURN ! ! *** END OF FUNCTION FUNCH6A ******************************************* ! END FUNCTION FUNCH6A !======================================================================= ! ! *** ISORROPIA CODE ! *** SUBROUTINE CALCH5 ! *** CASE H5 ! ! THE MAIN CHARACTERISTICS OF THIS REGIME ARE: ! 1. SULFATE POOR (SULRAT > 2.0) ; SODIUM RICH (SODRAT >= 2.0) ! 2. THERE IS BOTH A LIQUID & SOLID PHASE ! 3. SOLIDS POSSIBLE : (NH4)2SO4, NH4CL, NA2SO4 ! ! *** COPYRIGHT 1996-2006, UNIVERSITY OF MIAMI, CARNEGIE MELLON UNIVERSITY, ! *** GEORGIA INSTITUTE OF TECHNOLOGY ! *** WRITTEN BY ATHANASIOS NENES ! *** UPDATED BY CHRISTOS FOUNTOUKIS ! !======================================================================= ! SUBROUTINE CALCH5 USE ISRPIA USE SOLUT IMPLICIT DOUBLE PRECISION (A-H, O-Z) ! INCLUDE 'ISRPIA.INC' ! ! COMMON /SOLUT/ CHI1, CHI2, CHI3, CHI4, CHI5, CHI6, CHI7, CHI8, & ! PSI1, PSI2, PSI3, PSI4, PSI5, PSI6, PSI7, PSI8, & ! A1, A2, A3, A4, A5, A6, A7, A8 ! ! *** REGIME DEPENDS ON THE EXISTANCE OF NITRATES *********************** ! IF (W(4).LE.TINY .AND. W(5).LE.TINY) THEN SCASE = 'H5' CALL CALCH1A SCASE = 'H5' RETURN ENDIF ! ! *** SETUP PARAMETERS ************************************************ ! CALAOU = .TRUE. CHI1 = W(2) ! CNA2SO4 CHI2 = ZERO ! CNH42S4 CHI3 = ZERO ! CNH4CL FRNA = MAX (W(1)-2.D0*CHI1, ZERO) CHI8 = MIN (FRNA, W(4)) ! CNANO3 CHI4 = W(3) ! NH3(g) CHI5 = MAX (W(4)-CHI8, ZERO) ! HNO3(g) CHI7 = MIN (MAX(FRNA-CHI8, ZERO), W(5)) ! CNACL CHI6 = MAX (W(5)-CHI7, ZERO) ! HCL(g) ! PSI6LO = TINY PSI6HI = CHI6-TINY ! MIN(CHI6-TINY, CHI4) ! ! *** INITIAL VALUES FOR BISECTION ************************************ ! X1 = PSI6LO Y1 = FUNCH5A (X1) IF (ABS(Y1).LE.EPS .OR. CHI6.LE.TINY) GOTO 50 ! ! *** ROOT TRACKING ; FOR THE RANGE OF HI AND LO ********************** ! DX = (PSI6HI-PSI6LO)/REAL(NDIV) DO 10 I=1,NDIV X2 = X1+DX Y2 = FUNCH5A (X2) IF (SIGN(1.d0,Y1)*SIGN(1.d0,Y2).LT.ZERO) GOTO 20 ! (Y1*Y2.LT.ZERO) X1 = X2 Y1 = Y2 10 CONTINUE ! ! *** NO SUBDIVISION WITH SOLUTION; IF ABS(Y2) 2.0) ; SODIUM RICH (SODRAT >= 2.0) ! 2. THERE IS BOTH A LIQUID & SOLID PHASE ! 3. SOLIDS POSSIBLE : NONE ! ! *** COPYRIGHT 1996-2006, UNIVERSITY OF MIAMI, CARNEGIE MELLON UNIVERSITY, ! *** GEORGIA INSTITUTE OF TECHNOLOGY ! *** WRITTEN BY ATHANASIOS NENES ! *** UPDATED BY CHRISTOS FOUNTOUKIS ! !======================================================================= ! DOUBLE PRECISION FUNCTION FUNCH5A (X) USE ISRPIA USE SOLUT IMPLICIT DOUBLE PRECISION (A-H, O-Z) ! INCLUDE 'ISRPIA.INC' ! ! COMMON /SOLUT/ CHI1, CHI2, CHI3, CHI4, CHI5, CHI6, CHI7, CHI8, & ! PSI1, PSI2, PSI3, PSI4, PSI5, PSI6, PSI7, PSI8, & ! A1, A2, A3, A4, A5, A6, A7, A8 ! ! *** SETUP PARAMETERS ************************************************ ! PSI6 = X PSI1 = CHI1 PSI2 = ZERO PSI3 = ZERO PSI7 = CHI7 PSI8 = CHI8 FRST = .TRUE. CALAIN = .TRUE. ! ! *** SOLVE EQUATIONS ; WITH ITERATIONS FOR ACTIVITY COEF. ************ ! DO 10 I=1,NSWEEP ! A1 = XK5 *(WATER/GAMA(2))**3.0 A4 = (XK2/XKW)*R*TEMP*(GAMA(10)/GAMA(5))**2.0 A5 = XK4 *R*TEMP*(WATER/GAMA(10))**2.0 A6 = XK3 *R*TEMP*(WATER/GAMA(11))**2.0 A7 = XK8 *(WATER/GAMA(1))**2.0 A8 = XK9 *(WATER/GAMA(3))**2.0 A9 = XK1*WATER/GAMA(7)*(GAMA(8)/GAMA(7))**2. ! ! CALCULATE DISSOCIATION QUANTITIES ! PSI5 = CHI5*(PSI6+PSI7) - A6/A5*PSI8*(CHI6-PSI6-PSI3) PSI5 = PSI5/(A6/A5*(CHI6-PSI6-PSI3) + PSI6 + PSI7) PSI5 = MAX(PSI5, TINY) ! IF (W(3).GT.TINY .AND. WATER.GT.TINY) THEN ! First try 3rd order soln BB =-(CHI4 + PSI6 + PSI5 + 1.d0/A4) CC = CHI4*(PSI5+PSI6) DD = BB*BB-4.d0*CC PSI4 =0.5d0*(-BB - SQRT(DD)) PSI4 = MIN(PSI4,CHI4) ELSE PSI4 = TINY ENDIF ! IF (CHI1.GT.TINY .AND. WATER.GT.TINY) THEN ! NA2SO4 DISSOLUTION AA = PSI7+PSI8 BB = AA*AA CC =-A1/4.D0 CALL POLY3 (AA, BB, CC, PSI1, ISLV) IF (ISLV.EQ.0) THEN PSI1 = MIN (PSI1, CHI1) ELSE PSI1 = ZERO ENDIF ENDIF ! ! *** CALCULATE SPECIATION ******************************************** ! MOLAL (2) = PSI8 + PSI7 + 2.D0*PSI1 ! NAI MOLAL (3) = PSI4 ! NH4I MOLAL (4) = PSI6 + PSI7 ! CLI MOLAL (5) = PSI2 + PSI1 ! SO4I MOLAL (6) = ZERO MOLAL (7) = PSI5 + PSI8 ! NO3I ! SMIN = 2.d0*MOLAL(5)+MOLAL(7)+MOLAL(4)-MOLAL(2)-MOLAL(3) CALL CALCPH (SMIN, HI, OHI) MOLAL (1) = HI ! GNH3 = MAX(CHI4 - PSI4, TINY) GHNO3 = MAX(CHI5 - PSI5, TINY) GHCL = MAX(CHI6 - PSI6, TINY) ! CNH42S4 = ZERO CNH4NO3 = ZERO CNACL = MAX(CHI7 - PSI7, ZERO) CNANO3 = MAX(CHI8 - PSI8, ZERO) CNA2SO4 = MAX(CHI1 - PSI1, ZERO) ! CALL CALCMR ! Water content ! ! *** CALCULATE ACTIVITIES OR TERMINATE INTERNAL LOOP ***************** ! IF (FRST.AND.CALAOU .OR. .NOT.FRST.AND.CALAIN) THEN CALL CALCACT ELSE GOTO 20 ENDIF 10 CONTINUE ! ! *** CALCULATE FUNCTION VALUE FOR OUTER LOOP *************************** ! 20 FUNCH5A = MOLAL(3)*MOLAL(4)/GHCL/GNH3/A6/A4 - ONE ! RETURN ! ! *** END OF FUNCTION FUNCH5A ******************************************* ! END FUNCTION FUNCH5A !======================================================================= ! ! *** ISORROPIA CODE ! *** SUBROUTINE CALCH4 ! *** CASE H4 ! ! THE MAIN CHARACTERISTICS OF THIS REGIME ARE: ! 1. SULFATE POOR (SULRAT > 2.0) ; SODIUM RICH (SODRAT >= 2.0) ! 2. THERE IS BOTH A LIQUID & SOLID PHASE ! 3. SOLIDS POSSIBLE : NA2SO4 ! ! *** COPYRIGHT 1996-2006, UNIVERSITY OF MIAMI, CARNEGIE MELLON UNIVERSITY, ! *** GEORGIA INSTITUTE OF TECHNOLOGY ! *** WRITTEN BY ATHANASIOS NENES ! *** UPDATED BY CHRISTOS FOUNTOUKIS ! !======================================================================= ! SUBROUTINE CALCH4 USE ISRPIA USE SOLUT IMPLICIT DOUBLE PRECISION (A-H, O-Z) ! INCLUDE 'ISRPIA.INC' ! ! COMMON /SOLUT/ CHI1, CHI2, CHI3, CHI4, CHI5, CHI6, CHI7, CHI8, & ! PSI1, PSI2, PSI3, PSI4, PSI5, PSI6, PSI7, PSI8, & ! A1, A2, A3, A4, A5, A6, A7, A8 ! ! *** REGIME DEPENDS ON THE EXISTANCE OF NITRATES *********************** ! IF (W(4).LE.TINY .AND. W(5).LE.TINY) THEN SCASE = 'H4' CALL CALCH1A SCASE = 'H4' RETURN ENDIF ! ! *** SETUP PARAMETERS ************************************************ ! CALAOU = .TRUE. CHI1 = W(2) ! CNA2SO4 CHI2 = ZERO ! CNH42S4 CHI3 = ZERO ! CNH4CL FRNA = MAX (W(1)-2.D0*CHI1, ZERO) CHI8 = MIN (FRNA, W(4)) ! CNANO3 CHI4 = W(3) ! NH3(g) CHI5 = MAX (W(4)-CHI8, ZERO) ! HNO3(g) CHI7 = MIN (MAX(FRNA-CHI8, ZERO), W(5)) ! CNACL CHI6 = MAX (W(5)-CHI7, ZERO) ! HCL(g) ! PSI6LO = TINY PSI6HI = CHI6-TINY ! MIN(CHI6-TINY, CHI4) ! ! *** INITIAL VALUES FOR BISECTION ************************************ ! X1 = PSI6LO Y1 = FUNCH4A (X1) IF (ABS(Y1).LE.EPS .OR. CHI6.LE.TINY) GOTO 50 ! ! *** ROOT TRACKING ; FOR THE RANGE OF HI AND LO ********************** ! DX = (PSI6HI-PSI6LO)/REAL(NDIV) DO 10 I=1,NDIV X2 = X1+DX Y2 = FUNCH4A (X2) IF (SIGN(1.d0,Y1)*SIGN(1.d0,Y2).LT.ZERO) GOTO 20 ! (Y1*Y2.LT.ZERO) X1 = X2 Y1 = Y2 10 CONTINUE ! ! *** NO SUBDIVISION WITH SOLUTION; IF ABS(Y2) 2.0) ; SODIUM RICH (SODRAT >= 2.0) ! 2. THERE IS BOTH A LIQUID & SOLID PHASE ! 3. SOLIDS POSSIBLE : (NH4)2SO4, NH4CL, NA2SO4 ! ! *** COPYRIGHT 1996-2006, UNIVERSITY OF MIAMI, CARNEGIE MELLON UNIVERSITY, ! *** GEORGIA INSTITUTE OF TECHNOLOGY ! *** WRITTEN BY ATHANASIOS NENES ! *** UPDATED BY CHRISTOS FOUNTOUKIS ! !======================================================================= ! DOUBLE PRECISION FUNCTION FUNCH4A (X) USE ISRPIA USE SOLUT IMPLICIT DOUBLE PRECISION (A-H, O-Z) ! INCLUDE 'ISRPIA.INC' ! ! COMMON /SOLUT/ CHI1, CHI2, CHI3, CHI4, CHI5, CHI6, CHI7, CHI8, & ! PSI1, PSI2, PSI3, PSI4, PSI5, PSI6, PSI7, PSI8, & ! A1, A2, A3, A4, A5, A6, A7, A8 ! ! *** SETUP PARAMETERS ************************************************ ! PSI6 = X PSI1 = CHI1 PSI2 = ZERO PSI3 = ZERO PSI7 = CHI7 PSI8 = CHI8 FRST = .TRUE. CALAIN = .TRUE. ! ! *** SOLVE EQUATIONS ; WITH ITERATIONS FOR ACTIVITY COEF. ************ ! DO 10 I=1,NSWEEP ! A1 = XK5 *(WATER/GAMA(2))**3.0 A4 = (XK2/XKW)*R*TEMP*(GAMA(10)/GAMA(5))**2.0 A5 = XK4 *R*TEMP*(WATER/GAMA(10))**2.0 A6 = XK3 *R*TEMP*(WATER/GAMA(11))**2.0 A7 = XK8 *(WATER/GAMA(1))**2.0 A8 = XK9 *(WATER/GAMA(3))**2.0 A9 = XK1*WATER/GAMA(7)*(GAMA(8)/GAMA(7))**2. ! ! CALCULATE DISSOCIATION QUANTITIES ! PSI5 = CHI5*(PSI6+PSI7) - A6/A5*PSI8*(CHI6-PSI6-PSI3) PSI5 = PSI5/(A6/A5*(CHI6-PSI6-PSI3) + PSI6 + PSI7) PSI5 = MAX(PSI5, TINY) ! IF (W(3).GT.TINY .AND. WATER.GT.TINY) THEN ! First try 3rd order soln BB =-(CHI4 + PSI6 + PSI5 + 1.d0/A4) CC = CHI4*(PSI5+PSI6) DD = BB*BB-4.d0*CC PSI4 =0.5d0*(-BB - SQRT(DD)) PSI4 = MIN(PSI4,CHI4) ELSE PSI4 = TINY ENDIF ! IF (CHI1.GT.TINY .AND. WATER.GT.TINY) THEN ! NA2SO4 DISSOLUTION AA = PSI7+PSI8 BB = AA*AA CC =-A1/4.D0 CALL POLY3 (AA, BB, CC, PSI1, ISLV) IF (ISLV.EQ.0) THEN PSI1 = MIN (PSI1, CHI1) ELSE PSI1 = ZERO ENDIF ENDIF ! ! *** CALCULATE SPECIATION ******************************************** ! MOLAL (2) = PSI8 + PSI7 + 2.D0*PSI1 ! NAI MOLAL (3) = PSI4 ! NH4I MOLAL (4) = PSI6 + PSI7 ! CLI MOLAL (5) = PSI2 + PSI1 ! SO4I MOLAL (6) = ZERO MOLAL (7) = PSI5 + PSI8 ! NO3I ! SMIN = 2.d0*MOLAL(5)+MOLAL(7)+MOLAL(4)-MOLAL(2)-MOLAL(3) CALL CALCPH (SMIN, HI, OHI) MOLAL (1) = HI ! ! *** CALCULATE GAS / SOLID SPECIES (LIQUID IN MOLAL ALREADY) ********* ! GNH3 = MAX(CHI4 - PSI4, TINY) GHNO3 = MAX(CHI5 - PSI5, TINY) GHCL = MAX(CHI6 - PSI6, TINY) ! CNH42S4 = ZERO CNH4NO3 = ZERO CNACL = MAX(CHI7 - PSI7, ZERO) CNANO3 = MAX(CHI8 - PSI8, ZERO) CNA2SO4 = MAX(CHI1 - PSI1, ZERO) ! ! *** NH4Cl(s) calculations ! A3 = XK6 /(R*TEMP*R*TEMP) DELT = MIN(GNH3, GHCL) BB = -(GNH3+GHCL) CC = GNH3*GHCL-A3 DD = BB*BB - 4.D0*CC PSI31 = 0.5D0*(-BB + SQRT(DD)) PSI32 = 0.5D0*(-BB - SQRT(DD)) IF (DELT-PSI31.GT.ZERO .AND. PSI31.GT.ZERO) THEN PSI3 = PSI31 ELSEIF (DELT-PSI32.GT.ZERO .AND. PSI32.GT.ZERO) THEN PSI3 = PSI32 ELSE PSI3 = ZERO ENDIF ! ! *** CALCULATE GAS / SOLID SPECIES (LIQUID IN MOLAL ALREADY) ********* ! GNH3 = MAX(GNH3 - PSI3, TINY) GHCL = MAX(GHCL - PSI3, TINY) CNH4CL = PSI3 ! CALL CALCMR ! Water content ! ! *** CALCULATE ACTIVITIES OR TERMINATE INTERNAL LOOP ***************** ! IF (FRST.AND.CALAOU .OR. .NOT.FRST.AND.CALAIN) THEN CALL CALCACT ELSE GOTO 20 ENDIF 10 CONTINUE ! ! *** CALCULATE FUNCTION VALUE FOR OUTER LOOP *************************** ! 20 FUNCH4A = MOLAL(3)*MOLAL(4)/GHCL/GNH3/A6/A4 - ONE ! RETURN ! ! *** END OF FUNCTION FUNCH4A ******************************************* ! END FUNCTION FUNCH4A !======================================================================= ! ! *** ISORROPIA CODE ! *** SUBROUTINE CALCH3 ! *** CASE H3 ! ! THE MAIN CHARACTERISTICS OF THIS REGIME ARE: ! 1. SULFATE POOR (SULRAT > 2.0) ; SODIUM RICH (SODRAT >= 2.0) ! 2. THERE IS BOTH A LIQUID & SOLID PHASE ! 3. SOLIDS POSSIBLE : NH4CL, NA2SO4 ! ! *** COPYRIGHT 1996-2006, UNIVERSITY OF MIAMI, CARNEGIE MELLON UNIVERSITY, ! *** GEORGIA INSTITUTE OF TECHNOLOGY ! *** WRITTEN BY ATHANASIOS NENES ! *** UPDATED BY CHRISTOS FOUNTOUKIS ! !======================================================================= ! SUBROUTINE CALCH3 USE ISRPIA USE SOLUT IMPLICIT DOUBLE PRECISION (A-H, O-Z) ! INCLUDE 'ISRPIA.INC' ! ! COMMON /SOLUT/ CHI1, CHI2, CHI3, CHI4, CHI5, CHI6, CHI7, CHI8, & ! PSI1, PSI2, PSI3, PSI4, PSI5, PSI6, PSI7, PSI8, & ! A1, A2, A3, A4, A5, A6, A7, A8 ! ! *** REGIME DEPENDS ON THE EXISTANCE OF NITRATES *********************** ! IF (W(4).LE.TINY) THEN ! NO3 NOT EXIST, WATER NOT POSSIBLE SCASE = 'H3' CALL CALCH1A SCASE = 'H3' RETURN ENDIF ! ! *** SETUP PARAMETERS ************************************************ ! CALAOU = .TRUE. CHI1 = W(2) ! CNA2SO4 CHI2 = ZERO ! CNH42S4 CHI3 = ZERO ! CNH4CL FRNA = MAX (W(1)-2.D0*CHI1, ZERO) CHI8 = MIN (FRNA, W(4)) ! CNANO3 CHI4 = W(3) ! NH3(g) CHI5 = MAX (W(4)-CHI8, ZERO) ! HNO3(g) CHI7 = MIN (MAX(FRNA-CHI8, ZERO), W(5)) ! CNACL CHI6 = MAX (W(5)-CHI7, ZERO) ! HCL(g) ! PSI6LO = TINY PSI6HI = CHI6-TINY ! MIN(CHI6-TINY, CHI4) ! ! *** INITIAL VALUES FOR BISECTION ************************************ ! X1 = PSI6LO Y1 = FUNCH3A (X1) IF (ABS(Y1).LE.EPS .OR. CHI6.LE.TINY) GOTO 50 ! ! *** ROOT TRACKING ; FOR THE RANGE OF HI AND LO ********************** ! DX = (PSI6HI-PSI6LO)/REAL(NDIV) DO 10 I=1,NDIV X2 = X1+DX Y2 = FUNCH3A (X2) IF (SIGN(1.d0,Y1)*SIGN(1.d0,Y2).LT.ZERO) GOTO 20 ! (Y1*Y2.LT.ZERO) X1 = X2 Y1 = Y2 10 CONTINUE ! ! *** NO SUBDIVISION WITH SOLUTION; IF ABS(Y2) 2.0) ; SODIUM RICH (SODRAT >= 2.0) ! 2. THERE IS BOTH A LIQUID & SOLID PHASE ! 3. SOLIDS POSSIBLE : (NH4)2SO4, NH4CL, NA2SO4 ! ! *** COPYRIGHT 1996-2006, UNIVERSITY OF MIAMI, CARNEGIE MELLON UNIVERSITY, ! *** GEORGIA INSTITUTE OF TECHNOLOGY ! *** WRITTEN BY ATHANASIOS NENES ! *** UPDATED BY CHRISTOS FOUNTOUKIS ! !======================================================================= ! DOUBLE PRECISION FUNCTION FUNCH3A (X) USE ISRPIA USE SOLUT IMPLICIT DOUBLE PRECISION (A-H, O-Z) ! INCLUDE 'ISRPIA.INC' ! ! COMMON /SOLUT/ CHI1, CHI2, CHI3, CHI4, CHI5, CHI6, CHI7, CHI8, & ! PSI1, PSI2, PSI3, PSI4, PSI5, PSI6, PSI7, PSI8, & ! A1, A2, A3, A4, A5, A6, A7, A8 ! ! *** SETUP PARAMETERS ************************************************ ! PSI6 = X PSI1 = CHI1 PSI2 = ZERO PSI3 = ZERO PSI7 = CHI7 PSI8 = CHI8 FRST = .TRUE. CALAIN = .TRUE. ! ! *** SOLVE EQUATIONS ; WITH ITERATIONS FOR ACTIVITY COEF. ************ ! DO 10 I=1,NSWEEP ! A1 = XK5 *(WATER/GAMA(2))**3.0 A4 = (XK2/XKW)*R*TEMP*(GAMA(10)/GAMA(5))**2.0 A5 = XK4 *R*TEMP*(WATER/GAMA(10))**2.0 A6 = XK3 *R*TEMP*(WATER/GAMA(11))**2.0 A7 = XK8 *(WATER/GAMA(1))**2.0 A8 = XK9 *(WATER/GAMA(3))**2.0 A9 = XK1*WATER/GAMA(7)*(GAMA(8)/GAMA(7))**2. ! ! CALCULATE DISSOCIATION QUANTITIES ! PSI5 = CHI5*(PSI6+PSI7) - A6/A5*PSI8*(CHI6-PSI6-PSI3) PSI5 = PSI5/(A6/A5*(CHI6-PSI6-PSI3) + PSI6 + PSI7) PSI5 = MAX(PSI5, TINY) ! IF (W(3).GT.TINY .AND. WATER.GT.TINY) THEN ! First try 3rd order soln BB =-(CHI4 + PSI6 + PSI5 + 1.d0/A4) CC = CHI4*(PSI5+PSI6) DD = BB*BB-4.d0*CC PSI4 =0.5d0*(-BB - SQRT(DD)) PSI4 = MIN(PSI4,CHI4) ELSE PSI4 = TINY ENDIF ! IF (CHI7.GT.TINY .AND. WATER.GT.TINY) THEN ! NACL DISSOLUTION DIAK = (PSI8-PSI6)**2.D0 + 4.D0*A7 PSI7 = 0.5D0*( -(PSI8+PSI6) + SQRT(DIAK) ) PSI7 = MAX(MIN(PSI7, CHI7), ZERO) ENDIF ! IF (CHI1.GT.TINY .AND. WATER.GT.TINY) THEN ! NA2SO4 DISSOLUTION AA = PSI7+PSI8 BB = AA*AA CC =-A1/4.D0 CALL POLY3 (AA, BB, CC, PSI1, ISLV) IF (ISLV.EQ.0) THEN PSI1 = MIN (PSI1, CHI1) ELSE PSI1 = ZERO ENDIF ENDIF ! ! *** CALCULATE SPECIATION ******************************************** ! MOLAL (2) = PSI8 + PSI7 + 2.D0*PSI1 ! NAI MOLAL (3) = PSI4 ! NH4I MOLAL (4) = PSI6 + PSI7 ! CLI MOLAL (5) = PSI2 + PSI1 ! SO4I MOLAL (6) = ZERO MOLAL (7) = PSI5 + PSI8 ! NO3I ! SMIN = 2.d0*MOLAL(5)+MOLAL(7)+MOLAL(4)-MOLAL(2)-MOLAL(3) CALL CALCPH (SMIN, HI, OHI) MOLAL (1) = HI ! ! *** CALCULATE GAS / SOLID SPECIES (LIQUID IN MOLAL ALREADY) ********* ! GNH3 = MAX(CHI4 - PSI4, TINY) GHNO3 = MAX(CHI5 - PSI5, TINY) GHCL = MAX(CHI6 - PSI6, TINY) ! CNH42S4 = ZERO CNH4NO3 = ZERO CNACL = MAX(CHI7 - PSI7, ZERO) CNANO3 = MAX(CHI8 - PSI8, ZERO) CNA2SO4 = MAX(CHI1 - PSI1, ZERO) ! ! *** NH4Cl(s) calculations ! A3 = XK6 /(R*TEMP*R*TEMP) DELT = MIN(GNH3, GHCL) BB = -(GNH3+GHCL) CC = GNH3*GHCL-A3 DD = BB*BB - 4.D0*CC PSI31 = 0.5D0*(-BB + SQRT(DD)) PSI32 = 0.5D0*(-BB - SQRT(DD)) IF (DELT-PSI31.GT.ZERO .AND. PSI31.GT.ZERO) THEN PSI3 = PSI31 ELSEIF (DELT-PSI32.GT.ZERO .AND. PSI32.GT.ZERO) THEN PSI3 = PSI32 ELSE PSI3 = ZERO ENDIF ! ! *** CALCULATE GAS / SOLID SPECIES (LIQUID IN MOLAL ALREADY) ********* ! GNH3 = MAX(GNH3 - PSI3, TINY) GHCL = MAX(GHCL - PSI3, TINY) CNH4CL = PSI3 ! CALL CALCMR ! Water content ! ! *** CALCULATE ACTIVITIES OR TERMINATE INTERNAL LOOP ***************** ! IF (FRST.AND.CALAOU .OR. .NOT.FRST.AND.CALAIN) THEN CALL CALCACT ELSE GOTO 20 ENDIF 10 CONTINUE ! ! *** CALCULATE FUNCTION VALUE FOR OUTER LOOP *************************** ! 20 FUNCH3A = MOLAL(3)*MOLAL(4)/GHCL/GNH3/A6/A4 - ONE ! RETURN ! ! *** END OF FUNCTION FUNCH3A ******************************************* ! END FUNCTION FUNCH3A !======================================================================= ! ! *** ISORROPIA CODE ! *** SUBROUTINE CALCH2 ! *** CASE H2 ! ! THE MAIN CHARACTERISTICS OF THIS REGIME ARE: ! 1. SULFATE POOR (SULRAT > 2.0) ; SODIUM RICH (SODRAT >= 2.0) ! 2. SOLID & LIQUID AEROSOL POSSIBLE ! 3. SOLIDS POSSIBLE : NH4Cl, NA2SO4, NANO3, NACL ! ! THERE ARE THREE REGIMES IN THIS CASE: ! 1. NH4NO3(s) POSSIBLE. LIQUID & SOLID AEROSOL (SUBROUTINE CALCH2A) ! 2. NH4NO3(s) NOT POSSIBLE, AND RH < MDRH. SOLID AEROSOL ONLY ! 3. NH4NO3(s) NOT POSSIBLE, AND RH >= MDRH. (MDRH REGION) ! ! REGIMES 2. AND 3. ARE CONSIDERED TO BE THE SAME AS CASES H1A, H2B ! RESPECTIVELY (BECAUSE MDRH POINTS COINCIDE). ! ! *** COPYRIGHT 1996-2006, UNIVERSITY OF MIAMI, CARNEGIE MELLON UNIVERSITY, ! *** GEORGIA INSTITUTE OF TECHNOLOGY ! *** WRITTEN BY ATHANASIOS NENES ! *** UPDATED BY CHRISTOS FOUNTOUKIS ! !======================================================================= ! SUBROUTINE CALCH2 USE ISRPIA IMPLICIT DOUBLE PRECISION (A-H, O-Z) ! INCLUDE 'ISRPIA.INC' EXTERNAL CALCH1A, CALCH3 ! ! *** REGIME DEPENDS ON THE EXISTANCE OF NITRATES *********************** ! IF (W(4).GT.TINY) THEN ! NO3 EXISTS, WATER POSSIBLE SCASE = 'H2 ; SUBCASE 1' CALL CALCH2A SCASE = 'H2 ; SUBCASE 1' ELSE ! NO3 NON EXISTANT, WATER NOT POSSIBLE SCASE = 'H2 ; SUBCASE 1' CALL CALCH1A SCASE = 'H2 ; SUBCASE 1' ENDIF ! IF (WATER.LE.TINY .AND. RH.LT.DRMH2) THEN ! DRY AEROSOL SCASE = 'H2 ; SUBCASE 2' ! ELSEIF (WATER.LE.TINY .AND. RH.GE.DRMH2) THEN ! MDRH OF H2 SCASE = 'H2 ; SUBCASE 3' CALL CALCMDRH (RH, DRMH2, DRNANO3, CALCH1A, CALCH3) SCASE = 'H2 ; SUBCASE 3' ENDIF ! RETURN ! ! *** END OF SUBROUTINE CALCH2 ****************************************** ! END SUBROUTINE CALCH2 !======================================================================= ! ! *** ISORROPIA CODE ! *** SUBROUTINE CALCH2A ! *** CASE H2 ; SUBCASE 1 ! ! THE MAIN CHARACTERISTICS OF THIS REGIME ARE: ! 1. SULFATE POOR (SULRAT > 2.0) ; SODIUM RICH (SODRAT >= 2.0) ! 2. THERE IS BOTH A LIQUID & SOLID PHASE ! 3. SOLIDS POSSIBLE : NH4CL, NA2SO4 ! ! *** COPYRIGHT 1996-2006, UNIVERSITY OF MIAMI, CARNEGIE MELLON UNIVERSITY, ! *** GEORGIA INSTITUTE OF TECHNOLOGY ! *** WRITTEN BY ATHANASIOS NENES ! *** UPDATED BY CHRISTOS FOUNTOUKIS ! !======================================================================= ! SUBROUTINE CALCH2A USE ISRPIA USE SOLUT IMPLICIT DOUBLE PRECISION (A-H, O-Z) ! INCLUDE 'ISRPIA.INC' ! ! COMMON /SOLUT/ CHI1, CHI2, CHI3, CHI4, CHI5, CHI6, CHI7, CHI8, & ! PSI1, PSI2, PSI3, PSI4, PSI5, PSI6, PSI7, PSI8, & ! A1, A2, A3, A4, A5, A6, A7, A8 ! ! *** SETUP PARAMETERS ************************************************ ! CALAOU = .TRUE. CHI1 = W(2) ! CNA2SO4 CHI2 = ZERO ! CNH42S4 CHI3 = ZERO ! CNH4CL FRNA = MAX (W(1)-2.D0*CHI1, ZERO) CHI8 = MIN (FRNA, W(4)) ! CNANO3 CHI4 = W(3) ! NH3(g) CHI5 = MAX (W(4)-CHI8, ZERO) ! HNO3(g) CHI7 = MIN (MAX(FRNA-CHI8, ZERO), W(5)) ! CNACL CHI6 = MAX (W(5)-CHI7, ZERO) ! HCL(g) ! PSI6LO = TINY PSI6HI = CHI6-TINY ! MIN(CHI6-TINY, CHI4) ! ! *** INITIAL VALUES FOR BISECTION ************************************ ! X1 = PSI6LO Y1 = FUNCH2A (X1) IF (ABS(Y1).LE.EPS .OR. CHI6.LE.TINY) GOTO 50 ! ! *** ROOT TRACKING ; FOR THE RANGE OF HI AND LO ********************** ! DX = (PSI6HI-PSI6LO)/REAL(NDIV) DO 10 I=1,NDIV X2 = X1+DX Y2 = FUNCH2A (X2) IF (SIGN(1.d0,Y1)*SIGN(1.d0,Y2).LT.ZERO) GOTO 20 ! (Y1*Y2.LT.ZERO) X1 = X2 Y1 = Y2 10 CONTINUE ! ! *** NO SUBDIVISION WITH SOLUTION; IF ABS(Y2) 2.0) ; SODIUM RICH (SODRAT >= 2.0) ! 2. THERE IS BOTH A LIQUID & SOLID PHASE ! 3. SOLIDS POSSIBLE : (NH4)2SO4, NH4CL, NA2SO4 ! ! *** COPYRIGHT 1996-2006, UNIVERSITY OF MIAMI, CARNEGIE MELLON UNIVERSITY, ! *** GEORGIA INSTITUTE OF TECHNOLOGY ! *** WRITTEN BY ATHANASIOS NENES ! *** UPDATED BY CHRISTOS FOUNTOUKIS ! !======================================================================= ! DOUBLE PRECISION FUNCTION FUNCH2A (X) USE ISRPIA USE SOLUT IMPLICIT DOUBLE PRECISION (A-H, O-Z) ! INCLUDE 'ISRPIA.INC' ! ! COMMON /SOLUT/ CHI1, CHI2, CHI3, CHI4, CHI5, CHI6, CHI7, CHI8, & ! PSI1, PSI2, PSI3, PSI4, PSI5, PSI6, PSI7, PSI8, & ! A1, A2, A3, A4, A5, A6, A7, A8 ! ! *** SETUP PARAMETERS ************************************************ ! PSI6 = X PSI1 = CHI1 PSI2 = ZERO PSI3 = ZERO PSI7 = CHI7 PSI8 = CHI8 FRST = .TRUE. CALAIN = .TRUE. ! ! *** SOLVE EQUATIONS ; WITH ITERATIONS FOR ACTIVITY COEF. ************ ! DO 10 I=1,NSWEEP ! A1 = XK5 *(WATER/GAMA(2))**3.0 A4 = (XK2/XKW)*R*TEMP*(GAMA(10)/GAMA(5))**2.0 A5 = XK4 *R*TEMP*(WATER/GAMA(10))**2.0 A6 = XK3 *R*TEMP*(WATER/GAMA(11))**2.0 A7 = XK8 *(WATER/GAMA(1))**2.0 A8 = XK9 *(WATER/GAMA(3))**2.0 A64 = (XK3*XK2/XKW)*(GAMA(10)/GAMA(5)/GAMA(11))**2.0 A64 = A64*(R*TEMP*WATER)**2.0 A9 = XK1*WATER/GAMA(7)*(GAMA(8)/GAMA(7))**2. ! ! CALCULATE DISSOCIATION QUANTITIES ! PSI5 = CHI5*(PSI6+PSI7) - A6/A5*PSI8*(CHI6-PSI6-PSI3) PSI5 = PSI5/(A6/A5*(CHI6-PSI6-PSI3) + PSI6 + PSI7) PSI5 = MAX(PSI5, TINY) ! IF (W(3).GT.TINY .AND. WATER.GT.TINY) THEN ! First try 3rd order soln BB =-(CHI4 + PSI6 + PSI5 + 1.d0/A4) CC = CHI4*(PSI5+PSI6) DD = BB*BB-4.d0*CC PSI4 =0.5d0*(-BB - SQRT(DD)) PSI4 = MIN(PSI4,CHI4) ELSE PSI4 = TINY ENDIF ! IF (CHI7.GT.TINY .AND. WATER.GT.TINY) THEN ! NACL DISSOLUTION DIAK = (PSI8-PSI6)**2.D0 + 4.D0*A7 PSI7 = 0.5D0*( -(PSI8+PSI6) + SQRT(DIAK) ) PSI7 = MAX(MIN(PSI7, CHI7), ZERO) ENDIF ! IF (CHI8.GT.TINY .AND. WATER.GT.TINY) THEN ! NANO3 DISSOLUTION DIAK = (PSI7-PSI5)**2.D0 + 4.D0*A8 PSI8 = 0.5D0*( -(PSI7+PSI5) + SQRT(DIAK) ) PSI8 = MAX(MIN(PSI8, CHI8), ZERO) ENDIF ! IF (CHI1.GT.TINY .AND. WATER.GT.TINY) THEN ! NA2SO4 DISSOLUTION AA = PSI7+PSI8 BB = AA*AA CC =-A1/4.D0 CALL POLY3 (AA, BB, CC, PSI1, ISLV) IF (ISLV.EQ.0) THEN PSI1 = MIN (PSI1, CHI1) ELSE PSI1 = ZERO ENDIF ENDIF ! ! *** CALCULATE SPECIATION ******************************************** ! MOLAL (2) = PSI8 + PSI7 + 2.D0*PSI1 ! NAI MOLAL (3) = PSI4 ! NH4I MOLAL (4) = PSI6 + PSI7 ! CLI MOLAL (5) = PSI2 + PSI1 ! SO4I MOLAL (6) = ZERO ! HSO4I MOLAL (7) = PSI5 + PSI8 ! NO3I ! SMIN = 2.d0*MOLAL(5)+MOLAL(7)+MOLAL(4)-MOLAL(2)-MOLAL(3) CALL CALCPH (SMIN, HI, OHI) MOLAL (1) = HI ! GNH3 = MAX(CHI4 - PSI4, TINY) GHNO3 = MAX(CHI5 - PSI5, TINY) GHCL = MAX(CHI6 - PSI6, TINY) ! CNH42S4 = ZERO CNH4NO3 = ZERO CNACL = MAX(CHI7 - PSI7, ZERO) CNANO3 = MAX(CHI8 - PSI8, ZERO) CNA2SO4 = MAX(CHI1 - PSI1, ZERO) ! ! *** NH4Cl(s) calculations ! A3 = XK6 /(R*TEMP*R*TEMP) DELT = MIN(GNH3, GHCL) BB = -(GNH3+GHCL) CC = GNH3*GHCL-A3 DD = BB*BB - 4.D0*CC PSI31 = 0.5D0*(-BB + SQRT(DD)) PSI32 = 0.5D0*(-BB - SQRT(DD)) IF (DELT-PSI31.GT.ZERO .AND. PSI31.GT.ZERO) THEN PSI3 = PSI31 ELSEIF (DELT-PSI32.GT.ZERO .AND. PSI32.GT.ZERO) THEN PSI3 = PSI32 ELSE PSI3 = ZERO ENDIF ! ! *** CALCULATE GAS / SOLID SPECIES (LIQUID IN MOLAL ALREADY) ********* ! GNH3 = MAX(GNH3 - PSI3, TINY) GHCL = MAX(GHCL - PSI3, TINY) CNH4CL = PSI3 ! CALL CALCMR ! Water content ! ! *** CALCULATE ACTIVITIES OR TERMINATE INTERNAL LOOP ***************** ! IF (FRST.AND.CALAOU .OR. .NOT.FRST.AND.CALAIN) THEN CALL CALCACT ELSE GOTO 20 ENDIF 10 CONTINUE ! ! *** CALCULATE FUNCTION VALUE FOR OUTER LOOP *************************** ! 20 FUNCH2A = MOLAL(3)*MOLAL(4)/GHCL/GNH3/A64 - ONE ! RETURN ! ! *** END OF FUNCTION FUNCH2A ******************************************* ! END FUNCTION FUNCH2A !======================================================================= ! ! *** ISORROPIA CODE ! *** SUBROUTINE CALCH1 ! *** CASE H1 ! ! THE MAIN CHARACTERISTICS OF THIS REGIME ARE: ! 1. SULFATE POOR (SULRAT > 2.0) ; SODIUM RICH (SODRAT >= 2.0) ! 2. SOLID AEROSOL ONLY ! 3. SOLIDS POSSIBLE : NH4NO3, NH4CL, NA2SO4 ! ! THERE ARE TWO POSSIBLE REGIMES HERE, DEPENDING ON RELATIVE HUMIDITY: ! 1. WHEN RH >= MDRH ; LIQUID PHASE POSSIBLE (MDRH REGION) ! 2. WHEN RH < MDRH ; ONLY SOLID PHASE POSSIBLE (SUBROUTINE CALCH1A) ! ! *** COPYRIGHT 1996-2006, UNIVERSITY OF MIAMI, CARNEGIE MELLON UNIVERSITY, ! *** GEORGIA INSTITUTE OF TECHNOLOGY ! *** WRITTEN BY ATHANASIOS NENES ! *** UPDATED BY CHRISTOS FOUNTOUKIS ! !======================================================================= ! SUBROUTINE CALCH1 USE ISRPIA IMPLICIT DOUBLE PRECISION (A-H, O-Z) ! INCLUDE 'ISRPIA.INC' EXTERNAL CALCH1A, CALCH2A ! ! *** REGIME DEPENDS UPON THE AMBIENT RELATIVE HUMIDITY ***************** ! IF (RH.LT.DRMH1) THEN SCASE = 'H1 ; SUBCASE 1' CALL CALCH1A ! SOLID PHASE ONLY POSSIBLE SCASE = 'H1 ; SUBCASE 1' ELSE SCASE = 'H1 ; SUBCASE 2' ! LIQUID & SOLID PHASE POSSIBLE CALL CALCMDRH (RH, DRMH1, DRNH4NO3, CALCH1A, CALCH2A) SCASE = 'H1 ; SUBCASE 2' ENDIF ! RETURN ! ! *** END OF SUBROUTINE CALCH1 ****************************************** ! END SUBROUTINE CALCH1 !======================================================================= ! ! *** ISORROPIA CODE ! *** SUBROUTINE CALCH1A ! *** CASE H1 ; SUBCASE 1 ! ! THE MAIN CHARACTERISTICS OF THIS REGIME ARE: ! 1. SULFATE POOR (SULRAT > 2.0) ; SODIUM RICH (SODRAT >= 2.0) ! 2. SOLID AEROSOL ONLY ! 3. SOLIDS POSSIBLE : NH4NO3, NH4CL, NANO3, NA2SO4 ! ! *** COPYRIGHT 1996-2006, UNIVERSITY OF MIAMI, CARNEGIE MELLON UNIVERSITY, ! *** GEORGIA INSTITUTE OF TECHNOLOGY ! *** WRITTEN BY ATHANASIOS NENES ! *** UPDATED BY CHRISTOS FOUNTOUKIS ! !======================================================================= ! SUBROUTINE CALCH1A USE ISRPIA IMPLICIT DOUBLE PRECISION (A-H, O-Z) ! INCLUDE 'ISRPIA.INC' DOUBLE PRECISION LAMDA, LAMDA1, LAMDA2, KAPA, KAPA1, KAPA2, NAFR, & NO3FR ! ! *** CALCULATE NON VOLATILE SOLIDS *********************************** ! CNA2SO4 = W(2) CNH42S4 = ZERO NAFR = MAX (W(1)-2*CNA2SO4, ZERO) CNANO3 = MIN (NAFR, W(4)) NO3FR = MAX (W(4)-CNANO3, ZERO) CNACL = MIN (MAX(NAFR-CNANO3, ZERO), W(5)) CLFR = MAX (W(5)-CNACL, ZERO) ! ! *** CALCULATE VOLATILE SPECIES ************************************** ! ALF = W(3) ! FREE NH3 BET = CLFR ! FREE CL GAM = NO3FR ! FREE NO3 ! RTSQ = R*TEMP*R*TEMP A1 = XK6/RTSQ A2 = XK10/RTSQ ! THETA1 = GAM - BET*(A2/A1) THETA2 = A2/A1 ! ! QUADRATIC EQUATION SOLUTION ! BB = (THETA1-ALF-BET*(ONE+THETA2))/(ONE+THETA2) CC = (ALF*BET-A1-BET*THETA1)/(ONE+THETA2) DD = BB*BB - 4.0D0*CC IF (DD.LT.ZERO) GOTO 100 ! Solve each reaction seperately ! ! TWO ROOTS FOR KAPA, CHECK AND SEE IF ANY VALID ! SQDD = SQRT(DD) KAPA1 = 0.5D0*(-BB+SQDD) KAPA2 = 0.5D0*(-BB-SQDD) LAMDA1 = THETA1 + THETA2*KAPA1 LAMDA2 = THETA1 + THETA2*KAPA2 ! IF (KAPA1.GE.ZERO .AND. LAMDA1.GE.ZERO) THEN IF (ALF-KAPA1-LAMDA1.GE.ZERO .AND. & BET-KAPA1.GE.ZERO .AND. GAM-LAMDA1.GE.ZERO) THEN KAPA = KAPA1 LAMDA= LAMDA1 GOTO 200 ENDIF ENDIF ! IF (KAPA2.GE.ZERO .AND. LAMDA2.GE.ZERO) THEN IF (ALF-KAPA2-LAMDA2.GE.ZERO .AND. & BET-KAPA2.GE.ZERO .AND. GAM-LAMDA2.GE.ZERO) THEN KAPA = KAPA2 LAMDA= LAMDA2 GOTO 200 ENDIF ENDIF ! ! SEPERATE SOLUTION OF NH4CL & NH4NO3 EQUILIBRIA ! 100 KAPA = ZERO LAMDA = ZERO DD1 = (ALF+BET)*(ALF+BET) - 4.0D0*(ALF*BET-A1) DD2 = (ALF+GAM)*(ALF+GAM) - 4.0D0*(ALF*GAM-A2) ! ! NH4CL EQUILIBRIUM ! IF (DD1.GE.ZERO) THEN SQDD1 = SQRT(DD1) KAPA1 = 0.5D0*(ALF+BET + SQDD1) KAPA2 = 0.5D0*(ALF+BET - SQDD1) ! IF (KAPA1.GE.ZERO .AND. KAPA1.LE.MIN(ALF,BET)) THEN KAPA = KAPA1 ELSE IF (KAPA2.GE.ZERO .AND. KAPA2.LE.MIN(ALF,BET)) THEN KAPA = KAPA2 ELSE KAPA = ZERO ENDIF ENDIF ! ! NH4NO3 EQUILIBRIUM ! IF (DD2.GE.ZERO) THEN SQDD2 = SQRT(DD2) LAMDA1= 0.5D0*(ALF+GAM + SQDD2) LAMDA2= 0.5D0*(ALF+GAM - SQDD2) ! IF (LAMDA1.GE.ZERO .AND. LAMDA1.LE.MIN(ALF,GAM)) THEN LAMDA = LAMDA1 ELSE IF (LAMDA2.GE.ZERO .AND. LAMDA2.LE.MIN(ALF,GAM)) THEN LAMDA = LAMDA2 ELSE LAMDA = ZERO ENDIF ENDIF ! ! IF BOTH KAPA, LAMDA ARE > 0, THEN APPLY EXISTANCE CRITERION ! IF (KAPA.GT.ZERO .AND. LAMDA.GT.ZERO) THEN IF (BET .LT. LAMDA/THETA1) THEN KAPA = ZERO ELSE LAMDA= ZERO ENDIF ENDIF ! ! *** CALCULATE COMPOSITION OF VOLATILE SPECIES *********************** ! 200 CONTINUE CNH4NO3 = LAMDA CNH4CL = KAPA ! GNH3 = ALF - KAPA - LAMDA GHNO3 = GAM - LAMDA GHCL = BET - KAPA ! RETURN ! ! *** END OF SUBROUTINE CALCH1A ***************************************** ! END SUBROUTINE CALCH1A !======================================================================= ! ! *** ISORROPIA CODE ! *** SUBROUTINE CALCI6 ! *** CASE I6 ! ! THE MAIN CHARACTERISTICS OF THIS REGIME ARE: ! 1. SULFATE RICH, NO FREE ACID (1.0 <= SULRAT < 2.0) ! 2. SOLID & LIQUID AEROSOL POSSIBLE ! 3. SOLIDS POSSIBLE : (NH4)2SO4, NA2SO4 ! ! *** COPYRIGHT 1996-2006, UNIVERSITY OF MIAMI, CARNEGIE MELLON UNIVERSITY, ! *** GEORGIA INSTITUTE OF TECHNOLOGY ! *** WRITTEN BY ATHANASIOS NENES ! *** UPDATED BY CHRISTOS FOUNTOUKIS ! !======================================================================= ! SUBROUTINE CALCI6 USE ISRPIA USE SOLUT IMPLICIT DOUBLE PRECISION (A-H, O-Z) ! INCLUDE 'ISRPIA.INC' ! COMMON /SOLUT/ CHI1, CHI2, CHI3, CHI4, CHI5, CHI6, CHI7, CHI8, & ! PSI1, PSI2, PSI3, PSI4, PSI5, PSI6, PSI7, PSI8, & ! A1, A2, A3, A4, A5, A6, A7, A8 ! ! *** FIND DRY COMPOSITION ********************************************** ! CALL CALCI1A ! ! *** SETUP PARAMETERS ************************************************ ! CHI1 = CNH4HS4 ! Save from CALCI1 run CHI2 = CLC CHI3 = CNAHSO4 CHI4 = CNA2SO4 CHI5 = CNH42S4 ! PSI1 = CNH4HS4 ! ASSIGN INITIAL PSI's PSI2 = CLC PSI3 = CNAHSO4 PSI4 = CNA2SO4 PSI5 = CNH42S4 ! CALAOU = .TRUE. ! Outer loop activity calculation flag FRST = .TRUE. CALAIN = .TRUE. ! ! *** SOLVE EQUATIONS ; WITH ITERATIONS FOR ACTIVITY COEF. ************ ! DO 10 I=1,NSWEEP ! A6 = XK1 *WATER/GAMA(7)*(GAMA(8)/GAMA(7))**2. ! ! CALCULATE DISSOCIATION QUANTITIES ! BB = PSI2 + PSI4 + PSI5 + A6 ! PSI6 CC =-A6*(PSI2 + PSI3 + PSI1) DD = BB*BB - 4.D0*CC PSI6 = 0.5D0*(-BB + SQRT(DD)) ! ! *** CALCULATE SPECIATION ******************************************** ! MOLAL (1) = PSI6 ! HI MOLAL (2) = 2.D0*PSI4 + PSI3 ! NAI MOLAL (3) = 3.D0*PSI2 + 2.D0*PSI5 + PSI1 ! NH4I MOLAL (5) = PSI2 + PSI4 + PSI5 + PSI6 ! SO4I MOLAL (6) = PSI2 + PSI3 + PSI1 - PSI6 ! HSO4I CLC = ZERO CNAHSO4 = ZERO CNA2SO4 = CHI4 - PSI4 CNH42S4 = ZERO CNH4HS4 = ZERO CALL CALCMR ! Water content ! ! *** CALCULATE ACTIVITIES OR TERMINATE INTERNAL LOOP ***************** ! IF (FRST.AND.CALAOU .OR. .NOT.FRST.AND.CALAIN) THEN CALL CALCACT ELSE GOTO 20 ENDIF 10 CONTINUE ! 20 RETURN ! ! *** END OF SUBROUTINE CALCI6 ***************************************** ! END SUBROUTINE CALCI6 !======================================================================= ! ! *** ISORROPIA CODE ! *** SUBROUTINE CALCI5 ! *** CASE I5 ! ! THE MAIN CHARACTERISTICS OF THIS REGIME ARE: ! 1. SULFATE RICH, NO FREE ACID (1.0 <= SULRAT < 2.0) ! 2. SOLID & LIQUID AEROSOL POSSIBLE ! 3. SOLIDS POSSIBLE : (NH4)2SO4, NA2SO4 ! ! *** COPYRIGHT 1996-2006, UNIVERSITY OF MIAMI, CARNEGIE MELLON UNIVERSITY, ! *** GEORGIA INSTITUTE OF TECHNOLOGY ! *** WRITTEN BY ATHANASIOS NENES ! *** UPDATED BY CHRISTOS FOUNTOUKIS ! !======================================================================= ! SUBROUTINE CALCI5 USE ISRPIA USE SOLUT IMPLICIT DOUBLE PRECISION (A-H, O-Z) ! INCLUDE 'ISRPIA.INC' ! COMMON /SOLUT/ CHI1, CHI2, CHI3, CHI4, CHI5, CHI6, CHI7, CHI8, & ! PSI1, PSI2, PSI3, PSI4, PSI5, PSI6, PSI7, PSI8, & ! A1, A2, A3, A4, A5, A6, A7, A8 ! ! *** FIND DRY COMPOSITION ********************************************** ! CALL CALCI1A ! ! *** SETUP PARAMETERS ************************************************ ! CHI1 = CNH4HS4 ! Save from CALCI1 run CHI2 = CLC CHI3 = CNAHSO4 CHI4 = CNA2SO4 CHI5 = CNH42S4 ! PSI1 = CNH4HS4 ! ASSIGN INITIAL PSI's PSI2 = CLC PSI3 = CNAHSO4 PSI4 = ZERO PSI5 = CNH42S4 ! CALAOU =.TRUE. ! Outer loop activity calculation flag PSI4LO = ZERO ! Low limit PSI4HI = CHI4 ! High limit ! ! *** IF NA2SO4(S) =0, CALL FUNCI5B FOR Y4=0 *************************** ! IF (CHI4.LE.TINY) THEN Y1 = FUNCI5A (ZERO) GOTO 50 ENDIF ! ! *** INITIAL VALUES FOR BISECTION ************************************ ! X1 = PSI4HI Y1 = FUNCI5A (X1) YHI= Y1 ! Save Y-value at HI position ! ! *** YHI < 0.0 THE SOLUTION IS ALWAYS UNDERSATURATED WITH NA2SO4 ** ! IF (ABS(Y1).LE.EPS .OR. YHI.LT.ZERO) GOTO 50 ! ! *** ROOT TRACKING ; FOR THE RANGE OF HI AND LO ********************** ! DX = (PSI4HI-PSI4LO)/REAL(NDIV) DO 10 I=1,NDIV X2 = X1-DX Y2 = FUNCI5A (X2) IF (SIGN(1.d0,Y1)*SIGN(1.d0,Y2).LT.ZERO) GOTO 20 ! (Y1*Y2.LT.ZERO) X1 = X2 Y1 = Y2 10 CONTINUE ! ! *** { YLO, YHI } > 0.0 THE SOLUTION IS ALWAYS SUPERSATURATED WITH NH4CL ! YLO= Y1 ! Save Y-value at Hi position IF (YLO.GT.ZERO .AND. YHI.GT.ZERO) THEN Y3 = FUNCI5A (ZERO) GOTO 50 ELSE IF (ABS(Y2) .LT. EPS) THEN ! X2 IS A SOLUTION GOTO 50 ELSE CALL PUSHERR (0001, 'CALCI5') ! WARNING ERROR: NO SOLUTION GOTO 50 ENDIF ! ! *** PERFORM BISECTION *********************************************** ! 20 DO 30 I=1,MAXIT X3 = 0.5*(X1+X2) Y3 = FUNCI5A (X3) IF (SIGN(1.d0,Y1)*SIGN(1.d0,Y3) .LE. ZERO) THEN ! (Y1*Y3 .LE. ZERO) Y2 = Y3 X2 = X3 ELSE Y1 = Y3 X1 = X3 ENDIF IF (ABS(X2-X1) .LE. EPS*X1) GOTO 40 30 CONTINUE CALL PUSHERR (0002, 'CALCI5') ! WARNING ERROR: NO CONVERGENCE ! ! *** CONVERGED ; RETURN ********************************************** ! 40 X3 = 0.5*(X1+X2) Y3 = FUNCI5A (X3) ! 50 RETURN ! *** END OF SUBROUTINE CALCI5 ***************************************** ! END SUBROUTINE CALCI5 !======================================================================= ! ! *** ISORROPIA CODE ! *** SUBROUTINE FUNCI5A ! *** CASE I5 ! ! THE MAIN CHARACTERISTICS OF THIS REGIME ARE: ! 1. SULFATE RICH, NO FREE ACID (1.0 <= SULRAT < 2.0) ! 2. SOLID & LIQUID AEROSOL POSSIBLE ! 3. SOLIDS POSSIBLE : NA2SO4 ! ! *** COPYRIGHT 1996-2006, UNIVERSITY OF MIAMI, CARNEGIE MELLON UNIVERSITY, ! *** GEORGIA INSTITUTE OF TECHNOLOGY ! *** WRITTEN BY ATHANASIOS NENES ! *** UPDATED BY CHRISTOS FOUNTOUKIS ! !======================================================================= ! DOUBLE PRECISION FUNCTION FUNCI5A (P4) USE ISRPIA USE SOLUT IMPLICIT DOUBLE PRECISION (A-H, O-Z) ! INCLUDE 'ISRPIA.INC' ! COMMON /SOLUT/ CHI1, CHI2, CHI3, CHI4, CHI5, CHI6, CHI7, CHI8, & ! PSI1, PSI2, PSI3, PSI4, PSI5, PSI6, PSI7, PSI8, & ! A1, A2, A3, A4, A5, A6, A7, A8 ! ! *** SETUP PARAMETERS ************************************************ ! PSI4 = P4 ! PSI3 already assigned in FUNCI5A FRST = .TRUE. CALAIN = .TRUE. ! ! *** SOLVE EQUATIONS ; WITH ITERATIONS FOR ACTIVITY COEF. ************ ! DO 10 I=1,NSWEEP ! A4 = XK5 *(WATER/GAMA(2))**3.0 A5 = XK7 *(WATER/GAMA(4))**3.0 A6 = XK1 *WATER/GAMA(7)*(GAMA(8)/GAMA(7))**2. ! ! CALCULATE DISSOCIATION QUANTITIES ! BB = PSI2 + PSI4 + PSI5 + A6 ! PSI6 CC =-A6*(PSI2 + PSI3 + PSI1) DD = BB*BB - 4.D0*CC PSI6 = 0.5D0*(-BB + SQRT(DD)) ! ! *** CALCULATE SPECIATION ******************************************** ! MOLAL (1) = PSI6 ! HI MOLAL (2) = 2.D0*PSI4 + PSI3 ! NAI MOLAL (3) = 3.D0*PSI2 + 2.D0*PSI5 + PSI1 ! NH4I MOLAL (5) = PSI2 + PSI4 + PSI5 + PSI6 ! SO4I MOLAL (6) = PSI2 + PSI3 + PSI1 - PSI6 ! HSO4I CLC = ZERO CNAHSO4 = ZERO CNA2SO4 = CHI4 - PSI4 CNH42S4 = ZERO CNH4HS4 = ZERO CALL CALCMR ! Water content ! ! *** CALCULATE ACTIVITIES OR TERMINATE INTERNAL LOOP ***************** ! IF (FRST.AND.CALAOU .OR. .NOT.FRST.AND.CALAIN) THEN CALL CALCACT ELSE GOTO 20 ENDIF 10 CONTINUE ! ! *** CALCULATE OBJECTIVE FUNCTION ************************************ ! 20 A4 = XK5 *(WATER/GAMA(2))**3.0 FUNCI5A= MOLAL(5)*MOLAL(2)*MOLAL(2)/A4 - ONE RETURN ! ! *** END OF FUNCTION FUNCI5A ******************************************** ! END FUNCTION FUNCI5A !======================================================================= ! ! *** ISORROPIA CODE ! *** SUBROUTINE CALCI4 ! *** CASE I4 ! ! THE MAIN CHARACTERISTICS OF THIS REGIME ARE: ! 1. SULFATE RICH, NO FREE ACID (1.0 <= SULRAT < 2.0) ! 2. SOLID & LIQUID AEROSOL POSSIBLE ! 3. SOLIDS POSSIBLE : (NH4)2SO4, NA2SO4 ! ! *** COPYRIGHT 1996-2006, UNIVERSITY OF MIAMI, CARNEGIE MELLON UNIVERSITY, ! *** GEORGIA INSTITUTE OF TECHNOLOGY ! *** WRITTEN BY ATHANASIOS NENES ! *** UPDATED BY CHRISTOS FOUNTOUKIS ! !======================================================================= ! SUBROUTINE CALCI4 USE ISRPIA USE SOLUT IMPLICIT DOUBLE PRECISION (A-H, O-Z) ! INCLUDE 'ISRPIA.INC' ! COMMON /SOLUT/ CHI1, CHI2, CHI3, CHI4, CHI5, CHI6, CHI7, CHI8, & ! PSI1, PSI2, PSI3, PSI4, PSI5, PSI6, PSI7, PSI8, & ! A1, A2, A3, A4, A5, A6, A7, A8 ! ! *** FIND DRY COMPOSITION ********************************************** ! CALL CALCI1A ! ! *** SETUP PARAMETERS ************************************************ ! CHI1 = CNH4HS4 ! Save from CALCI1 run CHI2 = CLC CHI3 = CNAHSO4 CHI4 = CNA2SO4 CHI5 = CNH42S4 ! PSI1 = CNH4HS4 ! ASSIGN INITIAL PSI's PSI2 = CLC PSI3 = CNAHSO4 PSI4 = ZERO PSI5 = ZERO ! CALAOU = .TRUE. ! Outer loop activity calculation flag PSI4LO = ZERO ! Low limit PSI4HI = CHI4 ! High limit ! ! *** IF NA2SO4(S) =0, CALL FUNCI4B FOR Y4=0 *************************** ! IF (CHI4.LE.TINY) THEN Y1 = FUNCI4A (ZERO) GOTO 50 ENDIF ! ! *** INITIAL VALUES FOR BISECTION ************************************ ! X1 = PSI4HI Y1 = FUNCI4A (X1) YHI= Y1 ! Save Y-value at HI position ! ! *** YHI < 0.0 THE SOLUTION IS ALWAYS UNDERSATURATED WITH NA2SO4 ** ! IF (ABS(Y1).LE.EPS .OR. YHI.LT.ZERO) GOTO 50 ! ! *** ROOT TRACKING ; FOR THE RANGE OF HI AND LO ********************** ! DX = (PSI4HI-PSI4LO)/REAL(NDIV) DO 10 I=1,NDIV X2 = X1-DX Y2 = FUNCI4A (X2) IF (SIGN(1.d0,Y1)*SIGN(1.d0,Y2).LT.ZERO) GOTO 20 ! (Y1*Y2.LT.ZERO) X1 = X2 Y1 = Y2 10 CONTINUE ! ! *** { YLO, YHI } > 0.0 THE SOLUTION IS ALWAYS SUPERSATURATED WITH NH4CL ! YLO= Y1 ! Save Y-value at Hi position IF (YLO.GT.ZERO .AND. YHI.GT.ZERO) THEN Y3 = FUNCI4A (ZERO) GOTO 50 ELSE IF (ABS(Y2) .LT. EPS) THEN ! X2 IS A SOLUTION GOTO 50 ELSE CALL PUSHERR (0001, 'CALCI4') ! WARNING ERROR: NO SOLUTION GOTO 50 ENDIF ! ! *** PERFORM BISECTION *********************************************** ! 20 DO 30 I=1,MAXIT X3 = 0.5*(X1+X2) Y3 = FUNCI4A (X3) IF (SIGN(1.d0,Y1)*SIGN(1.d0,Y3) .LE. ZERO) THEN ! (Y1*Y3 .LE. ZERO) Y2 = Y3 X2 = X3 ELSE Y1 = Y3 X1 = X3 ENDIF IF (ABS(X2-X1) .LE. EPS*X1) GOTO 40 30 CONTINUE CALL PUSHERR (0002, 'CALCI4') ! WARNING ERROR: NO CONVERGENCE ! ! *** CONVERGED ; RETURN ********************************************** ! 40 X3 = 0.5*(X1+X2) Y3 = FUNCI4A (X3) ! 50 RETURN ! *** END OF SUBROUTINE CALCI4 ***************************************** ! END SUBROUTINE CALCI4 !======================================================================= ! ! *** ISORROPIA CODE ! *** SUBROUTINE FUNCI4A ! *** CASE I4 ! ! THE MAIN CHARACTERISTICS OF THIS REGIME ARE: ! 1. SULFATE RICH, NO FREE ACID (1.0 <= SULRAT < 2.0) ! 2. SOLID & LIQUID AEROSOL POSSIBLE ! 3. SOLIDS POSSIBLE : (NH4)2SO4, NA2SO4 ! ! *** COPYRIGHT 1996-2006, UNIVERSITY OF MIAMI, CARNEGIE MELLON UNIVERSITY, ! *** GEORGIA INSTITUTE OF TECHNOLOGY ! *** WRITTEN BY ATHANASIOS NENES ! *** UPDATED BY CHRISTOS FOUNTOUKIS ! !======================================================================= ! DOUBLE PRECISION FUNCTION FUNCI4A (P4) USE ISRPIA USE SOLUT IMPLICIT DOUBLE PRECISION (A-H, O-Z) ! INCLUDE 'ISRPIA.INC' ! COMMON /SOLUT/ CHI1, CHI2, CHI3, CHI4, CHI5, CHI6, CHI7, CHI8, & ! PSI1, PSI2, PSI3, PSI4, PSI5, PSI6, PSI7, PSI8, & ! A1, A2, A3, A4, A5, A6, A7, A8 ! ! *** SETUP PARAMETERS ************************************************ ! PSI4 = P4 ! PSI3 already assigned in FUNCI4A FRST = .TRUE. CALAIN = .TRUE. ! ! *** SOLVE EQUATIONS ; WITH ITERATIONS FOR ACTIVITY COEF. ************ ! DO 10 I=1,NSWEEP ! A4 = XK5 *(WATER/GAMA(2))**3.0 A5 = XK7 *(WATER/GAMA(4))**3.0 A6 = XK1 *WATER/GAMA(7)*(GAMA(8)/GAMA(7))**2. A7 = SQRT(A4/A5) ! ! CALCULATE DISSOCIATION QUANTITIES ! BB = PSI2 + PSI4 + PSI5 + A6 ! PSI6 CC =-A6*(PSI2 + PSI3 + PSI1) DD = BB*BB - 4.D0*CC PSI6 = 0.5D0*(-BB + SQRT(DD)) ! PSI5 = (PSI3 + 2.D0*PSI4 - A7*(3.D0*PSI2 + PSI1))/2.D0/A7 PSI5 = MIN (PSI5, CHI5) ! ! *** CALCULATE SPECIATION ******************************************** ! MOLAL (1) = PSI6 ! HI MOLAL (2) = 2.D0*PSI4 + PSI3 ! NAI MOLAL (3) = 3.D0*PSI2 + 2.D0*PSI5 + PSI1 ! NH4I MOLAL (5) = PSI2 + PSI4 + PSI5 + PSI6 ! SO4I MOLAL (6) = PSI2 + PSI3 + PSI1 - PSI6 ! HSO4I CLC = ZERO CNAHSO4 = ZERO CNA2SO4 = CHI4 - PSI4 CNH42S4 = CHI5 - PSI5 CNH4HS4 = ZERO CALL CALCMR ! Water content ! ! *** CALCULATE ACTIVITIES OR TERMINATE INTERNAL LOOP ***************** ! IF (FRST.AND.CALAOU .OR. .NOT.FRST.AND.CALAIN) THEN CALL CALCACT ELSE GOTO 20 ENDIF 10 CONTINUE ! ! *** CALCULATE OBJECTIVE FUNCTION ************************************ ! 20 A4 = XK5 *(WATER/GAMA(2))**3.0 FUNCI4A= MOLAL(5)*MOLAL(2)*MOLAL(2)/A4 - ONE RETURN ! ! *** END OF FUNCTION FUNCI4A ******************************************** ! END FUNCTION FUNCI4A !======================================================================= ! ! *** ISORROPIA CODE ! *** SUBROUTINE CALCI3 ! *** CASE I3 ! ! THE MAIN CHARACTERISTICS OF THIS REGIME ARE: ! 1. SULFATE RICH, NO FREE ACID (1.0 <= SULRAT < 2.0) ! 2. SOLID & LIQUID AEROSOL POSSIBLE ! 3. SOLIDS POSSIBLE : (NH4)2SO4, NA2SO4, NAHSO4, LC ! ! THERE ARE THREE REGIMES IN THIS CASE: ! 1.(NA,NH4)HSO4(s) POSSIBLE. LIQUID & SOLID AEROSOL (SUBROUTINE CALCI3A) ! 2.(NA,NH4)HSO4(s) NOT POSSIBLE, AND RH < MDRH. SOLID AEROSOL ONLY ! 3.(NA,NH4)HSO4(s) NOT POSSIBLE, AND RH >= MDRH. SOLID & LIQUID AEROSOL ! ! REGIMES 2. AND 3. ARE CONSIDERED TO BE THE SAME AS CASES I1A, I2B ! RESPECTIVELY ! ! *** COPYRIGHT 1996-2006, UNIVERSITY OF MIAMI, CARNEGIE MELLON UNIVERSITY, ! *** GEORGIA INSTITUTE OF TECHNOLOGY ! *** WRITTEN BY ATHANASIOS NENES ! *** UPDATED BY CHRISTOS FOUNTOUKIS ! !======================================================================= ! SUBROUTINE CALCI3 USE ISRPIA IMPLICIT DOUBLE PRECISION (A-H, O-Z) ! INCLUDE 'ISRPIA.INC' EXTERNAL CALCI1A, CALCI4 ! ! *** FIND DRY COMPOSITION ********************************************** ! CALL CALCI1A ! ! *** REGIME DEPENDS UPON THE POSSIBLE SOLIDS & RH ********************** ! IF (CNH4HS4.GT.TINY .OR. CNAHSO4.GT.TINY) THEN SCASE = 'I3 ; SUBCASE 1' CALL CALCI3A ! FULL SOLUTION SCASE = 'I3 ; SUBCASE 1' ENDIF ! IF (WATER.LE.TINY) THEN IF (RH.LT.DRMI3) THEN ! SOLID SOLUTION WATER = TINY DO 10 I=1,NIONS MOLAL(I) = ZERO 10 CONTINUE CALL CALCI1A SCASE = 'I3 ; SUBCASE 2' ! ELSEIF (RH.GE.DRMI3) THEN ! MDRH OF I3 SCASE = 'I3 ; SUBCASE 3' CALL CALCMDRH (RH, DRMI3, DRLC, CALCI1A, CALCI4) SCASE = 'I3 ; SUBCASE 3' ENDIF ENDIF ! RETURN ! ! *** END OF SUBROUTINE CALCI3 ****************************************** ! END SUBROUTINE CALCI3 !======================================================================= ! ! *** ISORROPIA CODE ! *** SUBROUTINE CALCI3A ! *** CASE I3 ; SUBCASE 1 ! ! THE MAIN CHARACTERISTICS OF THIS REGIME ARE: ! 1. SULFATE RICH, NO FREE ACID (1.0 <= SULRAT < 2.0) ! 2. SOLID & LIQUID AEROSOL POSSIBLE ! 3. SOLIDS POSSIBLE : (NH4)2SO4, NA2SO4, LC ! ! *** COPYRIGHT 1996-2006, UNIVERSITY OF MIAMI, CARNEGIE MELLON UNIVERSITY, ! *** GEORGIA INSTITUTE OF TECHNOLOGY ! *** WRITTEN BY ATHANASIOS NENES ! *** UPDATED BY CHRISTOS FOUNTOUKIS ! !======================================================================= ! SUBROUTINE CALCI3A USE ISRPIA USE SOLUT IMPLICIT DOUBLE PRECISION (A-H, O-Z) ! INCLUDE 'ISRPIA.INC' ! COMMON /SOLUT/ CHI1, CHI2, CHI3, CHI4, CHI5, CHI6, CHI7, CHI8, & ! PSI1, PSI2, PSI3, PSI4, PSI5, PSI6, PSI7, PSI8, & ! A1, A2, A3, A4, A5, A6, A7, A8 ! ! *** FIND DRY COMPOSITION ********************************************** ! CALL CALCI1A ! Needed when called from CALCMDRH ! ! *** SETUP PARAMETERS ************************************************ ! CHI1 = CNH4HS4 ! Save from CALCI1 run CHI2 = CLC CHI3 = CNAHSO4 CHI4 = CNA2SO4 CHI5 = CNH42S4 ! PSI1 = CNH4HS4 ! ASSIGN INITIAL PSI's PSI2 = ZERO PSI3 = CNAHSO4 PSI4 = ZERO PSI5 = ZERO ! CALAOU = .TRUE. ! Outer loop activity calculation flag PSI2LO = ZERO ! Low limit PSI2HI = CHI2 ! High limit ! ! *** INITIAL VALUES FOR BISECTION ************************************ ! X1 = PSI2HI Y1 = FUNCI3A (X1) YHI= Y1 ! Save Y-value at HI position ! ! *** YHI < 0.0 THE SOLUTION IS ALWAYS UNDERSATURATED WITH LC ********* ! IF (YHI.LT.EPS) GOTO 50 ! ! *** ROOT TRACKING ; FOR THE RANGE OF HI AND LO ********************** ! DX = (PSI2HI-PSI2LO)/REAL(NDIV) DO 10 I=1,NDIV X2 = MAX(X1-DX, PSI2LO) Y2 = FUNCI3A (X2) IF (SIGN(1.d0,Y1)*SIGN(1.d0,Y2).LT.ZERO) GOTO 20 ! (Y1*Y2.LT.ZERO) X1 = X2 Y1 = Y2 10 CONTINUE ! ! *** { YLO, YHI } > 0.0 THE SOLUTION IS ALWAYS SUPERSATURATED WITH LC ! IF (Y2.GT.EPS) Y2 = FUNCI3A (ZERO) GOTO 50 ! ! *** PERFORM BISECTION *********************************************** ! 20 DO 30 I=1,MAXIT X3 = 0.5*(X1+X2) Y3 = FUNCI3A (X3) IF (SIGN(1.d0,Y1)*SIGN(1.d0,Y3) .LE. ZERO) THEN ! (Y1*Y3 .LE. ZERO) Y2 = Y3 X2 = X3 ELSE Y1 = Y3 X1 = X3 ENDIF IF (ABS(X2-X1) .LE. EPS*X1) GOTO 40 30 CONTINUE CALL PUSHERR (0002, 'CALCI3A') ! WARNING ERROR: NO CONVERGENCE ! ! *** CONVERGED ; RETURN ********************************************** ! 40 X3 = 0.5*(X1+X2) Y3 = FUNCI3A (X3) ! 50 RETURN ! *** END OF SUBROUTINE CALCI3A ***************************************** ! END SUBROUTINE CALCI3A ! !======================================================================= ! ! *** ISORROPIA CODE ! *** SUBROUTINE FUNCI3A ! *** CASE I3 ; SUBCASE 1 ! ! THE MAIN CHARACTERISTICS OF THIS REGIME ARE: ! 1. SULFATE RICH, NO FREE ACID (1.0 <= SULRAT < 2.0) ! 2. SOLID & LIQUID AEROSOL POSSIBLE ! 3. SOLIDS POSSIBLE : (NH4)2SO4, NA2SO4, LC ! ! *** COPYRIGHT 1996-2006, UNIVERSITY OF MIAMI, CARNEGIE MELLON UNIVERSITY, ! *** GEORGIA INSTITUTE OF TECHNOLOGY ! *** WRITTEN BY ATHANASIOS NENES ! *** UPDATED BY CHRISTOS FOUNTOUKIS ! !======================================================================= ! DOUBLE PRECISION FUNCTION FUNCI3A (P2) USE ISRPIA USE SOLUT IMPLICIT DOUBLE PRECISION (A-H, O-Z) ! INCLUDE 'ISRPIA.INC' ! COMMON /SOLUT/ CHI1, CHI2, CHI3, CHI4, CHI5, CHI6, CHI7, CHI8, & ! PSI1, PSI2, PSI3, PSI4, PSI5, PSI6, PSI7, PSI8, & ! A1, A2, A3, A4, A5, A6, A7, A8 ! ! *** SETUP PARAMETERS ************************************************ ! PSI2 = P2 ! Save PSI2 in COMMON BLOCK PSI4LO = ZERO ! Low limit for PSI4 PSI4HI = CHI4 ! High limit for PSI4 ! ! *** IF NH3 =0, CALL FUNCI3B FOR Y4=0 ******************************** ! IF (CHI4.LE.TINY) THEN FUNCI3A = FUNCI3B (ZERO) GOTO 50 ENDIF ! ! *** INITIAL VALUES FOR BISECTION ************************************ ! X1 = PSI4HI Y1 = FUNCI3B (X1) IF (ABS(Y1).LE.EPS) GOTO 50 YHI= Y1 ! Save Y-value at HI position ! ! *** YHI < 0.0 THE SOLUTION IS ALWAYS UNDERSATURATED WITH NA2SO4 ***** ! IF (YHI.LT.ZERO) GOTO 50 ! ! *** ROOT TRACKING ; FOR THE RANGE OF HI AND LO ********************** ! DX = (PSI4HI-PSI4LO)/REAL(NDIV) DO 10 I=1,NDIV X2 = MAX(X1-DX, PSI4LO) Y2 = FUNCI3B (X2) IF (SIGN(1.d0,Y1)*SIGN(1.d0,Y2).LT.ZERO) GOTO 20 ! (Y1*Y2.LT.ZERO) X1 = X2 Y1 = Y2 10 CONTINUE ! ! *** { YLO, YHI } > 0.0 THE SOLUTION IS ALWAYS SUPERSATURATED WITH NA2SO4 ! IF (Y2.GT.EPS) Y2 = FUNCI3B (PSI4LO) GOTO 50 ! ! *** PERFORM BISECTION *********************************************** ! 20 DO 30 I=1,MAXIT X3 = 0.5*(X1+X2) Y3 = FUNCI3B (X3) IF (SIGN(1.d0,Y1)*SIGN(1.d0,Y3) .LE. ZERO) THEN ! (Y1*Y3 .LE. ZERO) Y2 = Y3 X2 = X3 ELSE Y1 = Y3 X1 = X3 ENDIF IF (ABS(X2-X1) .LE. EPS*X1) GOTO 40 30 CONTINUE CALL PUSHERR (0004, 'FUNCI3A') ! WARNING ERROR: NO CONVERGENCE ! ! *** INNER LOOP CONVERGED ********************************************** ! 40 X3 = 0.5*(X1+X2) Y3 = FUNCI3B (X3) ! ! *** CALCULATE FUNCTION VALUE FOR OUTER LOOP *************************** ! 50 A2 = XK13*(WATER/GAMA(13))**5.0 FUNCI3A = MOLAL(5)*MOLAL(6)*MOLAL(3)**3.D0/A2 - ONE RETURN ! ! *** END OF FUNCTION FUNCI3A ******************************************* ! END FUNCTION FUNCI3A !======================================================================= ! ! *** ISORROPIA CODE ! *** FUNCTION FUNCI3B ! *** CASE I3 ; SUBCASE 2 ! ! THE MAIN CHARACTERISTICS OF THIS REGIME ARE: ! 1. SULFATE RICH, NO FREE ACID (1.0 <= SULRAT < 2.0) ! 2. SOLID & LIQUID AEROSOL POSSIBLE ! 3. SOLIDS POSSIBLE : (NH4)2SO4, NA2SO4, LC ! ! SOLUTION IS SAVED IN COMMON BLOCK /CASE/ ! !======================================================================= ! DOUBLE PRECISION FUNCTION FUNCI3B (P4) USE ISRPIA USE SOLUT IMPLICIT DOUBLE PRECISION (A-H, O-Z) ! INCLUDE 'ISRPIA.INC' ! ! COMMON /SOLUT/ CHI1, CHI2, CHI3, CHI4, CHI5, CHI6, CHI7, CHI8, & ! PSI1, PSI2, PSI3, PSI4, PSI5, PSI6, PSI7, PSI8, & ! A1, A2, A3, A4, A5, A6, A7, A8 ! ! *** SETUP PARAMETERS ************************************************ ! PSI4 = P4 ! ! *** SETUP PARAMETERS ************************************************ ! FRST = .TRUE. CALAIN = .TRUE. ! ! *** SOLVE EQUATIONS ; WITH ITERATIONS FOR ACTIVITY COEF. ************ ! DO 10 I=1,NSWEEP ! A4 = XK5*(WATER/GAMA(2))**3.0 A5 = XK7*(WATER/GAMA(4))**3.0 A6 = XK1*WATER/GAMA(7)*(GAMA(8)/GAMA(7))**2. A7 = SQRT(A4/A5) ! ! CALCULATE DISSOCIATION QUANTITIES ! BB = PSI2 + PSI4 + PSI5 + A6 ! PSI6 CC =-A6*(PSI2 + PSI3 + PSI1) DD = BB*BB - 4.D0*CC PSI6 = 0.5D0*(-BB + SQRT(DD)) ! PSI5 = (PSI3 + 2.D0*PSI4 - A7*(3.D0*PSI2 + PSI1))/2.D0/A7 PSI5 = MIN (PSI5, CHI5) ! ! *** CALCULATE SPECIATION ******************************************** ! MOLAL(1) = PSI6 ! HI MOLAL(2) = 2.D0*PSI4 + PSI3 ! NAI MOLAL(3) = 3.D0*PSI2 + 2.D0*PSI5 + PSI1 ! NH4I MOLAL(5) = PSI2 + PSI4 + PSI5 + PSI6 ! SO4I MOLAL(6) = MAX(PSI2 + PSI3 + PSI1 - PSI6, TINY) ! HSO4I CLC = MAX(CHI2 - PSI2, ZERO) CNAHSO4 = ZERO CNA2SO4 = MAX(CHI4 - PSI4, ZERO) CNH42S4 = MAX(CHI5 - PSI5, ZERO) CNH4HS4 = ZERO CALL CALCMR ! Water content ! ! *** CALCULATE ACTIVITIES OR TERMINATE INTERNAL LOOP ***************** ! IF (FRST.AND.CALAOU .OR. .NOT.FRST.AND.CALAIN) THEN CALL CALCACT ELSE GOTO 20 ENDIF 10 CONTINUE ! ! *** CALCULATE OBJECTIVE FUNCTION ************************************ ! 20 A4 = XK5 *(WATER/GAMA(2))**3.0 FUNCI3B= MOLAL(5)*MOLAL(2)*MOLAL(2)/A4 - ONE RETURN ! ! *** END OF FUNCTION FUNCI3B ******************************************** ! END FUNCTION FUNCI3B !======================================================================= ! ! *** ISORROPIA CODE ! *** SUBROUTINE CALCI2 ! *** CASE I2 ! ! THE MAIN CHARACTERISTICS OF THIS REGIME ARE: ! 1. SULFATE RICH, NO FREE ACID (1.0 <= SULRAT < 2.0) ! 2. SOLID & LIQUID AEROSOL POSSIBLE ! 3. SOLIDS POSSIBLE : (NH4)2SO4, NA2SO4, NAHSO4, LC ! ! THERE ARE THREE REGIMES IN THIS CASE: ! 1. NH4HSO4(s) POSSIBLE. LIQUID & SOLID AEROSOL (SUBROUTINE CALCI2A) ! 2. NH4HSO4(s) NOT POSSIBLE, AND RH < MDRH. SOLID AEROSOL ONLY ! 3. NH4HSO4(s) NOT POSSIBLE, AND RH >= MDRH. SOLID & LIQUID AEROSOL ! ! REGIMES 2. AND 3. ARE CONSIDERED TO BE THE SAME AS CASES I1A, I2B ! RESPECTIVELY ! ! *** COPYRIGHT 1996-2006, UNIVERSITY OF MIAMI, CARNEGIE MELLON UNIVERSITY, ! *** GEORGIA INSTITUTE OF TECHNOLOGY ! *** WRITTEN BY ATHANASIOS NENES ! *** UPDATED BY CHRISTOS FOUNTOUKIS ! !======================================================================= ! SUBROUTINE CALCI2 USE ISRPIA IMPLICIT DOUBLE PRECISION (A-H, O-Z) ! INCLUDE 'ISRPIA.INC' EXTERNAL CALCI1A, CALCI3A ! ! *** FIND DRY COMPOSITION ********************************************** ! CALL CALCI1A ! ! *** REGIME DEPENDS UPON THE POSSIBLE SOLIDS & RH ********************** ! IF (CNH4HS4.GT.TINY) THEN SCASE = 'I2 ; SUBCASE 1' CALL CALCI2A SCASE = 'I2 ; SUBCASE 1' ENDIF ! IF (WATER.LE.TINY) THEN IF (RH.LT.DRMI2) THEN ! SOLID SOLUTION ONLY WATER = TINY DO 10 I=1,NIONS MOLAL(I) = ZERO 10 CONTINUE CALL CALCI1A SCASE = 'I2 ; SUBCASE 2' ! ELSEIF (RH.GE.DRMI2) THEN ! MDRH OF I2 SCASE = 'I2 ; SUBCASE 3' CALL CALCMDRH (RH, DRMI2, DRNAHSO4, CALCI1A, CALCI3A) SCASE = 'I2 ; SUBCASE 3' ENDIF ENDIF ! RETURN ! ! *** END OF SUBROUTINE CALCI2 ****************************************** ! END SUBROUTINE CALCI2 !======================================================================= ! ! *** ISORROPIA CODE ! *** SUBROUTINE CALCI2A ! *** CASE I2 ; SUBCASE A ! ! THE MAIN CHARACTERISTICS OF THIS REGIME ARE: ! 1. SULFATE RICH, NO FREE ACID (1.0 <= SULRAT < 2.0) ! 2. SOLID & LIQUID AEROSOL POSSIBLE ! 3. SOLIDS POSSIBLE : (NH4)2SO4, NA2SO4, NAHSO4, LC ! ! *** COPYRIGHT 1996-2006, UNIVERSITY OF MIAMI, CARNEGIE MELLON UNIVERSITY, ! *** GEORGIA INSTITUTE OF TECHNOLOGY ! *** WRITTEN BY ATHANASIOS NENES ! *** UPDATED BY CHRISTOS FOUNTOUKIS ! !======================================================================= ! SUBROUTINE CALCI2A USE ISRPIA USE SOLUT IMPLICIT DOUBLE PRECISION (A-H, O-Z) ! INCLUDE 'ISRPIA.INC' ! COMMON /SOLUT/ CHI1, CHI2, CHI3, CHI4, CHI5, CHI6, CHI7, CHI8, & ! PSI1, PSI2, PSI3, PSI4, PSI5, PSI6, PSI7, PSI8, & ! A1, A2, A3, A4, A5, A6, A7, A8 ! ! *** FIND DRY COMPOSITION ********************************************** ! CALL CALCI1A ! Needed when called from CALCMDRH ! ! *** SETUP PARAMETERS ************************************************ ! CHI1 = CNH4HS4 ! Save from CALCI1 run CHI2 = CLC CHI3 = CNAHSO4 CHI4 = CNA2SO4 CHI5 = CNH42S4 ! PSI1 = CNH4HS4 ! ASSIGN INITIAL PSI's PSI2 = ZERO PSI3 = ZERO PSI4 = ZERO PSI5 = ZERO ! CALAOU = .TRUE. ! Outer loop activity calculation flag PSI2LO = ZERO ! Low limit PSI2HI = CHI2 ! High limit ! ! *** INITIAL VALUES FOR BISECTION ************************************ ! X1 = PSI2HI Y1 = FUNCI2A (X1) YHI= Y1 ! Save Y-value at HI position ! ! *** YHI < 0.0 THE SOLUTION IS ALWAYS UNDERSATURATED WITH LC ********* ! IF (YHI.LT.EPS) GOTO 50 ! ! *** ROOT TRACKING ; FOR THE RANGE OF HI AND LO ********************** ! DX = (PSI2HI-PSI2LO)/REAL(NDIV) DO 10 I=1,NDIV X2 = MAX(X1-DX, PSI2LO) Y2 = FUNCI2A (X2) IF (SIGN(1.d0,Y1)*SIGN(1.d0,Y2).LT.ZERO) GOTO 20 ! (Y1*Y2.LT.ZERO) X1 = X2 Y1 = Y2 10 CONTINUE ! ! *** { YLO, YHI } > 0.0 THE SOLUTION IS ALWAYS SUPERSATURATED WITH LC ! IF (Y2.GT.EPS) Y2 = FUNCI3A (ZERO) GOTO 50 ! ! *** PERFORM BISECTION *********************************************** ! 20 DO 30 I=1,MAXIT X3 = 0.5*(X1+X2) Y3 = FUNCI2A (X3) IF (SIGN(1.d0,Y1)*SIGN(1.d0,Y3) .LE. ZERO) THEN ! (Y1*Y3 .LE. ZERO) Y2 = Y3 X2 = X3 ELSE Y1 = Y3 X1 = X3 ENDIF IF (ABS(X2-X1) .LE. EPS*X1) GOTO 40 30 CONTINUE CALL PUSHERR (0002, 'CALCI2A') ! WARNING ERROR: NO CONVERGENCE ! ! *** CONVERGED ; RETURN ********************************************** ! 40 X3 = 0.5*(X1+X2) Y3 = FUNCI2A (X3) ! 50 RETURN ! *** END OF SUBROUTINE CALCI2A ***************************************** ! END SUBROUTINE CALCI2A !======================================================================= ! ! *** ISORROPIA CODE ! *** SUBROUTINE FUNCI2A ! *** CASE I2 ; SUBCASE 1 ! ! THE MAIN CHARACTERISTICS OF THIS REGIME ARE: ! 1. SULFATE RICH, NO FREE ACID (1.0 <= SULRAT < 2.0) ! 2. SOLID & LIQUID AEROSOL POSSIBLE ! 3. SOLIDS POSSIBLE : (NH4)2SO4, NA2SO4, NAHSO4, LC ! ! *** COPYRIGHT 1996-2006, UNIVERSITY OF MIAMI, CARNEGIE MELLON UNIVERSITY, ! *** GEORGIA INSTITUTE OF TECHNOLOGY ! *** WRITTEN BY ATHANASIOS NENES ! *** UPDATED BY CHRISTOS FOUNTOUKIS ! !======================================================================= ! DOUBLE PRECISION FUNCTION FUNCI2A (P2) USE ISRPIA USE SOLUT IMPLICIT DOUBLE PRECISION (A-H, O-Z) ! INCLUDE 'ISRPIA.INC' ! COMMON /SOLUT/ CHI1, CHI2, CHI3, CHI4, CHI5, CHI6, CHI7, CHI8, & ! PSI1, PSI2, PSI3, PSI4, PSI5, PSI6, PSI7, PSI8, & ! A1, A2, A3, A4, A5, A6, A7, A8 ! ! *** SETUP PARAMETERS ************************************************ ! FRST = .TRUE. CALAIN = .TRUE. PSI2 = P2 ! Save PSI2 in COMMON BLOCK PSI3 = CHI3 PSI4 = CHI4 PSI5 = CHI5 PSI6 = ZERO ! ! *** SOLVE EQUATIONS ; WITH ITERATIONS FOR ACTIVITY COEF. ************ ! DO 10 I=1,NSWEEP ! A3 = XK11*(WATER/GAMA(12))**2.0 A4 = XK5 *(WATER/GAMA(2))**3.0 A5 = XK7 *(WATER/GAMA(4))**3.0 A6 = XK1 *WATER/GAMA(7)*(GAMA(8)/GAMA(7))**2. A7 = SQRT(A4/A5) ! ! CALCULATE DISSOCIATION QUANTITIES ! IF (CHI5.GT.TINY .AND. WATER.GT.TINY) THEN PSI5 = (PSI3 + 2.D0*PSI4 - A7*(3.D0*PSI2 + PSI1))/2.D0/A7 PSI5 = MAX(MIN (PSI5, CHI5), TINY) ENDIF ! IF (CHI4.GT.TINY .AND. WATER.GT.TINY) THEN AA = PSI2+PSI5+PSI6+PSI3 BB = PSI3*AA CC = 0.25D0*(PSI3*PSI3*(PSI2+PSI5+PSI6)-A4) CALL POLY3 (AA, BB, CC, PSI4, ISLV) IF (ISLV.EQ.0) THEN PSI4 = MIN (PSI4, CHI4) ELSE PSI4 = ZERO ENDIF ENDIF ! IF (CHI3.GT.TINY .AND. WATER.GT.TINY) THEN AA = 2.D0*PSI4 + PSI2 + PSI1 - PSI6 BB = 2.D0*PSI4*(PSI2 + PSI1 - PSI6) - A3 CC = ZERO CALL POLY3 (AA, BB, CC, PSI3, ISLV) IF (ISLV.EQ.0) THEN PSI3 = MIN (PSI3, CHI3) ELSE PSI3 = ZERO ENDIF ENDIF ! BB = PSI2 + PSI4 + PSI5 + A6 ! PSI6 CC =-A6*(PSI2 + PSI3 + PSI1) DD = BB*BB - 4.D0*CC PSI6 = 0.5D0*(-BB + SQRT(DD)) ! ! *** CALCULATE SPECIATION ******************************************** ! MOLAL (1) = PSI6 ! HI MOLAL (2) = 2.D0*PSI4 + PSI3 ! NAI MOLAL (3) = 3.D0*PSI2 + 2.D0*PSI5 + PSI1 ! NH4I MOLAL (5) = PSI2 + PSI4 + PSI5 + PSI6 ! SO4I MOLAL (6) = PSI2 + PSI3 + PSI1 - PSI6 ! HSO4I CLC = CHI2 - PSI2 CNAHSO4 = CHI3 - PSI3 CNA2SO4 = CHI4 - PSI4 CNH42S4 = CHI5 - PSI5 CNH4HS4 = ZERO CALL CALCMR ! Water content ! ! *** CALCULATE ACTIVITIES OR TERMINATE INTERNAL LOOP ***************** ! IF (FRST.AND.CALAOU .OR. .NOT.FRST.AND.CALAIN) THEN CALL CALCACT ELSE GOTO 20 ENDIF 10 CONTINUE ! ! *** CALCULATE FUNCTION VALUE FOR OUTER LOOP *************************** ! 20 A2 = XK13*(WATER/GAMA(13))**5.0 FUNCI2A = MOLAL(5)*MOLAL(6)*MOLAL(3)**3.D0/A2 - ONE RETURN ! ! *** END OF FUNCTION FUNCI2A ******************************************* ! END FUNCTION FUNCI2A !======================================================================= ! ! *** ISORROPIA CODE ! *** SUBROUTINE CALCI1 ! *** CASE I1 ! ! THE MAIN CHARACTERISTICS OF THIS REGIME ARE: ! 1. SULFATE RICH, NO FREE ACID (1.0 <= SULRAT < 2.0) ! 2. SOLID AEROSOL ONLY ! 3. SOLIDS POSSIBLE : NH4NO3, NH4CL, NA2SO4 ! ! THERE ARE TWO POSSIBLE REGIMES HERE, DEPENDING ON RELATIVE HUMIDITY: ! 1. WHEN RH >= MDRH ; LIQUID PHASE POSSIBLE (MDRH REGION) ! 2. WHEN RH < MDRH ; ONLY SOLID PHASE POSSIBLE (SUBROUTINE CALCI1A) ! ! *** COPYRIGHT 1996-2006, UNIVERSITY OF MIAMI, CARNEGIE MELLON UNIVERSITY, ! *** GEORGIA INSTITUTE OF TECHNOLOGY ! *** WRITTEN BY ATHANASIOS NENES ! *** UPDATED BY CHRISTOS FOUNTOUKIS ! !======================================================================= ! SUBROUTINE CALCI1 USE ISRPIA IMPLICIT DOUBLE PRECISION (A-H, O-Z) ! INCLUDE 'ISRPIA.INC' EXTERNAL CALCI1A, CALCI2A ! ! *** REGIME DEPENDS UPON THE AMBIENT RELATIVE HUMIDITY ***************** ! IF (RH.LT.DRMI1) THEN SCASE = 'I1 ; SUBCASE 1' CALL CALCI1A ! SOLID PHASE ONLY POSSIBLE SCASE = 'I1 ; SUBCASE 1' ELSE SCASE = 'I1 ; SUBCASE 2' ! LIQUID & SOLID PHASE POSSIBLE CALL CALCMDRH (RH, DRMI1, DRNH4HS4, CALCI1A, CALCI2A) SCASE = 'I1 ; SUBCASE 2' ENDIF ! ! *** AMMONIA IN GAS PHASE ********************************************** ! ! CALL CALCNH3 ! RETURN ! ! *** END OF SUBROUTINE CALCI1 ****************************************** ! END SUBROUTINE CALCI1 !======================================================================= ! ! *** ISORROPIA CODE ! *** SUBROUTINE CALCI1A ! *** CASE I1 ; SUBCASE 1 ! ! THE MAIN CHARACTERISTICS OF THIS REGIME ARE: ! 1. SULFATE RICH, NO FREE ACID (1.0 <= SULRAT < 2.0) ! 2. SOLID AEROSOL ONLY ! 3. SOLIDS POSSIBLE : NH4HSO4, NAHSO4, (NH4)2SO4, NA2SO4, LC ! ! *** COPYRIGHT 1996-2006, UNIVERSITY OF MIAMI, CARNEGIE MELLON UNIVERSITY, ! *** GEORGIA INSTITUTE OF TECHNOLOGY ! *** WRITTEN BY ATHANASIOS NENES ! *** UPDATED BY CHRISTOS FOUNTOUKIS ! !======================================================================= ! SUBROUTINE CALCI1A USE ISRPIA IMPLICIT DOUBLE PRECISION (A-H, O-Z) ! INCLUDE 'ISRPIA.INC' ! ! *** CALCULATE NON VOLATILE SOLIDS *********************************** ! CNA2SO4 = 0.5D0*W(1) CNH4HS4 = ZERO CNAHSO4 = ZERO CNH42S4 = ZERO FRSO4 = MAX(W(2)-CNA2SO4, ZERO) ! CLC = MIN(W(3)/3.D0, FRSO4/2.D0) FRSO4 = MAX(FRSO4-2.D0*CLC, ZERO) FRNH4 = MAX(W(3)-3.D0*CLC, ZERO) ! IF (FRSO4.LE.TINY) THEN CLC = MAX(CLC - FRNH4, ZERO) CNH42S4 = 2.D0*FRNH4 ELSEIF (FRNH4.LE.TINY) THEN CNH4HS4 = 3.D0*MIN(FRSO4, CLC) CLC = MAX(CLC-FRSO4, ZERO) IF (CNA2SO4.GT.TINY) THEN FRSO4 = MAX(FRSO4-CNH4HS4/3.D0, ZERO) CNAHSO4 = 2.D0*FRSO4 CNA2SO4 = MAX(CNA2SO4-FRSO4, ZERO) ENDIF ENDIF ! ! *** CALCULATE GAS SPECIES ********************************************* ! GHNO3 = W(4) GHCL = W(5) GNH3 = ZERO ! RETURN ! ! *** END OF SUBROUTINE CALCI1A ***************************************** ! END SUBROUTINE CALCI1A !======================================================================= ! ! *** ISORROPIA CODE ! *** SUBROUTINE CALCJ3 ! *** CASE J3 ! ! THE MAIN CHARACTERISTICS OF THIS REGIME ARE: ! 1. SULFATE RICH, FREE ACID (SULRAT < 1.0) ! 2. THERE IS ONLY A LIQUID PHASE ! ! *** COPYRIGHT 1996-2006, UNIVERSITY OF MIAMI, CARNEGIE MELLON UNIVERSITY, ! *** GEORGIA INSTITUTE OF TECHNOLOGY ! *** WRITTEN BY ATHANASIOS NENES ! *** UPDATED BY CHRISTOS FOUNTOUKIS ! !======================================================================= ! SUBROUTINE CALCJ3 USE ISRPIA IMPLICIT DOUBLE PRECISION (A-H, O-Z) ! INCLUDE 'ISRPIA.INC' ! DOUBLE PRECISION LAMDA, KAPA ! ! *** SETUP PARAMETERS ************************************************ ! CALAOU = .TRUE. ! Outer loop activity calculation flag FRST = .TRUE. CALAIN = .TRUE. ! LAMDA = MAX(W(2) - W(3) - W(1), TINY) ! FREE H2SO4 CHI1 = W(1) ! NA TOTAL as NaHSO4 CHI2 = W(3) ! NH4 TOTAL as NH4HSO4 PSI1 = CHI1 PSI2 = CHI2 ! ALL NH4HSO4 DELIQUESCED ! ! *** SOLVE EQUATIONS ; WITH ITERATIONS FOR ACTIVITY COEF. ************ ! DO 10 I=1,NSWEEP ! A3 = XK1 *WATER/GAMA(7)*(GAMA(8)/GAMA(7))**2.0 ! ! CALCULATE DISSOCIATION QUANTITIES ! BB = A3+LAMDA ! KAPA CC =-A3*(LAMDA + PSI1 + PSI2) DD = BB*BB-4.D0*CC KAPA = 0.5D0*(-BB+SQRT(DD)) ! ! *** CALCULATE SPECIATION ******************************************** ! MOLAL (1) = LAMDA + KAPA ! HI MOLAL (2) = PSI1 ! NAI MOLAL (3) = PSI2 ! NH4I MOLAL (4) = ZERO ! CLI MOLAL (5) = KAPA ! SO4I MOLAL (6) = LAMDA + PSI1 + PSI2 - KAPA ! HSO4I MOLAL (7) = ZERO ! NO3I ! CNAHSO4 = ZERO CNH4HS4 = ZERO ! CALL CALCMR ! Water content ! ! *** CALCULATE ACTIVITIES OR TERMINATE INTERNAL LOOP ***************** ! IF (FRST.AND.CALAOU .OR. .NOT.FRST.AND.CALAIN) THEN CALL CALCACT ELSE GOTO 50 ENDIF 10 CONTINUE ! 50 RETURN ! ! *** END OF SUBROUTINE CALCJ3 ****************************************** ! END SUBROUTINE CALCJ3 !======================================================================= ! ! *** ISORROPIA CODE ! *** SUBROUTINE CALCJ2 ! *** CASE J2 ! ! THE MAIN CHARACTERISTICS OF THIS REGIME ARE: ! 1. SULFATE RICH, FREE ACID (SULRAT < 1.0) ! 2. THERE IS BOTH A LIQUID & SOLID PHASE ! 3. SOLIDS POSSIBLE : NAHSO4 ! ! *** COPYRIGHT 1996-2006, UNIVERSITY OF MIAMI, CARNEGIE MELLON UNIVERSITY, ! *** GEORGIA INSTITUTE OF TECHNOLOGY ! *** WRITTEN BY ATHANASIOS NENES ! *** UPDATED BY CHRISTOS FOUNTOUKIS ! !======================================================================= ! SUBROUTINE CALCJ2 USE ISRPIA USE CASEJ IMPLICIT DOUBLE PRECISION (A-H, O-Z) ! INCLUDE 'ISRPIA.INC' ! ! DOUBLE PRECISION LAMDA, KAPA ! COMMON /CASEJ/ CHI1, CHI2, CHI3, LAMDA, KAPA, PSI1, PSI2, PSI3, & ! A1, A2, A3 ! ! *** SETUP PARAMETERS ************************************************ ! CALAOU = .TRUE. ! Outer loop activity calculation flag CHI1 = W(1) ! NA TOTAL CHI2 = W(3) ! NH4 TOTAL PSI1LO = TINY ! Low limit PSI1HI = CHI1 ! High limit ! ! *** INITIAL VALUES FOR BISECTION ************************************ ! X1 = PSI1HI Y1 = FUNCJ2 (X1) YHI= Y1 ! Save Y-value at HI position ! ! *** YHI < 0.0 THE SOLUTION IS ALWAYS UNDERSATURATED WITH NH42SO4 **** ! IF (ABS(Y1).LE.EPS .OR. YHI.LT.ZERO) GOTO 50 ! ! *** ROOT TRACKING ; FOR THE RANGE OF HI AND LO ********************** ! DX = (PSI1HI-PSI1LO)/REAL(NDIV) DO 10 I=1,NDIV X2 = X1-DX Y2 = FUNCJ2 (X2) IF (SIGN(1.d0,Y1)*SIGN(1.d0,Y2).LT.ZERO) GOTO 20 ! (Y1*Y2.LT.ZERO) X1 = X2 Y1 = Y2 10 CONTINUE ! ! *** { YLO, YHI } > 0.0 THE SOLUTION IS ALWAYS SUPERSATURATED WITH NH42SO4 ! YLO= Y1 ! Save Y-value at Hi position IF (YLO.GT.ZERO .AND. YHI.GT.ZERO) THEN Y3 = FUNCJ2 (ZERO) GOTO 50 ELSE IF (ABS(Y2) .LT. EPS) THEN ! X2 IS A SOLUTION GOTO 50 ELSE CALL PUSHERR (0001, 'CALCJ2') ! WARNING ERROR: NO SOLUTION GOTO 50 ENDIF ! ! *** PERFORM BISECTION *********************************************** ! 20 DO 30 I=1,MAXIT X3 = 0.5*(X1+X2) Y3 = FUNCJ2 (X3) IF (SIGN(1.d0,Y1)*SIGN(1.d0,Y3) .LE. ZERO) THEN ! (Y1*Y3 .LE. ZERO) Y2 = Y3 X2 = X3 ELSE Y1 = Y3 X1 = X3 ENDIF IF (ABS(X2-X1) .LE. EPS*X1) GOTO 40 30 CONTINUE CALL PUSHERR (0002, 'CALCJ2') ! WARNING ERROR: NO CONVERGENCE ! ! *** CONVERGED ; RETURN ********************************************** ! 40 X3 = 0.5*(X1+X2) Y3 = FUNCJ2 (X3) ! 50 RETURN ! ! *** END OF SUBROUTINE CALCJ2 ****************************************** ! END SUBROUTINE CALCJ2 !======================================================================= ! ! *** ISORROPIA CODE ! *** SUBROUTINE FUNCJ2 ! *** CASE J2 ! ! THE MAIN CHARACTERISTICS OF THIS REGIME ARE: ! 1. SULFATE RICH, FREE ACID (SULRAT < 1.0) ! 2. THERE IS BOTH A LIQUID & SOLID PHASE ! 3. SOLIDS POSSIBLE : (NH4)2SO4, NH4CL, NA2SO4 ! ! *** COPYRIGHT 1996-2006, UNIVERSITY OF MIAMI, CARNEGIE MELLON UNIVERSITY, ! *** GEORGIA INSTITUTE OF TECHNOLOGY ! *** WRITTEN BY ATHANASIOS NENES ! *** UPDATED BY CHRISTOS FOUNTOUKIS ! !======================================================================= ! DOUBLE PRECISION FUNCTION FUNCJ2 (P1) USE CASEJ USE ISRPIA IMPLICIT DOUBLE PRECISION (A-H, O-Z) ! INCLUDE 'ISRPIA.INC' ! ! DOUBLE PRECISION LAMDA, KAPA ! COMMON /CASEJ/ CHI1, CHI2, CHI3, LAMDA, KAPA, PSI1, PSI2, PSI3, & ! A1, A2, A3 ! ! *** SETUP PARAMETERS ************************************************ ! FRST = .TRUE. CALAIN = .TRUE. ! LAMDA = MAX(W(2) - W(3) - W(1), TINY) ! FREE H2SO4 PSI1 = P1 PSI2 = CHI2 ! ALL NH4HSO4 DELIQUESCED ! ! *** SOLVE EQUATIONS ; WITH ITERATIONS FOR ACTIVITY COEF. ************ ! DO 10 I=1,NSWEEP ! A1 = XK11 *(WATER/GAMA(12))**2.0 A3 = XK1 *WATER/GAMA(7)*(GAMA(8)/GAMA(7))**2.0 ! ! CALCULATE DISSOCIATION QUANTITIES ! BB = A3+LAMDA ! KAPA CC =-A3*(LAMDA + PSI1 + PSI2) DD = BB*BB-4.D0*CC KAPA = 0.5D0*(-BB+SQRT(DD)) ! ! *** CALCULATE SPECIATION ******************************************** ! MOLAL (1) = LAMDA + KAPA ! HI MOLAL (2) = PSI1 ! NAI MOLAL (3) = PSI2 ! NH4I MOLAL (4) = ZERO ! CLI MOLAL (5) = KAPA ! SO4I MOLAL (6) = LAMDA + PSI1 + PSI2 - KAPA ! HSO4I MOLAL (7) = ZERO ! NO3I ! CNAHSO4 = MAX(CHI1-PSI1,ZERO) CNH4HS4 = ZERO ! CALL CALCMR ! Water content ! ! *** CALCULATE ACTIVITIES OR TERMINATE INTERNAL LOOP ***************** ! IF (FRST.AND.CALAOU .OR. .NOT.FRST.AND.CALAIN) THEN CALL CALCACT ELSE GOTO 20 ENDIF 10 CONTINUE ! ! *** CALCULATE OBJECTIVE FUNCTION ************************************ ! 20 FUNCJ2 = MOLAL(2)*MOLAL(6)/A1 - ONE ! ! *** END OF FUNCTION FUNCJ2 ******************************************* ! END FUNCTION FUNCJ2 !======================================================================= ! ! *** ISORROPIA CODE ! *** SUBROUTINE CALCJ1 ! *** CASE J1 ! ! THE MAIN CHARACTERISTICS OF THIS REGIME ARE: ! 1. SULFATE RICH, FREE ACID (SULRAT < 1.0) ! 2. THERE IS BOTH A LIQUID & SOLID PHASE ! 3. SOLIDS POSSIBLE : NH4HSO4, NAHSO4 ! ! *** COPYRIGHT 1996-2006, UNIVERSITY OF MIAMI, CARNEGIE MELLON UNIVERSITY, ! *** GEORGIA INSTITUTE OF TECHNOLOGY ! *** WRITTEN BY ATHANASIOS NENES ! *** UPDATED BY CHRISTOS FOUNTOUKIS ! !======================================================================= ! SUBROUTINE CALCJ1 USE ISRPIA USE CASEJ IMPLICIT DOUBLE PRECISION (A-H, O-Z) ! INCLUDE 'ISRPIA.INC' ! ! DOUBLE PRECISION LAMDA, KAPA ! COMMON /CASEJ/ CHI1, CHI2, CHI3, LAMDA, KAPA, PSI1, PSI2, PSI3, & ! A1, A2, A3 ! ! *** SETUP PARAMETERS ************************************************ ! CALAOU =.TRUE. ! Outer loop activity calculation flag CHI1 = W(1) ! Total NA initially as NaHSO4 CHI2 = W(3) ! Total NH4 initially as NH4HSO4 ! PSI1LO = TINY ! Low limit PSI1HI = CHI1 ! High limit ! ! *** INITIAL VALUES FOR BISECTION ************************************ ! X1 = PSI1HI Y1 = FUNCJ1 (X1) YHI= Y1 ! Save Y-value at HI position ! ! *** YHI < 0.0 THE SOLUTION IS ALWAYS UNDERSATURATED WITH NH42SO4 **** ! IF (ABS(Y1).LE.EPS .OR. YHI.LT.ZERO) GOTO 50 ! ! *** ROOT TRACKING ; FOR THE RANGE OF HI AND LO ********************** ! DX = (PSI1HI-PSI1LO)/REAL(NDIV) DO 10 I=1,NDIV X2 = X1-DX Y2 = FUNCJ1 (X2) IF (SIGN(1.d0,Y1)*SIGN(1.d0,Y2).LT.ZERO) GOTO 20 ! (Y1*Y2.LT.ZERO) X1 = X2 Y1 = Y2 10 CONTINUE ! ! *** { YLO, YHI } > 0.0 THE SOLUTION IS ALWAYS SUPERSATURATED WITH NH42SO4 ! YLO= Y1 ! Save Y-value at Hi position IF (YLO.GT.ZERO .AND. YHI.GT.ZERO) THEN Y3 = FUNCJ1 (ZERO) GOTO 50 ELSE IF (ABS(Y2) .LT. EPS) THEN ! X2 IS A SOLUTION GOTO 50 ELSE CALL PUSHERR (0001, 'CALCJ1') ! WARNING ERROR: NO SOLUTION GOTO 50 ENDIF ! ! *** PERFORM BISECTION *********************************************** ! 20 DO 30 I=1,MAXIT X3 = 0.5*(X1+X2) Y3 = FUNCJ1 (X3) IF (SIGN(1.d0,Y1)*SIGN(1.d0,Y3) .LE. ZERO) THEN ! (Y1*Y3 .LE. ZERO) Y2 = Y3 X2 = X3 ELSE Y1 = Y3 X1 = X3 ENDIF IF (ABS(X2-X1) .LE. EPS*X1) GOTO 40 30 CONTINUE CALL PUSHERR (0002, 'CALCJ1') ! WARNING ERROR: NO CONVERGENCE ! ! *** CONVERGED ; RETURN ********************************************** ! 40 X3 = 0.5*(X1+X2) Y3 = FUNCJ1 (X3) ! 50 RETURN ! ! *** END OF SUBROUTINE CALCJ1 ****************************************** ! END SUBROUTINE CALCJ1 !======================================================================= ! ! *** ISORROPIA CODE ! *** SUBROUTINE FUNCJ1 ! *** CASE J1 ! ! THE MAIN CHARACTERISTICS OF THIS REGIME ARE: ! 1. SULFATE RICH, FREE ACID (SULRAT < 1.0) ! 2. THERE IS BOTH A LIQUID & SOLID PHASE ! 3. SOLIDS POSSIBLE : (NH4)2SO4, NH4CL, NA2SO4 ! ! *** COPYRIGHT 1996-2006, UNIVERSITY OF MIAMI, CARNEGIE MELLON UNIVERSITY, ! *** GEORGIA INSTITUTE OF TECHNOLOGY ! *** WRITTEN BY ATHANASIOS NENES ! *** UPDATED BY CHRISTOS FOUNTOUKIS ! !======================================================================= ! DOUBLE PRECISION FUNCTION FUNCJ1 (P1) USE ISRPIA USE CASEJ IMPLICIT DOUBLE PRECISION (A-H, O-Z) ! INCLUDE 'ISRPIA.INC' ! DOUBLE PRECISION LAMDA, KAPA ! COMMON /CASEJ/ CHI1, CHI2, CHI3, LAMDA, KAPA, PSI1, PSI2, PSI3, & ! A1, A2, A3 ! ! *** SETUP PARAMETERS ************************************************ ! FRST = .TRUE. CALAIN = .TRUE. ! LAMDA = MAX(W(2) - W(3) - W(1), TINY) ! FREE H2SO4 PSI1 = P1 ! ! *** SOLVE EQUATIONS ; WITH ITERATIONS FOR ACTIVITY COEF. ************ ! DO 10 I=1,NSWEEP ! A1 = XK11 *(WATER/GAMA(12))**2.0 A2 = XK12 *(WATER/GAMA(09))**2.0 A3 = XK1 *WATER/GAMA(7)*(GAMA(8)/GAMA(7))**2.0 ! PSI2 = 0.5*(-(LAMDA+PSI1) + SQRT((LAMDA+PSI1)**2.D0+4.D0*A2)) ! PSI2 PSI2 = MIN (PSI2, CHI2) ! BB = A3+LAMDA ! KAPA CC =-A3*(LAMDA + PSI2 + PSI1) DD = BB*BB-4.D0*CC KAPA = 0.5D0*(-BB+SQRT(DD)) ! ! *** SAVE CONCENTRATIONS IN MOLAL ARRAY ****************************** ! MOLAL (1) = LAMDA + KAPA ! HI MOLAL (2) = PSI1 ! NAI MOLAL (3) = PSI2 ! NH4I MOLAL (4) = ZERO MOLAL (5) = KAPA ! SO4I MOLAL (6) = LAMDA + PSI1 + PSI2 - KAPA ! HSO4I MOLAL (7) = ZERO ! CNAHSO4 = MAX(CHI1-PSI1,ZERO) CNH4HS4 = MAX(CHI2-PSI2,ZERO) ! CALL CALCMR ! Water content ! ! *** CALCULATE ACTIVITIES OR TERMINATE INTERNAL LOOP ***************** ! IF (FRST.AND.CALAOU .OR. .NOT.FRST.AND.CALAIN) THEN CALL CALCACT ELSE GOTO 20 ENDIF 10 CONTINUE ! ! *** CALCULATE OBJECTIVE FUNCTION ************************************ ! 20 FUNCJ1 = MOLAL(2)*MOLAL(6)/A1 - ONE ! ! *** END OF FUNCTION FUNCJ1 ******************************************* ! END FUNCTION FUNCJ1