Changeset 15297 for NEMO/branches/2021/dev_r14383_PISCES_NEWDEV_PISCO/src/TOP/PISCES/SED/sedinorg.F90

Timestamp:

2021-09-28T11:20:56+02:00 (3 years ago)

Author:

aumont

Message:

major update of the sediment module

File:

• NEMO/branches/2021/dev_r14383_PISCES_NEWDEV_PISCO/src/TOP/PISCES/SED/sedinorg.F90 (modified) (4 diffs)

NEMO/branches/2021/dev_r14383_PISCES_NEWDEV_PISCO/src/TOP/PISCES/SED/sedinorg.F90

@@ -7 +7 @@
    !! * Modules used
    USE sed     ! sediment global variable
+   USE sed_oce
    USE sedini
+   USE sedmat
    USE lib_mpp         ! distribued memory computing library
    USE lib_fortran
@@ -45 +47 @@
       ! --- local variables
       INTEGER   ::  ji,jk          ! dummy looop indices
-      REAL(wp)  ::  zsieq
-      REAL(wp)  ::  zsolid1, zreasat
+      REAL(wp), DIMENSION(jpoce) ::  zsieq, reac_silf
+      REAL(wp)  ::  zsolid1, zreasat, zco3sat
       REAL(wp)  ::  zsatur, zsatur2, znusil, zsolcpcl, zsolcpsi, zexcess
+      REAL(wp), DIMENSION(jpoce,jpksed) :: zundsat, zrearat, psms
       !!
       !!----------------------------------------------------------------------
@@ -71 +74 @@
          END DO
          zsolcpsi = MAX( zsolcpsi, rtrn )
-         zsieq = sieqs(ji) * MAX(0.25, 1.0 - (0.045 * zsolcpcl / zsolcpsi )**0.58 )
+         zsieq(ji) = sieqs(ji) * MAX(0.2, 1.0 - (0.045 * zsolcpcl / zsolcpsi )**0.58 )
+         reac_silf(ji) = reac_sil * ( 0.05 + 0.055 * ( 1.64 * ( zsolcpcl / zsolcpsi + 0.01 ) )**(-0.75) ) / 1.25 
+      END DO
 
-         !----------------------------------------------------------
-         ! 5.  Beginning of  Pore Water diffusion and solid reaction
-         !---------------------------------------------------------
       
-         !-----------------------------------------------------------------------------
-         ! For jk=2,jpksed, and for couple 
-         !  1 : jwsil/jsopal  ( SI/Opal )
-         !  2 : jsclay/jsclay ( clay/clay ) 
-         !  3 : jwoxy/jspoc   ( O2/POC )
-         !  reaction rate is a function of solid=concentration in solid reactif in [mol/l] 
-         !  and undersaturation in [mol/l].
-         !  Solid weight fractions should be in ie [mol/l])
-         !  second member and solution are in zundsat variable
-         !-------------------------------------------------------------------------
+      DO ji = 1, jpoce
          DO jk = 2, jpksed
             zsolid1 = volc(ji,jk,jsopal) * solcp(ji,jk,jsopal)
-            zsatur = MAX(0., ( zsieq - pwcp(ji,jk,jwsil) ) / zsieq )
+            zsatur = MAX(0., ( zsieq(ji) - pwcp(ji,jk,jwsil) ) / zsieq(ji) )
             zsatur2 = (1.0 + temp(ji) / 400.0 )**37
             znusil = ( 0.225 * ( 1.0 + temp(ji) / 15.) * zsatur + 0.775 * zsatur2 * zsatur**9.25 )
-            solcp(ji,jk,jsopal) = solcp(ji,jk,jsopal) - reac_sil * znusil * dtsed * solcp(ji,jk,jsopal)
-            pwcp(ji,jk,jwsil) = pwcp(ji,jk,jwsil) + reac_sil * znusil * dtsed * zsolid1
+            solcp(ji,jk,jsopal) = solcp(ji,jk,jsopal) - reac_silf(ji) * znusil * dtsed * solcp(ji,jk,jsopal)
+            pwcp(ji,jk,jwsil) = pwcp(ji,jk,jwsil) + reac_silf(ji) * znusil * dtsed * zsolid1
          END DO
       END DO
@@ -104 +97 @@
       ! *densSW(l)**2 converts aksps [mol2/kg sol2] into [mol2/l2] to get [undsat] in [mol/l]
       DO ji = 1, jpoce
-         saturco3(ji,:) = 1.0 - co3por(ji,:) * calcon2(ji) / ( aksps(ji) * densSW(ji) * densSW(ji) + rtrn ) 
-         DO jk = 2, jpksed
+         zco3sat = aksps(ji) * densSW(ji) * densSW(ji) / ( calcon2(ji) + rtrn )
+         saturco3(ji,:) = 1.0 - co3por(ji,:) / ( rtrn + zco3sat )
+         DO jk = 1, jpksed
             zsolid1 = volc(ji,jk,jscal) * solcp(ji,jk,jscal)
-            zexcess = MAX( 0., saturco3(ji,jk) ) 
-            zreasat = reac_cal * dtsed * zexcess * zsolid1
-            solcp(ji,jk,jscal) = solcp(ji,jk,jscal) - zreasat / volc(ji,jk,jscal)
+            zundsat(ji,jk) = MAX( 0., zco3sat - co3por(ji,jk) )
+            zrearat(ji,jk) = ( reac_cal * zsolid1 / ( zco3sat + rtrn ) ) / &
+            &                ( 1. + reac_cal * dtsed * zundsat(ji,jk) / ( zco3sat + rtrn ) )
+         END DO
+      END DO
+
+      psms(:,:) = 0.0
+      ! solves tridiagonal system
+      CALL sed_mat_dsri( jwdic, -zrearat(:,:), psms(:,:), dtsed, zundsat )
+
+      ! New solid concentration values (jk=2 to jksed) for cacO3
+      DO jk = 2, jpksed
+         DO ji = 1, jpoce
+            zreasat = zrearat(ji,jk) * dtsed * zundsat(ji,jk) / ( volc(ji,jk,jscal) + rtrn )
+            solcp(ji,jk,jscal) = solcp(ji,jk,jscal) - zreasat
+            zreasat = zrearat(ji,jk) * dtsed * 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
-         END DO
-      END DO
+            pwcp(ji,jk,jwalk)  = pwcp(ji,jk,jwalk) + 2.* zreasat
+         ENDDO
+      ENDDO
+
 
       IF( ln_timing )  CALL timing_stop('sed_inorg')