Index: NEMO/branches/2019/dev_r11708_aumont_PISCES_QUOTA/src/TOP/PISCES/SED/sedinorg.F90
===================================================================
--- NEMO/branches/2019/dev_r11708_aumont_PISCES_QUOTA/src/TOP/PISCES/SED/sedinorg.F90	(revision 12360)
+++ NEMO/branches/2019/dev_r11708_aumont_PISCES_QUOTA/src/TOP/PISCES/SED/sedinorg.F90	(revision 14323)
@@ -23,5 +23,5 @@
 CONTAINS
    
-   SUBROUTINE sed_inorg( kt ) 
+   SUBROUTINE sed_inorg( kt, knt ) 
       !!----------------------------------------------------------------------
       !!                   ***  ROUTINE sed_inorg  ***
@@ -46,5 +46,5 @@
       !!----------------------------------------------------------------------
       !! Arguments
-      INTEGER, INTENT(in) ::   kt       ! number of iteration
+      INTEGER, INTENT(in) ::   kt, knt       ! number of iteration
       ! --- local variables
       INTEGER :: ji, jk, js, jw         ! dummy looop indices
@@ -54,5 +54,5 @@
       REAL(wp), DIMENSION(jpoce,jpksed,jpsol) :: zvolc    ! temp. variables
       REAL(wp), DIMENSION(jpoce) :: zsieq
-      REAL(wp)  ::  zsolid1, zvolw, zreasat
+      REAL(wp)  ::  zsolid1, zvolw, zreasat, zbeta, zgamma
       REAL(wp)  ::  zsatur, zsatur2, znusil, zsolcpcl, zsolcpsi
       !!
@@ -61,5 +61,5 @@
       IF( ln_timing )  CALL timing_start('sed_inorg')
 !
-      IF( kt == nitsed000 ) THEN
+      IF( kt == nitsed000 .AND. knt == 1 ) THEN
          IF (lwp) THEN
             WRITE(numsed,*) ' sed_inorg : Dissolution reaction '
@@ -73,5 +73,5 @@
       
       zrearat1(:,:) = 0.    ;   zundsat(:,:)  = 0. 
-      zrearat2(:,:) = 0.    ;   zrearat2(:,:) = 0.
+      zrearat2(:,:) = 0.
       zco3eq(:)     = rtrn
       zvolc(:,:,:)  = 0.
@@ -86,6 +86,6 @@
          zsolcpsi = 0.0
          DO jk = 1, jpksed
-            zsolcpsi = zsolcpsi + solcp(ji,jk,jsopal) * dz(jk)
-            zsolcpcl = zsolcpcl + solcp(ji,jk,jsclay) * dz(jk)
+            zsolcpsi = zsolcpsi + solcp(ji,jk,jsopal) * vols3d(ji,jk)
+            zsolcpcl = zsolcpcl + solcp(ji,jk,jsclay) * vols3d(ji,jk)
          END DO
          zsolcpsi = MAX( zsolcpsi, rtrn )
@@ -135,28 +135,15 @@
       DO jk = 2, jpksed
          DO ji = 1, jpoce
