From 28a379304e6958be68e7bcf629819b7b878d6116 Mon Sep 17 00:00:00 2001 From: Philippe WAUTELET <philippe.wautelet@aero.obs-mip.fr> Date: Tue, 13 Feb 2018 16:56:36 +0100 Subject: [PATCH] Philippe 13/02/2018: use ifdef MNH_REAL to prevent problems with intrinsics on Blue Gene/Q --- src/MNH/ares.f | 4 + src/MNH/bl89.f90 | 5 + src/MNH/ch_f77.fx90 | 20 +++ src/MNH/isocom.f | 39 +++-- src/MNH/isofwd.f | 39 +++-- src/MNH/isorev.f | 250 +++++++++++++++--------------- src/MNH/mode_RBK90_Integrator.f90 | 4 + src/include/isrpia.inc | 4 +- 8 files changed, 215 insertions(+), 150 deletions(-) diff --git a/src/MNH/ares.f b/src/MNH/ares.f index dc5770c78..68ac0c894 100644 --- a/src/MNH/ares.f +++ b/src/MNH/ares.f @@ -1122,7 +1122,11 @@ C NOW HERE WE HAVE ONLY ONE REAL ROOT part2=abs(rr) part3=(part1+part2)**one3rd crutes(1) = +#if (MNH_REAL == 8) + & -sign(1.0E0,rr) * ( part3 + (qq/part3) ) - a2/3. +#else & -sign(1.0D0,rr) * ( part3 + (qq/part3) ) - a2/3. +#endif crutes(2)=0. crutes(3)=0. nr=1 diff --git a/src/MNH/bl89.f90 b/src/MNH/bl89.f90 index c710f8c34..28225f342 100644 --- a/src/MNH/bl89.f90 +++ b/src/MNH/bl89.f90 @@ -293,8 +293,13 @@ DO JK=IKTB,IKTE !* 7. final mixing length ! DO J1D=1,IIU*IJU +#if (MNH_REAL == 8) + ZLWORK1=MAX(ZLMDN(J1D,JK),1.E-10) + ZLWORK2=MAX(ZLWORK(J1D),1.E-10) +#else ZLWORK1=MAX(ZLMDN(J1D,JK),1.D-10) ZLWORK2=MAX(ZLWORK(J1D),1.D-10) +#endif ZPOTE = ZLWORK1 / ZLWORK2 ZLWORK2=1.d0 + ZPOTE**(2./3.) ZLM(J1D,JK) = Z2SQRT2*ZLWORK1/(ZLWORK2*SQRT(ZLWORK2)) diff --git a/src/MNH/ch_f77.fx90 b/src/MNH/ch_f77.fx90 index 8f9ba7acb..08341fb19 100644 --- a/src/MNH/ch_f77.fx90 +++ b/src/MNH/ch_f77.fx90 @@ -8652,13 +8652,21 @@ c ** 51545.0 + 2.4e6 = noon 1 jan 2000 c ** force mean longitude between 0 and 360 degs MNLONG = 280.460 + 0.9856474*TIME +#if (MNH_REAL == 8) + MNLONG = MOD( MNLONG, 360.E0 ) +#else MNLONG = MOD( MNLONG, 360.D0 ) +#endif IF( MNLONG.LT.0. ) MNLONG = MNLONG + 360. c ** mean anomaly in radians between 0 and 2*pi MNANOM = 357.528 + 0.9856003*TIME +#if (MNH_REAL == 8) + MNANOM = MOD( MNANOM, 360.E0 ) +#else MNANOM = MOD( MNANOM, 360.D0 ) +#endif IF( MNANOM.LT.0.) MNANOM = MNANOM + 360. MNANOM = MNANOM*RPD @@ -8667,7 +8675,11 @@ c ** ecliptic longitude and obliquity c ** of ecliptic in radians ECLONG = MNLONG + 1.915*SIN( MNANOM ) + 0.020*SIN( 2.*MNANOM ) +#if (MNH_REAL == 8) + ECLONG = MOD( ECLONG, 360.E0 ) +#else ECLONG = MOD( ECLONG, 360.D0 ) +#endif IF( ECLONG.LT.0. ) ECLONG = ECLONG + 360. OBLQEC = 23.439 - 0.0000004*TIME @@ -8699,13 +8711,21 @@ c ** Greenwich mean sidereal time in hours c ** Hour not changed to sidereal time since c ** 'time' includes the fractional day +#if (MNH_REAL == 8) + GMST = MOD( GMST, 24.E0) +#else GMST = MOD( GMST, 24.D0) +#endif IF( GMST.LT.0. ) GMST = GMST + 24. c ** local mean sidereal time in radians LMST = GMST + LONG / 15. +#if (MNH_REAL == 8) + LMST = MOD( LMST, 24.E0 ) +#else LMST = MOD( LMST, 24.D0 ) +#endif IF( LMST.LT.0. ) LMST = LMST + 24. LMST = LMST*15.*RPD diff --git a/src/MNH/isocom.f b/src/MNH/isocom.f index 27049a2a5..814eedca4 100644 --- a/src/MNH/isocom.f +++ b/src/MNH/isocom.f @@ -374,14 +374,25 @@ C C C *** DEFAULT VALUES ************************************************* C +#if (MNH_REAL == 8) + DATA TEMP/298.E0/, R/82.0567E-6/, RH/0.9E0/, EPS/1E-6/, MAXIT/100/, + & TINY/1E-20/, GREAT/1E10/, ZERO/0.0E0/, ONE/1.0E0/,NSWEEP/4/, + & TINY2/1E-11/,NDIV/5/, TWO/2.0E0/, HALF/0.5E0/, FOUR/4.0E0/ +C + DATA MOLAL/NIONS*0.0E0/, MOLALR/NPAIR*0.0E0/, GAMA/NPAIR*0.1E0/, + & GAMOU/NPAIR*1E10/, GAMIN/NPAIR*1E10/, CALAIN/.TRUE./, + & CALAOU/.TRUE./, EPSACT/5E-2/, ICLACT/0/, + & IACALC/1/, WFTYP/2/ +#else DATA TEMP/298.0/, R/82.0567D-6/, RH/0.9D0/, EPS/1D-6/, MAXIT/100/, & TINY/1D-20/, GREAT/1D10/, ZERO/0.0D0/, ONE/1.0D0/,NSWEEP/4/, - & TINY2/1D-11/,NDIV/5/ + & TINY2/1D-11/,NDIV/5/, TWO/2.0D0/, HALF/0.5E0/, FOUR/4.0E0/ C DATA MOLAL/NIONS*0.0D0/, MOLALR/NPAIR*0.0D0/, GAMA/NPAIR*0.1D0/, & GAMOU/NPAIR*1D10/, GAMIN/NPAIR*1D10/, CALAIN/.TRUE./, & CALAOU/.TRUE./, EPSACT/5D-2/, ICLACT/0/, & IACALC/1/, WFTYP/2/ +#endif C DATA ERRSTK/NERRMX*0/, ERRMSG/NERRMX*' '/, NOFER/0/, & STKOFL/.FALSE./ @@ -1543,7 +1554,7 @@ CC ENDIF C C *** CALCULATE HCL SPECIATION IN THE GAS PHASE ************************* C - GHCL = MAX(X-DELT, 0.0d0) ! GAS HCL + GHCL = MAX(X-DELT, ZERO) ! GAS HCL C C *** CALCULATE HCL SPECIATION IN THE LIQUID PHASE ********************** C @@ -1640,7 +1651,7 @@ CC ENDIF C C *** CALCULATE HNO3 SPECIATION IN THE GAS PHASE ************************ C - GHNO3 = MAX(X-DELT, 0.0d0) ! GAS HNO3 + GHNO3 = MAX(X-DELT, ZERO) ! GAS HNO3 C C *** CALCULATE HNO3 SPECIATION IN THE LIQUID PHASE ********************* C @@ -2323,7 +2334,7 @@ C MOLALR(2) = 0.5*MOLAL(2) ! NA2SO4 TOTS4 = MOLAL(5)+MOLAL(6) ! Total SO4 MOLALR(4) = MAX(TOTS4 - MOLALR(2), ZERO) ! (NH4)2SO4 - FRNH4 = MAX(MOLAL(3) - 2.D0*MOLALR(4), ZERO) + FRNH4 = MAX(MOLAL(3) - TWO*MOLALR(4), ZERO) MOLALR(5) = MIN(MOLAL(7),FRNH4) ! NH4NO3 FRNH4 = MAX(FRNH4 - MOLALR(5), ZERO) MOLALR(6) = MIN(MOLAL(4), FRNH4) ! NH4CL @@ -2774,13 +2785,13 @@ C BB =-GG CC =-AKW DD = BB*BB - 4.D0*CC - HI = MAX(0.5D0*(-BB + SQRT(DD)),CN) + HI = MAX(HALF*(-BB + SQRT(DD)),CN) OHI= AKW/HI ELSE ! OH- in excess BB = GG CC =-AKW DD = BB*BB - 4.D0*CC - OHI= MAX(0.5D0*(-BB + SQRT(DD)),CN) + OHI= MAX(HALF*(-BB + SQRT(DD)),CN) HI = AKW/OHI ENDIF C @@ -2831,7 +2842,11 @@ C DO 30 I=1,NIONS IONIC=IONIC + MOLAL(I)*Z(I)*Z(I) 30 CONTINUE - IONIC = MAX(MIN(0.5*IONIC/WATER,20.d0), TINY) +#if (MNH_REAL == 8) + IONIC = MAX(MIN(HALF*IONIC/WATER,20.e0), TINY) +#else + IONIC = MAX(MIN(HALF*IONIC/WATER,20.d0), TINY) +#endif C C *** CALCULATE BINARY ACTIVITY COEFFICIENTS *************************** C @@ -3811,7 +3826,7 @@ C D = 0, THREE REAL (ONE DOUBLE) ROOTS C ELSE IF (D.LE.EPS) THEN ! -EPS <= D <= EPS : D = ZERO IX = 2 - SSIG = SIGN (1.D0, R) + SSIG = SIGN (ONE, R) S = SSIG*(ABS(R))**EXPON X(1) = 2.D0*S - EXPON*A1 X(2) = -S - EXPON*A1 @@ -3821,8 +3836,8 @@ C ELSE ! D > EPS : D > ZERO IX = 1 SQD = SQRT(D) - SSIG = SIGN (1.D0, R+SQD) ! TRANSFER SIGN TO SSIG - TSIG = SIGN (1.D0, R-SQD) + SSIG = SIGN (ONE, R+SQD) ! TRANSFER SIGN TO SSIG + TSIG = SIGN (ONE, R-SQD) S = SSIG*(ABS(R+SQD))**EXPON ! EXPONENTIATE ABS() T = TSIG*(ABS(R-SQD))**EXPON X(1) = S + T - EXPON*A1 @@ -3886,7 +3901,7 @@ C DO 10 I=1,NDIV X2 = X1+DX Y2 = FUNC (X2) - IF (SIGN(1.d0,Y1)*SIGN(1.d0,Y2) .LT. ZERO) GOTO 20 ! (Y1*Y2.LT.ZERO) + IF (SIGN(ONE,Y1)*SIGN(ONE,Y2) .LT. ZERO) GOTO 20 ! (Y1*Y2.LT.ZERO) X1 = X2 Y1 = Y2 10 CONTINUE @@ -3906,7 +3921,7 @@ C 20 DO 30 I=1,MAXIT X3 = 0.5*(X1+X2) Y3 = FUNC (X3) - IF (SIGN(1.d0,Y1)*SIGN(1.d0,Y3) .LE. ZERO) THEN ! (Y1*Y3 .LE. ZERO) + IF (SIGN(ONE,Y1)*SIGN(ONE,Y3) .LE. ZERO) THEN ! (Y1*Y3 .LE. ZERO) Y2 = Y3 X2 = X3 ELSE diff --git a/src/MNH/isofwd.f b/src/MNH/isofwd.f index 8e030d139..1fabfbf06 100644 --- a/src/MNH/isofwd.f +++ b/src/MNH/isofwd.f @@ -260,8 +260,13 @@ C C C *** ADJUST FOR TOO LITTLE AMMONIUM AND CHLORIDE *********************** C +#if (MNH_REAL == 8) + WI(3) = MAX (WI(3), 1.E-10) ! NH4+ : 1e-4 umoles/m3 + WI(5) = MAX (WI(5), 1.E-10) ! Cl- : 1e-4 umoles/m3 +#else 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 +#endif C C *** ADJUST FOR TOO LITTLE SODIUM, SULFATE AND NITRATE COMBINED ******** C @@ -2189,7 +2194,7 @@ C CHI1 = 0.5*W(1) CHI2 = MAX (W(2)-CHI1, ZERO) CHI3 = ZERO - CHI4 = MAX (W(3)-2.D0*CHI2, ZERO) + CHI4 = MAX (W(3)-TWO*CHI2, ZERO) CHI5 = W(4) CHI6 = W(5) C @@ -2318,7 +2323,7 @@ CCC IF(CHI4.GT.TINY) THEN IF(W(2).GT.TINY) THEN ! Accounts for NH3 evaporation BB =-(CHI4 + PSI6 + PSI5 + 1./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 + DD = MAX(BB*BB-FOUR*CC,ZERO) ! Patch proposed by Uma Shankar, 19/11/01 PSI4 =0.5d0*(-BB - SQRT(DD)) ELSE PSI4 = TINY @@ -2397,7 +2402,7 @@ C CHI1 = 0.5*W(1) CHI2 = MAX (W(2)-CHI1, ZERO) CHI3 = ZERO - CHI4 = MAX (W(3)-2.D0*CHI2, ZERO) + CHI4 = MAX (W(3)-TWO*CHI2, ZERO) CHI5 = W(4) CHI6 = W(5) C @@ -2525,7 +2530,7 @@ CCC IF(CHI4.GT.TINY) THEN IF(W(2).GT.TINY) THEN ! Accounts for NH3 evaporation BB =-(CHI4 + PSI6 + PSI5 + 1./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 + DD = MAX(BB*BB-FOUR*CC,ZERO) ! Patch proposed by Uma shankar, 19/11/2001 PSI4 =0.5d0*(-BB - SQRT(DD)) ELSE PSI4 = TINY @@ -2683,7 +2688,7 @@ C CHI1 = 0.5*W(1) CHI2 = MAX (W(2)-CHI1, ZERO) CHI3 = ZERO - CHI4 = MAX (W(3)-2.D0*CHI2, ZERO) + CHI4 = MAX (W(3)-TWO*CHI2, ZERO) CHI5 = W(4) CHI6 = W(5) C @@ -2830,7 +2835,7 @@ CCC IF(CHI4.GT.TINY) THEN IF(W(2).GT.TINY) THEN ! Accounts for NH3 evaporation BB =-(CHI4 + PSI6 + PSI5 + 1./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 + DD = MAX(BB*BB-FOUR*CC,ZERO) ! Patch proposed by Uma Shankar, 19/11/01 PSI4 =0.5d0*(-BB - SQRT(DD)) ELSE PSI4 = TINY @@ -2984,7 +2989,7 @@ C CHI1 = 0.5*W(1) CHI2 = MAX (W(2)-CHI1, ZERO) CHI3 = ZERO - CHI4 = MAX (W(3)-2.D0*CHI2, ZERO) + CHI4 = MAX (W(3)-TWO*CHI2, ZERO) CHI5 = W(4) CHI6 = W(5) C @@ -3422,7 +3427,7 @@ C CHI1 = W(2) ! CNA2SO4 CHI2 = ZERO ! CNH42S4 CHI3 = ZERO ! CNH4CL - FRNA = MAX (W(1)-2.D0*CHI1, ZERO) + FRNA = MAX (W(1)-TWO*CHI1, ZERO) CHI8 = MIN (FRNA, W(4)) ! CNANO3 CHI4 = W(3) ! NH3(g) CHI5 = MAX (W(4)-CHI8, ZERO) ! HNO3(g) @@ -3639,7 +3644,7 @@ C CHI1 = W(2) ! CNA2SO4 CHI2 = ZERO ! CNH42S4 CHI3 = ZERO ! CNH4CL - FRNA = MAX (W(1)-2.D0*CHI1, ZERO) + FRNA = MAX (W(1)-TWO*CHI1, ZERO) CHI8 = MIN (FRNA, W(4)) ! CNANO3 CHI4 = W(3) ! NH3(g) CHI5 = MAX (W(4)-CHI8, ZERO) ! HNO3(g) @@ -3868,7 +3873,7 @@ C CHI1 = W(2) ! CNA2SO4 CHI2 = ZERO ! CNH42S4 CHI3 = ZERO ! CNH4CL - FRNA = MAX (W(1)-2.D0*CHI1, ZERO) + FRNA = MAX (W(1)-TWO*CHI1, ZERO) CHI8 = MIN (FRNA, W(4)) ! CNANO3 CHI4 = W(3) ! NH3(g) CHI5 = MAX (W(4)-CHI8, ZERO) ! HNO3(g) @@ -4122,7 +4127,7 @@ C CHI1 = W(2) ! CNA2SO4 CHI2 = ZERO ! CNH42S4 CHI3 = ZERO ! CNH4CL - FRNA = MAX (W(1)-2.D0*CHI1, ZERO) + FRNA = MAX (W(1)-TWO*CHI1, ZERO) CHI8 = MIN (FRNA, W(4)) ! CNANO3 CHI4 = W(3) ! NH3(g) CHI5 = MAX (W(4)-CHI8, ZERO) ! HNO3(g) @@ -4431,7 +4436,7 @@ C CHI1 = W(2) ! CNA2SO4 CHI2 = ZERO ! CNH42S4 CHI3 = ZERO ! CNH4CL - FRNA = MAX (W(1)-2.D0*CHI1, ZERO) + FRNA = MAX (W(1)-TWO*CHI1, ZERO) CHI8 = MIN (FRNA, W(4)) ! CNANO3 CHI4 = W(3) ! NH3(g) CHI5 = MAX (W(4)-CHI8, ZERO) ! HNO3(g) @@ -6009,9 +6014,15 @@ C CNH42S4 = ZERO FRSO4 = MAX(W(2)-CNA2SO4, ZERO) C +#if (MNH_REAL == 8) + CLC = MIN(W(3)/3.E0, FRSO4/2.E0) + FRSO4 = MAX(FRSO4-2.E0*CLC, ZERO) + FRNH4 = MAX(W(3)-3.E0*CLC, ZERO) +#else CLC = MIN(W(3)/3.D0, FRSO4/2.D0) FRSO4 = MAX(FRSO4-2.D0*CLC, ZERO) FRNH4 = MAX(W(3)-3.D0*CLC, ZERO) +#endif C IF (FRSO4.LE.TINY) THEN CLC = MAX(CLC - FRNH4, ZERO) @@ -6021,7 +6032,11 @@ C CNH4HS4 = 3.D0*MIN(FRSO4, CLC) CLC = MAX(CLC-FRSO4, ZERO) IF (CNA2SO4.GT.TINY) THEN +#if (MNH_REAL == 8) + FRSO4 = MAX(FRSO4-CNH4HS4/3.E0, ZERO) +#else FRSO4 = MAX(FRSO4-CNH4HS4/3.D0, ZERO) +#endif CNAHSO4 = 2.D0*FRSO4 CNA2SO4 = MAX(CNA2SO4-FRSO4, ZERO) ENDIF diff --git a/src/MNH/isorev.f b/src/MNH/isorev.f index f4588f11c..f774dcd20 100644 --- a/src/MNH/isorev.f +++ b/src/MNH/isorev.f @@ -36,7 +36,7 @@ C IF (RH.GE.DRNH42S4) THEN ! WET AEROSOL, NEED NH4 AT SRATIO=2.0 SULRATW = GETASR(WAER(2), RHI) ! AEROSOL SULFATE RATIO ELSE - SULRATW = 2.0D0 ! DRY AEROSOL SULFATE RATIO + SULRATW = TWO ! DRY AEROSOL SULFATE RATIO ENDIF SULRAT = WAER(3)/WAER(2) ! SULFATE RATIO C @@ -160,7 +160,7 @@ C IF (TRYLIQ .AND. RH.GE.DRNH4NO3) THEN ! *** WET AEROSOL SULRATW = GETASR(WAER(2), RHI) ! LIMITING SULFATE RATIO ELSE - SULRATW = 2.0D0 ! *** DRY AEROSOL + SULRATW = TWO ! *** DRY AEROSOL ENDIF SULRAT = WAER(3)/WAER(2) C @@ -309,7 +309,7 @@ C ccC ccC *** CHECK IF TOO MUCH SODIUM ; ADJUST AND ISSUE ERROR MESSAGE ********* ccC -cc REST = 2.D0*WAER(2) + WAER(4) + WAER(5) +cc REST = TWO*WAER(2) + WAER(4) + WAER(5) cc IF (WAER(1).GT.REST) THEN ! NA > 2*SO4+CL+NO3 ? cc WAER(1) = (ONE-1D-6)*REST ! Adjust Na amount cc CALL PUSHERR (0050, 'ISRP3R') ! Warning error: Na adjusted @@ -318,13 +318,13 @@ C C *** CALCULATE SULFATE & SODIUM RATIOS ********************************* C IF (TRYLIQ .AND. RH.GE.DRNH4NO3) THEN ! ** WET AEROSOL - FRSO4 = WAER(2) - WAER(1)/2.0D0 ! SULFATE UNBOUND BY SODIUM + FRSO4 = WAER(2) - WAER(1)/TWO ! SULFATE UNBOUND BY SODIUM FRSO4 = MAX(FRSO4, TINY) SRI = GETASR(FRSO4, RHI) ! SULFATE RATIO FOR NH4+ SULRATW = (WAER(1)+FRSO4*SRI)/WAER(2) ! LIMITING SULFATE RATIO - SULRATW = MIN (SULRATW, 2.0D0) + SULRATW = MIN (SULRATW, TWO) ELSE - SULRATW = 2.0D0 ! ** DRY AEROSOL + SULRATW = TWO ! ** DRY AEROSOL ENDIF SULRAT = (WAER(1)+WAER(3))/WAER(2) SODRAT = WAER(1)/WAER(2) @@ -527,7 +527,7 @@ C C C *** CALCULATE WATER CONTENT ***************************************** C - MOLALR(4)= MIN(WAER(2), 0.5d0*WAER(3)) + MOLALR(4)= MIN(WAER(2), HALF*WAER(3)) WATER = MOLALR(4)/M0(4) ! ZSR correlation C C *** SOLVE EQUATIONS ; WITH ITERATIONS FOR ACTIVITY COEF. ************ @@ -541,7 +541,7 @@ C SO4I = WAER(2) HSO4I= ZERO C - CALL CALCPH (2.D0*SO4I - NH4I, HI, OHI) ! Get pH + CALL CALCPH (TWO*SO4I - NH4I, HI, OHI) ! Get pH C NH3AQ = ZERO ! AMMONIA EQUILIBRIUM IF (HI.LT.OHI) THEN @@ -606,11 +606,11 @@ C SUBROUTINE CALCK1 INCLUDE 'isrpia.inc' C - CNH42S4 = MIN(WAER(2),0.5d0*WAER(3)) ! For bad input problems + CNH42S4 = MIN(WAER(2),HALF*WAER(3)) ! For bad input problems GNH3 = ZERO C W(2) = CNH42S4 - W(3) = 2.D0*CNH42S4 + GNH3 + W(3) = TWO*CNH42S4 + GNH3 C RETURN C @@ -649,8 +649,8 @@ C C C *** AEROSOL WATER CONTENT C - MOLALR(4) = MIN(WAER(2),0.5d0*WAER(3)) ! (NH4)2SO4 - AML5 = MAX(WAER(3)-2.D0*MOLALR(4),ZERO) ! "free" NH4 + MOLALR(4) = MIN(WAER(2),HALF*WAER(3)) ! (NH4)2SO4 + AML5 = MAX(WAER(3)-TWO*MOLALR(4),ZERO) ! "free" NH4 MOLALR(5) = MAX(MIN(AML5,WAER(4)), ZERO) ! NH4NO3=MIN("free",NO3) WATER = MOLALR(4)/M0(4) + MOLALR(5)/M0(5) WATER = MAX(WATER, TINY) @@ -671,13 +671,13 @@ C SO4I = WAER(2) HSO4I = ZERO C - CALL CALCPH (2.D0*SO4I + NO3I - NH4I, HI, OHI) + CALL CALCPH (TWO*SO4I + NO3I - NH4I, HI, OHI) C C AMMONIA ASSOCIATION EQUILIBRIUM C NH3AQ = ZERO NO3AQ = ZERO - GG = 2.D0*SO4I + NO3I - NH4I + GG = TWO*SO4I + NO3I - NH4I IF (HI.LT.OHI) THEN CALL CALCAMAQ2 (-GG, NH4I, OHI, NH3AQ) HI = AKW/OHI @@ -752,8 +752,8 @@ C C C *** SETUP PARAMETERS ************************************************ C - CHI1 = MIN(WAER(2),0.5d0*WAER(3)) ! (NH4)2SO4 - CHI2 = MAX(WAER(3) - 2.D0*CHI1, ZERO) ! "Free" NH4+ + CHI1 = MIN(WAER(2),HALF*WAER(3)) ! (NH4)2SO4 + CHI2 = MAX(WAER(3) - TWO*CHI1, ZERO) ! "Free" NH4+ CHI3 = MAX(WAER(4) - CHI2, ZERO) ! "Free" NO3 C PSI2 = CHI2 @@ -869,18 +869,18 @@ CC A21 = XK21*WATER*R*TEMP C C ION CONCENTRATIONS C - NH4I = 2.D0*PSI1 + PSI2 + NH4I = TWO*PSI1 + PSI2 NO3I = PSI2 + PSI3 SO4I = PSI1 HSO4I = ZERO C - CALL CALCPH (2.D0*SO4I + NO3I - NH4I, HI, OHI) + CALL CALCPH (TWO*SO4I + NO3I - NH4I, HI, OHI) C C AMMONIA ASSOCIATION EQUILIBRIUM C NH3AQ = ZERO NO3AQ = ZERO - GG = 2.D0*SO4I + NO3I - NH4I + GG = TWO*SO4I + NO3I - NH4I IF (HI.LT.OHI) THEN CALL CALCAMAQ2 (-GG, NH4I, OHI, NH3AQ) HI = AKW/OHI @@ -1005,13 +1005,13 @@ CCC A1 = XK10/R/TEMP/R/TEMP C C *** CALCULATE AEROSOL COMPOSITION ************************************ C -CCC CHI1 = 2.D0*WAER(4) ! Free parameter ; arbitrary value. +CCC CHI1 = TWO*WAER(4) ! Free parameter ; arbitrary value. PSI1 = WAER(4) C C *** The following statment is here to avoid negative NH4+ values in C CALCN? routines that call CALCN1A C - PSI2 = MAX(MIN(WAER(2),0.5d0*(WAER(3)-PSI1)),TINY) + PSI2 = MAX(MIN(WAER(2),HALF*(WAER(3)-PSI1)),TINY) C CNH4NO3 = PSI1 CNH42S4 = PSI2 @@ -1090,18 +1090,18 @@ C C C SOLUTION ACIDIC OR BASIC? C - GG = 2.D0*SO4I + NO3I + CLI - NAI - NH4I + GG = TWO*SO4I + NO3I + CLI - NAI - NH4I IF (GG.GT.TINY) THEN ! H+ in excess BB =-GG CC =-AKW - DD = BB*BB - 4.D0*CC - HI = 0.5D0*(-BB + SQRT(DD)) + DD = BB*BB - FOUR*CC + HI = HALF*(-BB + SQRT(DD)) OHI= AKW/HI ELSE ! OH- in excess BB = GG CC =-AKW - DD = BB*BB - 4.D0*CC - OHI= 0.5D0*(-BB + SQRT(DD)) + DD = BB*BB - FOUR*CC + OHI= HALF*(-BB + SQRT(DD)) HI = AKW/OHI ENDIF C @@ -1112,7 +1112,7 @@ C HI = AKW/OHI HSO4I = ZERO ELSE - GGNO3 = MAX(2.D0*SO4I + NO3I - NAI - NH4I, ZERO) + GGNO3 = MAX(TWO*SO4I + NO3I - NAI - NH4I, ZERO) GGCL = MAX(GG-GGNO3, ZERO) IF (GGCL .GT.TINY) CALL CALCCLAQ2 (GGCL, CLI, HI, CLAQ) ! HCl IF (GGNO3.GT.TINY) THEN @@ -1255,7 +1255,7 @@ C C C ION CONCENTRATIONS ; CORRECTIONS C - NAI = WAER(1) - 2.D0*ROOT3 + NAI = WAER(1) - TWO*ROOT3 SO4I= WAER(2) - ROOT3 NH4I = WAER(3) NO3I = WAER(4) @@ -1263,18 +1263,18 @@ C C C SOLUTION ACIDIC OR BASIC? C - GG = 2.D0*SO4I + NO3I + CLI - NAI - NH4I + GG = TWO*SO4I + NO3I + CLI - NAI - NH4I IF (GG.GT.TINY) THEN ! H+ in excess BB =-GG CC =-AKW - DD = BB*BB - 4.D0*CC - HI = 0.5D0*(-BB + SQRT(DD)) + DD = BB*BB - FOUR*CC + HI = HALF*(-BB + SQRT(DD)) OHI= AKW/HI ELSE ! OH- in excess BB = GG CC =-AKW - DD = BB*BB - 4.D0*CC - OHI= 0.5D0*(-BB + SQRT(DD)) + DD = BB*BB - FOUR*CC + OHI= HALF*(-BB + SQRT(DD)) HI = AKW/OHI ENDIF C @@ -1285,7 +1285,7 @@ C HI = AKW/OHI HSO4I = ZERO ELSE - GGNO3 = MAX(2.D0*SO4I + NO3I - NAI - NH4I, ZERO) + GGNO3 = MAX(TWO*SO4I + NO3I - NAI - NH4I, ZERO) GGCL = MAX(GG-GGNO3, ZERO) IF (GGCL .GT.TINY) CALL CALCCLAQ2 (GGCL, CLI, HI, CLAQ) ! HCl IF (GGNO3.GT.TINY) THEN @@ -1493,8 +1493,8 @@ C AMMONIUM SULFATE C IF (NH4I*NH4I*SO4I .GT. A4) THEN BB =-(WAER(2)+WAER(3)-ROOT3) - CC = WAER(3)*(WAER(2)-ROOT3+0.5D0*WAER(3)) - DD =-((WAER(2)-ROOT3)*WAER(3)**2.D0 + A4)/4.D0 + CC = WAER(3)*(WAER(2)-ROOT3+HALF*WAER(3)) + DD =-((WAER(2)-ROOT3)*WAER(3)**TWO + A4)/FOUR CALL POLY3(BB, CC, DD, ROOT1, ISLV) IF (ISLV.NE.0) ROOT1 = TINY ROOT1 = MIN(ROOT1, WAER(3), WAER(2)-ROOT3, CHI6) @@ -1506,26 +1506,26 @@ C C C ION CONCENTRATIONS C - NAI = WAER(1) - 2.D0*ROOT3 + NAI = WAER(1) - TWO*ROOT3 SO4I= WAER(2) - ROOT1 - ROOT3 - NH4I= WAER(3) - 2.D0*ROOT1 + NH4I= WAER(3) - TWO*ROOT1 NO3I= WAER(4) CLI = WAER(5) C C SOLUTION ACIDIC OR BASIC? C - GG = 2.D0*SO4I + NO3I + CLI - NAI - NH4I + GG = TWO*SO4I + NO3I + CLI - NAI - NH4I IF (GG.GT.TINY) THEN ! H+ in excess BB =-GG CC =-AKW - DD = BB*BB - 4.D0*CC - HI = 0.5D0*(-BB + SQRT(DD)) + DD = BB*BB - FOUR*CC + HI = HALF*(-BB + SQRT(DD)) OHI= AKW/HI ELSE ! OH- in excess BB = GG CC =-AKW - DD = BB*BB - 4.D0*CC - OHI= 0.5D0*(-BB + SQRT(DD)) + DD = BB*BB - FOUR*CC + OHI= HALF*(-BB + SQRT(DD)) HI = AKW/OHI ENDIF C @@ -1536,7 +1536,7 @@ C HI = AKW/OHI HSO4I = ZERO ELSE - GGNO3 = MAX(2.D0*SO4I + NO3I - NAI - NH4I, ZERO) + GGNO3 = MAX(TWO*SO4I + NO3I - NAI - NH4I, ZERO) GGCL = MAX(GG-GGNO3, ZERO) IF (GGCL .GT.TINY) CALL CALCCLAQ2 (GGCL, CLI, HI, CLAQ) ! HCl IF (GGNO3.GT.TINY) THEN @@ -1742,21 +1742,21 @@ C C AMMONIUM CHLORIDE C IF (NH4I*CLI .GT. A14) THEN - BB =-(WAER(3) + WAER(5) - 2.D0*ROOT1) - CC = WAER(5)*(WAER(3) - 2.D0*ROOT1) - A14 - DD = BB*BB - 4.D0*CC + BB =-(WAER(3) + WAER(5) - TWO*ROOT1) + CC = WAER(5)*(WAER(3) - TWO*ROOT1) - A14 + DD = BB*BB - FOUR*CC IF (DD.LT.ZERO) THEN ROOT2 = ZERO ELSE DD = SQRT(DD) - ROOT2A= 0.5D0*(-BB+DD) - ROOT2B= 0.5D0*(-BB-DD) + ROOT2A= HALF*(-BB+DD) + ROOT2B= HALF*(-BB-DD) IF (ZERO.LE.ROOT2A) THEN ROOT2 = ROOT2A ELSE ROOT2 = ROOT2B ENDIF - ROOT2 = MIN(ROOT2, WAER(5), WAER(3) - 2.D0*ROOT1, CHI4) + ROOT2 = MIN(ROOT2, WAER(5), WAER(3) - TWO*ROOT1, CHI4) ROOT2 = MAX(ROOT2, ZERO) PSI4 = CHI4 - ROOT2 ENDIF @@ -1783,8 +1783,8 @@ C AMMONIUM SULFATE C IF (NH4I*NH4I*SO4I .GT. A4) THEN BB =-(WAER(2)+WAER(3)-ROOT2-ROOT3) - CC = (WAER(3)-ROOT2)*(WAER(2)-ROOT3+0.5D0*(WAER(3)-ROOT2)) - DD =-((WAER(2)-ROOT3)*(WAER(3)-ROOT2)**2.D0 + A4)/4.D0 + CC = (WAER(3)-ROOT2)*(WAER(2)-ROOT3+HALF*(WAER(3)-ROOT2)) + DD =-((WAER(2)-ROOT3)*(WAER(3)-ROOT2)**TWO + A4)/FOUR CALL POLY3(BB, CC, DD, ROOT1, ISLV) IF (ISLV.NE.0) ROOT1 = TINY ROOT1 = MIN(ROOT1, WAER(3)-ROOT2, WAER(2)-ROOT3, CHI6) @@ -1796,26 +1796,26 @@ C C C ION CONCENTRATIONS C - NAI = WAER(1) - 2.D0*ROOT3 + NAI = WAER(1) - TWO*ROOT3 SO4I= WAER(2) - ROOT1 - ROOT3 - NH4I= WAER(3) - ROOT2 - 2.D0*ROOT1 + NH4I= WAER(3) - ROOT2 - TWO*ROOT1 NO3I= WAER(4) CLI = WAER(5) - ROOT2 C C SOLUTION ACIDIC OR BASIC? C - GG = 2.D0*SO4I + NO3I + CLI - NAI - NH4I + GG = TWO*SO4I + NO3I + CLI - NAI - NH4I IF (GG.GT.TINY) THEN ! H+ in excess BB =-GG CC =-AKW - DD = BB*BB - 4.D0*CC - HI = 0.5D0*(-BB + SQRT(DD)) + DD = BB*BB - FOUR*CC + HI = HALF*(-BB + SQRT(DD)) OHI= AKW/HI ELSE ! OH- in excess BB = GG CC =-AKW - DD = BB*BB - 4.D0*CC - OHI= 0.5D0*(-BB + SQRT(DD)) + DD = BB*BB - FOUR*CC + OHI= HALF*(-BB + SQRT(DD)) HI = AKW/OHI ENDIF C @@ -1826,7 +1826,7 @@ C HI = AKW/OHI HSO4I = ZERO ELSE - GGNO3 = MAX(2.D0*SO4I + NO3I - NAI - NH4I, ZERO) + GGNO3 = MAX(TWO*SO4I + NO3I - NAI - NH4I, ZERO) GGCL = MAX(GG-GGNO3, ZERO) IF (GGCL .GT.TINY) CALL CALCCLAQ2 (GGCL, CLI, HI, CLAQ) ! HCl IF (GGNO3.GT.TINY) THEN @@ -1992,11 +1992,11 @@ C C C *** CALCULATE SOLIDS ************************************************** C - CNA2SO4 = 0.5d0*WAER(1) + CNA2SO4 = HALF*WAER(1) FRSO4 = MAX (WAER(2)-CNA2SO4, ZERO) C - CNH42S4 = MAX (MIN(FRSO4,0.5d0*WAER(3)), TINY) - FRNH3 = MAX (WAER(3)-2.D0*CNH42S4, ZERO) + CNH42S4 = MAX (MIN(FRSO4,HALF*WAER(3)), TINY) + FRNH3 = MAX (WAER(3)-TWO*CNH42S4, ZERO) C CNH4NO3 = MIN (FRNH3, WAER(4)) CCC FRNO3 = MAX (WAER(4)-CNH4NO3, ZERO) @@ -2080,18 +2080,18 @@ C C C SOLUTION ACIDIC OR BASIC? C - GG = 2.D0*WAER(2) + NO3I + CLI - NAI - NH4I + GG = TWO*WAER(2) + NO3I + CLI - NAI - NH4I IF (GG.GT.TINY) THEN ! H+ in excess BB =-GG CC =-AKW - DD = BB*BB - 4.D0*CC - HI = 0.5D0*(-BB + SQRT(DD)) + DD = BB*BB - FOUR*CC + HI = HALF*(-BB + SQRT(DD)) OHI= AKW/HI ELSE ! OH- in excess BB = GG CC =-AKW - DD = BB*BB - 4.D0*CC - OHI= 0.5D0*(-BB + SQRT(DD)) + DD = BB*BB - FOUR*CC + OHI= HALF*(-BB + SQRT(DD)) HI = AKW/OHI ENDIF C @@ -2101,7 +2101,7 @@ C CALL CALCAMAQ2 (-GG, NH4I, OHI, NH3AQ) HI = AKW/OHI ELSE - GGNO3 = MAX(2.D0*SO4I + NO3I - NAI - NH4I, ZERO) + GGNO3 = MAX(TWO*SO4I + NO3I - NAI - NH4I, ZERO) GGCL = MAX(GG-GGNO3, ZERO) IF (GGCL .GT.TINY) CALL CALCCLAQ2 (GGCL, CLI, HI, CLAQ) ! HCl IF (GGNO3.GT.TINY) THEN @@ -2254,7 +2254,7 @@ C C C ION CONCENTRATIONS C - NAI = WAER(1) - 2.D0*ROOT + NAI = WAER(1) - TWO*ROOT SO4I = WAER(2) - ROOT NH4I = WAER(3) NO3I = WAER(4) @@ -2262,18 +2262,18 @@ C C C SOLUTION ACIDIC OR BASIC? C - GG = 2.D0*SO4I + NO3I + CLI - NAI - NH4I + GG = TWO*SO4I + NO3I + CLI - NAI - NH4I IF (GG.GT.TINY) THEN ! H+ in excess BB =-GG CC =-AKW - DD = BB*BB - 4.D0*CC - HI = 0.5D0*(-BB + SQRT(DD)) + DD = BB*BB - FOUR*CC + HI = HALF*(-BB + SQRT(DD)) OHI= AKW/HI ELSE ! OH- in excess BB = GG CC =-AKW - DD = BB*BB - 4.D0*CC - OHI= 0.5D0*(-BB + SQRT(DD)) + DD = BB*BB - FOUR*CC + OHI= HALF*(-BB + SQRT(DD)) HI = AKW/OHI ENDIF C @@ -2283,7 +2283,7 @@ C CALL CALCAMAQ2 (-GG, NH4I, OHI, NH3AQ) HI = AKW/OHI ELSE - GGNO3 = MAX(2.D0*SO4I + NO3I - NAI - NH4I, ZERO) + GGNO3 = MAX(TWO*SO4I + NO3I - NAI - NH4I, ZERO) GGCL = MAX(GG-GGNO3, ZERO) IF (GGCL .GT.TINY) CALL CALCCLAQ2 (GGCL, CLI, HI, CLAQ) ! HCl IF (GGNO3.GT.TINY) THEN @@ -2485,7 +2485,7 @@ C IF (ISLV.NE.0) ROOT = TINY ROOT = MIN (MAX(ROOT,ZERO), CHI1) PSI1 = CHI1-ROOT - NAI = WAER(1) - 2.D0*ROOT + NAI = WAER(1) - TWO*ROOT SO4I = WAER(2) - ROOT ENDIF PSCONV1 = ABS(PSI1-PSIO1) .LE. EPS*PSIO1 @@ -2497,8 +2497,8 @@ C IF (NH4I*CLI .GT. A14) THEN BB =-(NH4I + CLI) CC =-A14 + NH4I*CLI - DD = BB*BB - 4.D0*CC - ROOT = 0.5D0*(-BB-SQRT(DD)) + DD = BB*BB - FOUR*CC + ROOT = HALF*(-BB-SQRT(DD)) IF (ROOT.GT.TINY) THEN ROOT = MIN(ROOT, CHI4) PSI4 = CHI4 - ROOT @@ -2513,18 +2513,18 @@ C C C SOLUTION ACIDIC OR BASIC? C - GG = 2.D0*SO4I + NO3I + CLI - NAI - NH4I + GG = TWO*SO4I + NO3I + CLI - NAI - NH4I IF (GG.GT.TINY) THEN ! H+ in excess BB =-GG CC =-AKW - DD = BB*BB - 4.D0*CC - HI = 0.5D0*(-BB + SQRT(DD)) + DD = BB*BB - FOUR*CC + HI = HALF*(-BB + SQRT(DD)) OHI= AKW/HI ELSE ! OH- in excess BB = GG CC =-AKW - DD = BB*BB - 4.D0*CC - OHI= 0.5D0*(-BB + SQRT(DD)) + DD = BB*BB - FOUR*CC + OHI= HALF*(-BB + SQRT(DD)) HI = AKW/OHI ENDIF C @@ -2534,7 +2534,7 @@ C CALL CALCAMAQ2 (-GG, NH4I, OHI, NH3AQ) HI = AKW/OHI ELSE - GGNO3 = MAX(2.D0*SO4I + NO3I - NAI - NH4I, ZERO) + GGNO3 = MAX(TWO*SO4I + NO3I - NAI - NH4I, ZERO) GGCL = MAX(GG-GGNO3, ZERO) IF (GGCL .GT.TINY) CALL CALCCLAQ2 (GGCL, CLI, HI, CLAQ) ! HCl IF (GGNO3.GT.TINY) THEN @@ -2762,9 +2762,9 @@ C IF (NH4I*CLI .GT. A14) THEN BB =-(WAER(3) + WAER(5) - ROOT3) CC =-A14 + NH4I*(WAER(5) - ROOT3) - DD = MAX(BB*BB - 4.D0*CC, ZERO) - ROOT2A= 0.5D0*(-BB+SQRT(DD)) - ROOT2B= 0.5D0*(-BB-SQRT(DD)) + DD = MAX(BB*BB - FOUR*CC, ZERO) + ROOT2A= HALF*(-BB+SQRT(DD)) + ROOT2B= HALF*(-BB-SQRT(DD)) IF (ZERO.LE.ROOT2A) THEN ROOT2 = ROOT2A ELSE @@ -2781,8 +2781,8 @@ C SODIUM SULFATE C IF (NAI*NAI*SO4I .GT. A5) THEN BB =-(CHI1 + WAER(1) - ROOT3) - CC = 0.25D0*(WAER(1) - ROOT3)*(4.D0*CHI1+WAER(1)-ROOT3) - DD =-0.25D0*(CHI1*(WAER(1)-ROOT3)**2.D0 - A5) + CC = 0.25D0*(WAER(1) - ROOT3)*(FOUR*CHI1+WAER(1)-ROOT3) + DD =-0.25D0*(CHI1*(WAER(1)-ROOT3)**TWO - A5) CALL POLY3(BB, CC, DD, ROOT1, ISLV) IF (ISLV.NE.0) ROOT1 = TINY ROOT1 = MIN (MAX(ROOT1,ZERO), MAX(WAER(1)-ROOT3,ZERO), @@ -2794,7 +2794,7 @@ C C C ION CONCENTRATIONS C - NAI = WAER(1) - (2.D0*ROOT1 + ROOT3) + NAI = WAER(1) - (TWO*ROOT1 + ROOT3) SO4I= WAER(2) - ROOT1 NH4I= WAER(3) - ROOT2 CLI = WAER(5) - (ROOT3 + ROOT2) @@ -2803,11 +2803,11 @@ C C SODIUM CHLORIDE ; To obtain new value for ROOT3 C IF (NAI*CLI .GT. A8) THEN - BB =-((CHI1-2.D0*ROOT1) + (WAER(5) - ROOT2)) - CC = (CHI1-2.D0*ROOT1)*(WAER(5) - ROOT2) - A8 - DD = SQRT(MAX(BB*BB - 4.D0*CC, TINY)) - ROOT3A= 0.5D0*(-BB-SQRT(DD)) - ROOT3B= 0.5D0*(-BB+SQRT(DD)) + BB =-((CHI1-TWO*ROOT1) + (WAER(5) - ROOT2)) + CC = (CHI1-TWO*ROOT1)*(WAER(5) - ROOT2) - A8 + DD = SQRT(MAX(BB*BB - FOUR*CC, TINY)) + ROOT3A= HALF*(-BB-SQRT(DD)) + ROOT3B= HALF*(-BB+SQRT(DD)) IF (ZERO.LE.ROOT3A) THEN ROOT3 = ROOT3A ELSE @@ -2821,18 +2821,18 @@ C C C SOLUTION ACIDIC OR BASIC? C - GG = 2.D0*SO4I + NO3I + CLI - NAI - NH4I + GG = TWO*SO4I + NO3I + CLI - NAI - NH4I IF (GG.GT.TINY) THEN ! H+ in excess BB =-GG CC =-AKW - DD = BB*BB - 4.D0*CC - HI = 0.5D0*(-BB + SQRT(DD)) + DD = BB*BB - FOUR*CC + HI = HALF*(-BB + SQRT(DD)) OHI= AKW/HI ELSE ! OH- in excess BB = GG CC =-AKW - DD = BB*BB - 4.D0*CC - OHI= 0.5D0*(-BB + SQRT(DD)) + DD = BB*BB - FOUR*CC + OHI= HALF*(-BB + SQRT(DD)) HI = AKW/OHI ENDIF C @@ -2842,7 +2842,7 @@ C CALL CALCAMAQ2 (-GG, NH4I, OHI, NH3AQ) HI = AKW/OHI ELSE - GGNO3 = MAX(2.D0*SO4I + NO3I - NAI - NH4I, ZERO) + GGNO3 = MAX(TWO*SO4I + NO3I - NAI - NH4I, ZERO) GGCL = MAX(GG-GGNO3, ZERO) IF (GGCL .GT.TINY) CALL CALCCLAQ2 (GGCL, CLI, HI, CLAQ) ! HCl IF (GGNO3.GT.TINY) THEN @@ -3116,10 +3116,10 @@ C IF (NH4I*CLI .GT. A14) THEN BB =-(WAER(3) + WAER(5) - ROOT3) CC = NH4I*(WAER(5) - ROOT3) - A14 - DD = MAX(BB*BB - 4.D0*CC, ZERO) + DD = MAX(BB*BB - FOUR*CC, ZERO) DD = SQRT(DD) - ROOT2A= 0.5D0*(-BB+DD) - ROOT2B= 0.5D0*(-BB-DD) + ROOT2A= HALF*(-BB+DD) + ROOT2B= HALF*(-BB-DD) IF (ZERO.LE.ROOT2A) THEN ROOT2 = ROOT2A ELSE @@ -3135,10 +3135,10 @@ C SODIUM SULFATE C IF (NAI*NAI*SO4I .GT. A5) THEN BB =-(WAER(2) + WAER(1) - ROOT3 - ROOT4) - CC = WAER(1)*(2.D0*ROOT3 + 2.D0*ROOT4 - 4.D0*WAER(2) - ONE) - & -(ROOT3 + ROOT4)**2.0 + 4.D0*WAER(2)*(ROOT3 + ROOT4) + CC = WAER(1)*(TWO*ROOT3 + TWO*ROOT4 - FOUR*WAER(2) - ONE) + & -(ROOT3 + ROOT4)**2.0 + FOUR*WAER(2)*(ROOT3 + ROOT4) CC =-0.25*CC - DD = WAER(1)*WAER(2)*(ONE - 2.D0*ROOT3 - 2.D0*ROOT4) + + DD = WAER(1)*WAER(2)*(ONE - TWO*ROOT3 - TWO*ROOT4) + & WAER(2)*(ROOT3 + ROOT4)**2.0 - A5 DD =-0.25*DD CALL POLY3(BB, CC, DD, ROOT1, ISLV) @@ -3152,11 +3152,11 @@ C C SODIUM NITRATE C IF (NAI*NO3I .GT. A9) THEN - BB =-(WAER(4) + WAER(1) - 2.D0*ROOT1 - ROOT3) - CC = WAER(4)*(WAER(1) - 2.D0*ROOT1 - ROOT3) - A9 - DD = SQRT(MAX(BB*BB - 4.D0*CC, TINY)) - ROOT4A= 0.5D0*(-BB-DD) - ROOT4B= 0.5D0*(-BB+DD) + BB =-(WAER(4) + WAER(1) - TWO*ROOT1 - ROOT3) + CC = WAER(4)*(WAER(1) - TWO*ROOT1 - ROOT3) - A9 + DD = SQRT(MAX(BB*BB - FOUR*CC, TINY)) + ROOT4A= HALF*(-BB-DD) + ROOT4B= HALF*(-BB+DD) IF (ZERO.LE.ROOT4A) THEN ROOT4 = ROOT4A ELSE @@ -3170,7 +3170,7 @@ C C C ION CONCENTRATIONS C - NAI = WAER(1) - (2.D0*ROOT1 + ROOT3 + ROOT4) + NAI = WAER(1) - (TWO*ROOT1 + ROOT3 + ROOT4) SO4I= WAER(2) - ROOT1 NH4I= WAER(3) - ROOT2 NO3I= WAER(4) - ROOT4 @@ -3179,11 +3179,11 @@ C C SODIUM CHLORIDE ; To obtain new value for ROOT3 C IF (NAI*CLI .GT. A8) THEN - BB =-(WAER(1) - 2.D0*ROOT1 + WAER(5) - ROOT2 - ROOT4) - CC = (WAER(5) + ROOT2)*(WAER(1) - 2.D0*ROOT1 - ROOT4) - A8 - DD = SQRT(MAX(BB*BB - 4.D0*CC, TINY)) - ROOT3A= 0.5D0*(-BB-DD) - ROOT3B= 0.5D0*(-BB+DD) + BB =-(WAER(1) - TWO*ROOT1 + WAER(5) - ROOT2 - ROOT4) + CC = (WAER(5) + ROOT2)*(WAER(1) - TWO*ROOT1 - ROOT4) - A8 + DD = SQRT(MAX(BB*BB - FOUR*CC, TINY)) + ROOT3A= HALF*(-BB-DD) + ROOT3B= HALF*(-BB+DD) IF (ZERO.LE.ROOT3A) THEN ROOT3 = ROOT3A ELSE @@ -3197,18 +3197,18 @@ C C C SOLUTION ACIDIC OR BASIC? C - GG = 2.D0*SO4I + NO3I + CLI - NAI - NH4I + GG = TWO*SO4I + NO3I + CLI - NAI - NH4I IF (GG.GT.TINY) THEN ! H+ in excess BB =-GG CC =-AKW - DD = BB*BB - 4.D0*CC - HI = 0.5D0*(-BB + SQRT(DD)) + DD = BB*BB - FOUR*CC + HI = HALF*(-BB + SQRT(DD)) OHI= AKW/HI ELSE ! OH- in excess BB = GG CC =-AKW - DD = BB*BB - 4.D0*CC - OHI= 0.5D0*(-BB + SQRT(DD)) + DD = BB*BB - FOUR*CC + OHI= HALF*(-BB + SQRT(DD)) HI = AKW/OHI ENDIF C @@ -3218,7 +3218,7 @@ C CALL CALCAMAQ2 (-GG, NH4I, OHI, NH3AQ) HI = AKW/OHI ELSE - GGNO3 = MAX(2.D0*SO4I + NO3I - NAI - NH4I, ZERO) + GGNO3 = MAX(TWO*SO4I + NO3I - NAI - NH4I, ZERO) GGCL = MAX(GG-GGNO3, ZERO) IF (GGCL .GT.TINY) CALL CALCCLAQ2 (GGCL, CLI, HI, CLAQ) ! HCl IF (GGNO3.GT.TINY) THEN diff --git a/src/MNH/mode_RBK90_Integrator.f90 b/src/MNH/mode_RBK90_Integrator.f90 index 4e1c88a42..8fdf8bc5e 100644 --- a/src/MNH/mode_RBK90_Integrator.f90 +++ b/src/MNH/mode_RBK90_Integrator.f90 @@ -733,7 +733,11 @@ Stage: DO istage = 1, ros_S END DO Err = SQRT(Err/N) +#if (MNH_REAL == 8) + ros_ErrorNorm = MAX(Err,1.0e-10) +#else ros_ErrorNorm = MAX(Err,1.0d-10) +#endif END FUNCTION ros_ErrorNorm diff --git a/src/include/isrpia.inc b/src/include/isrpia.inc index 8bf43f8bd..d60e7fbf4 100644 --- a/src/include/isrpia.inc +++ b/src/include/isrpia.inc @@ -87,7 +87,9 @@ C C *** GENERIC VARIABLES ************************************************ C CHARACTER VERSION*14 - COMMON /CGEN/ GREAT, TINY, TINY2, ZERO, ONE, VERSION + COMMON /CGEN/ GREAT, TINY, TINY2, + & ZERO, HALF, ONE, TWO, FOUR, + & VERSION C C *** END OF INCLUDE FILE ********************************************** C -- GitLab