[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

<none>



Friends,
I may be stupid but i didn't understand all of Bob Keys programs.
why so complicated
The easiest way to compute the chemistry is to insert CO3-- in the
expression of alk (A=... below). then alk is a function of hydrogen
ions only which has to be determined from given alk and Sco2.
for positive hydrogen (highly desireable) the function is strictly
monotonous ( negative first derivative) and the second derivative
is positive everywhere.
The following example (obviously vectorizeable) works in course of time integration.
For initialisation i recommend to set Hi deliberately too small,
say 10**-9. Due to the positive curvature the procedure converges
within no more than 8 iterations without overshooting.
If anybody insists to include also PO4 and Si.. (far below any precision
of measurements) it has to be treated like BT.
Cheers
Ernst
       DO 43 KI=1,2

       DO 43 K=1,KE
       DO 43 J=1,JE
       DO 43 I=1,IE

        H=HI(I,J,K)
        R=CO3(I,J,K)
        C=SCO212(I,J,K)
      T1=h/ak13(i,j,k)
      T2=h/ak23(i,j,k)
      AKW=AKW3(I,J,K)
      BT=BORAT(I,J,K)
      AKB=AKB3(I,J,K)
      alk=alkali(i,j,k)
      A=C*(2.+t2)/(1.+t2+t2*t1)  +AKW/H-H+BT/(1.+H/AKB)-ALK

      DADH=C*( 1./(AK2*(1.+T2+T2*T1))- (2.+T2)*(1./AK2+2.*T1/AK2)/
     1 (1.+T2+T2*T1)**2)
     1 -AKW/H**2-1.-(BT/AKB)/(1.+H/AKB)**2
      dddhhh=a/dadh
    

      H=H-dddhhh
  
      HI(I,J,K)=HI(I,J,K)-dddhhh
      co3(i,j,k)=c/(1.+hi(i,j,k)*(1.+hi(i,j,k)/ak13(i,j,k))/ak23(i,j,k))
43      CONTINUE