-            zsolid1 = zvolc(ji,jk,jsopal) * solcp(ji,jk,jsopal)
-            zsatur = MAX(0., zundsat(ji,jk) / zsieq(ji) )
-            zsatur2 = (1.0 + temp(ji) / 400.0 )**37
-            znusil = ( 0.225 * ( 1.0 + temp(ji) / 15.) + 0.775 * zsatur2 * zsatur**2.25 ) / zsieq(ji)
-            zrearat1(ji,jk)  = ( reac_sil * znusil * dtsed * zsolid1 ) / &
-               &                ( 1. + reac_sil * znusil * dtsed * zundsat(ji,jk) )
-         ENDDO
-      ENDDO
-
-      CALL sed_mat( jwsil, jpoce, jpksed, zrearat1, zrearat2, zundsat, dtsed )
-
-      ! New solid concentration values (jk=2 to jksed) for each couple 
-      DO jk = 2, jpksed
-         DO ji = 1, jpoce
-            zreasat = zrearat1(ji,jk) * zundsat(ji,jk) / ( zvolc(ji,jk,jsopal) )
-            solcp(ji,jk,jsopal) = solcp(ji,jk,jsopal) - zreasat
-         ENDDO
-      ENDDO
-
-      ! New pore water concentrations    
-      DO jk = 1, jpksed
-         DO ji = 1, jpoce
-            ! Acid Silicic 
-            pwcp(ji,jk,jwsil)  = zsieq(ji) - zundsat(ji,jk)
+            IF ( zundsat(ji,jk) > 0. ) THEN
+               zsolid1 = zvolc(ji,jk,jsopal) * solcp(ji,jk,jsopal)
+               zsatur = zundsat(ji,jk) / zsieq(ji)
+               zsatur2 = (1.0 + temp(ji) / 400.0 )**37
+               znusil = ( 0.225 * ( 1.0 + temp(ji) / 15.) + 0.775 * zsatur2 * zsatur**8.25 ) / zsieq(ji)
+               zbeta = 1.0 - reac_sil * znusil * dtsed2 * zundsat(ji,jk) + reac_sil * znusil * dtsed2 * zsolid1
+               zgamma = - reac_sil * znusil * dtsed2 * zundsat(ji,jk)
+               zundsat(ji,jk) = ( - zbeta + SQRT( zbeta**2 - 4.0 * zgamma ) ) / ( 2.0 * reac_sil * znusil * dtsed2 )
+               solcp(ji,jk,jsopal ) = zsolid1 / ( 1. + reac_sil * znusil * dtsed2 * zundsat(ji,jk) ) / zvolc(ji,jk,jsopal)
+               pwcp(ji,jk,jwsil) = zsieq(ji) - zundsat(ji,jk) 
+            ENDIF
          ENDDO
       ENDDO
@@ -175,5 +162,5 @@
             zco3eq(ji)     = aksps(ji) * densSW(ji) * densSW(ji) / ( calcon2(ji) + rtrn )
             zco3eq(ji)     = MAX( rtrn, zco3eq(ji) ) 
-            zundsat(ji,jk) = MAX(0., zco3eq(ji) - co3por(ji,jk) )
+            zundsat(ji,jk) = ( zco3eq(ji) - co3por(ji,jk) ) / zco3eq(ji)
          ENDDO
       ENDDO
@@ -181,117 +168,19 @@
       DO jk = 2, jpksed
          DO ji = 1, jpoce
-            zsolid1 = zvolc(ji,jk,jscal) * solcp(ji,jk,jscal)
-            zrearat1(ji,jk) = ( reac_cal * dtsed * zsolid1 / zco3eq(ji) ) / &
-                  &               ( 1. + reac_cal * dtsed * zundsat(ji,jk) / zco3eq(ji) )
+            IF ( zundsat(ji,jk) > 0. ) THEN
+               zsolid1 = zvolc(ji,jk,jscal) * solcp(ji,jk,jscal)
+               zbeta = 1.0 - reac_cal * dtsed2 * zundsat(ji,jk) + reac_cal * dtsed2 * zsolid1
+               zgamma = -reac_cal * dtsed2 * zundsat(ji,jk)
+               zundsat(ji,jk) = ( -zbeta + SQRT( zbeta**2 - 4.0 * zgamma ) ) / ( 2.0 * reac_sil * znusil * dtsed2 )
+               saturco3(ji,jk) = zundsat(ji,jk)
+               zreasat = reac_cal * dtsed2 * zundsat(ji,jk) * zsolid1 / ( 1. + reac_cal * dtsed2 * zundsat(ji,jk) )
+               solcp(ji,jk,jscal) = solcp(ji,jk,jscal) - zreasat / zvolc(ji,jk,jscal)
+               ! For DIC
+               pwcp(ji,jk,jwdic)  = pwcp(ji,jk,jwdic) + zreasat
+               ! For alkalinity
+               pwcp(ji,jk,jwalk)  = pwcp(ji,jk,jwalk) + 2.0 * zreasat
+            ENDIF
          END DO
       END DO
-
-      ! solves tridiagonal system
-      CALL sed_mat( jwdic, jpoce, jpksed, zrearat1, zrearat2, zundsat, dtsed )
-
-      ! New solid concentration values (jk=2 to jksed) for cacO3
-      DO jk = 2, jpksed
-         DO ji = 1, jpoce
-            zreasat = zrearat1(ji,jk) * zundsat(ji,jk) / zvolc(ji,jk,jscal)
-            solcp(ji,jk,jscal) = solcp(ji,jk,jscal) - zreasat
-         ENDDO
-      ENDDO
-
-      ! New dissolved concentrations
-      DO jk = 1, jpksed
-         DO ji = 1, jpoce
-            zreasat = zrearat1(ji,jk) * zundsat(ji,jk)    
-            ! For DIC
-            pwcp(ji,jk,jwdic)  = pwcp(ji,jk,jwdic) + zreasat
-            ! For alkalinity
-            pwcp(ji,jk,jwalk)  = pwcp(ji,jk,jwalk) + 2.0 * zreasat 
-         ENDDO
-      ENDDO
-
-      !-------------------------------------------------
-      ! Beginning DIC, Alkalinity 
-      !-------------------------------------------------
-      
-      DO jk = 1, jpksed
-         DO ji = 1, jpoce      
-            zundsat(ji,jk)   = pwcp(ji,jk,jwdic)
-            zrearat1(ji,jk)  = 0.
-         ENDDO
-      ENDDO
-
-      ! solves tridiagonal system
-      CALL sed_mat( jwdic, jpoce, jpksed, zrearat1, zrearat2, zundsat, dtsed )
-
-     ! New dissolved concentrations      
-      DO jk = 1, jpksed
-         DO ji = 1, jpoce                      
-            pwcp(ji,jk,jwdic) = zundsat(ji,jk)
-         ENDDO
-      ENDDO            
-
-      !-------------------------------------------------
-      ! Beginning DIC, Alkalinity 
-      !-------------------------------------------------
-
-      DO jk = 1, jpksed
-         DO ji = 1, jpoce
-            zundsat(ji,jk) = pwcp(ji,jk,jwalk)
-            zrearat1(ji,jk) = 0.
-         ENDDO
-      ENDDO
-!
-!      ! solves tridiagonal system
-      CALL sed_mat( jwalk, jpoce, jpksed, zrearat1, zrearat2, zundsat, dtsed )
-!
-!      ! New dissolved concentrations      
-      DO jk = 1, jpksed
-         DO ji = 1, jpoce
-            pwcp(ji,jk,jwalk) = zundsat(ji,jk)
-         ENDDO
-      ENDDO
-      
-      !----------------------------------
-      !   Back to initial geometry
-      !-----------------------------
-      
-      !---------------------------------------------------------------------
-      !   1/ Compensation for ajustement of the bottom water concentrations
-      !      (see note n° 1 about *por(2))
-      !--------------------------------------------------------------------
-      DO jw = 1, jpwat
-         DO ji = 1, jpoce
-            pwcp(ji,1,jw) = pwcp(ji,1,jw) + &
-               &            pwcp(ji,2,jw) * dzdep(ji) * por(2) / dzkbot(ji)
-         END DO
-      ENDDO
-      
-      !-----------------------------------------------------------------------
-      !    2/ Det of new rainrg taking account of the new weight fraction obtained 
-      !      in dz3d(2) after diffusion/reaction (react/diffu are also in dzdep!)
-      !      This new rain (rgntg rm) will be used in advection/burial routine
-      !------------------------------------------------------------------------
-      DO js = 1, jpsol
-         DO ji = 1, jpoce
-            rainrg(ji,js) = raintg(ji) * solcp(ji,2,js)
-            rainrm(ji,js) = rainrg(ji,js) / mol_wgt(js)
-         END DO
-      ENDDO
-
-      !  New raintg
-      raintg(:) = 0.
-      DO js = 1, jpsol
-         DO ji = 1, jpoce
-            raintg(ji) = raintg(ji) + rainrg(ji,js)
-         END DO
-      ENDDO
-      
-      !--------------------------------
-      !    3/ back to initial geometry
-      !--------------------------------
-      DO ji = 1, jpoce
-         dz3d  (ji,2) = dz(2)
-         volw3d(ji,2) = dz3d(ji,2) * por(2)
-         vols3d(ji,2) = dz3d(ji,2) * por1(2)
-      ENDDO
 
       IF( ln_timing )  CALL timing_stop('sed_inorg')
