New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
sedinorg.F90 in NEMO/branches/2019/dev_r11708_aumont_PISCES_QUOTA/src/TOP/PISCES/SED – NEMO

source: NEMO/branches/2019/dev_r11708_aumont_PISCES_QUOTA/src/TOP/PISCES/SED/sedinorg.F90 @ 14325

Last change on this file since 14325 was 14325, checked in by aumont, 3 years ago

bug in the CacO3 dissolution scheme
File size: 7.7 KB
Line 
1MODULE sedinorg
2   !!======================================================================
3   !!              ***  MODULE  sedinorg  ***
4   !!    Sediment : dissolution and reaction in pore water of
5   !!               inorganic species
6   !!=====================================================================
7   !! * Modules used
8   USE sed     ! sediment global variable
9   USE sed_oce
10   USE sedmat  ! linear system of equations
11   USE sedco3  ! carbonate ion and proton concentration
12   USE sedini
13   USE seddsr
14   USE lib_mpp         ! distribued memory computing library
15   USE lib_fortran
16
17   IMPLICIT NONE
18   PRIVATE
19
20   PUBLIC sed_inorg
21
22   !! $Id: seddsr.F90 5215 2015-04-15 16:11:56Z nicolasmartin $
23CONTAINS
24   
25   SUBROUTINE sed_inorg( kt, knt ) 
26      !!----------------------------------------------------------------------
27      !!                   ***  ROUTINE sed_inorg  ***
28      !!
29      !!  ** Purpose :  computes pore water dissolution and reaction
30      !!
31      !!  ** Methode :  implicit simultaneous computation of undersaturation
32      !!               resulting from diffusive pore water transport and chemical
33      !!               pore water reactions. Solid material is consumed according
34      !!               to redissolution and remineralisation
35      !!
36      !!  ** Remarks :
37      !!              - undersaturation : deviation from saturation concentration
38      !!              - reaction rate   : sink of undersaturation from dissolution
39      !!                                 of solid material
40      !!
41      !!   History :
42      !!        !  98-08 (E. Maier-Reimer, Christoph Heinze )  Original code
43      !!        !  04-10 (N. Emprin, M. Gehlen ) f90
44      !!        !  06-04 (C. Ethe)  Re-organization
45      !!        !  19-08 (O. Aumont) Debugging and improvement of the model
46      !!----------------------------------------------------------------------
47      !! Arguments
48      INTEGER, INTENT(in) ::   kt, knt       ! number of iteration
49      ! --- local variables
50      INTEGER :: ji, jk, js, jw         ! dummy looop indices
51      REAL(wp), DIMENSION(jpoce,jpksed) :: zrearat1, zrearat2    ! reaction rate in pore water
52      REAL(wp), DIMENSION(jpoce,jpksed) :: zundsat    ! undersaturation ; indice jpwatp1 is for calcite   
53      REAL(wp), DIMENSION(jpoce) :: zco3eq
54      REAL(wp), DIMENSION(jpoce,jpksed,jpsol) :: zvolc    ! temp. variables
55      REAL(wp), DIMENSION(jpoce) :: zsieq
56      REAL(wp)  ::  zsolid1, zvolw, zreasat, zbeta, zgamma
57      REAL(wp)  ::  zsatur, zsatur2, znusil, zsolcpcl, zsolcpsi
58      !!
59      !!----------------------------------------------------------------------
60
61      IF( ln_timing )  CALL timing_start('sed_inorg')
62!
63      IF( kt == nitsed000 .AND. knt == 1 ) THEN
64         IF (lwp) THEN
65            WRITE(numsed,*) ' sed_inorg : Dissolution reaction '
66            WRITE(numsed,*) ' '
67         ENDIF
68!         !
69      ENDIF
70
71     ! Initializations
72     !----------------------
73     
74      zrearat1(:,:) = 0.    ;   zundsat(:,:)  = 0. 
75      zrearat2(:,:) = 0.
76      zco3eq(:)     = rtrn
77      zvolc(:,:,:)  = 0.
78
79      ! -----------------------------------------------
80      ! Computation of Si solubility
81      ! Param of Ridgwell et al. 2002
82      ! -----------------------------------------------
83
84      DO ji = 1, jpoce
85         zsolcpcl = 0.0
86         zsolcpsi = 0.0
87         DO jk = 1, jpksed
88            zsolcpsi = zsolcpsi + solcp(ji,jk,jsopal) * vols3d(ji,jk)
89            zsolcpcl = zsolcpcl + solcp(ji,jk,jsclay) * vols3d(ji,jk)
90         END DO
91         zsolcpsi = MAX( zsolcpsi, rtrn )
92         zsieq(ji) = sieqs(ji) * MAX(0.25, 1.0 - (0.045 * zsolcpcl / zsolcpsi )**0.58 )
93         zsieq(ji) = MAX( rtrn, sieqs(ji) )
94      END DO
95
96      DO js = 1, jpsol
97         DO jk = 1, jpksed
98            DO ji = 1, jpoce   
99               zvolc(ji,jk,js) = ( vols3d(ji,jk) * dens_mol_wgt(js) ) /  &
100                  &              ( volw3d(ji,jk) * 1.e-3  )     
101            ENDDO
102         ENDDO
103      ENDDO
104
105      !----------------------------------------------------------
106      ! 5.  Beginning of  Pore Water diffusion and solid reaction
107      !---------------------------------------------------------
108     
109      !-----------------------------------------------------------------------------
110      ! For jk=2,jpksed, and for couple
111      !  1 : jwsil/jsopal  ( SI/Opal )
112      !  2 : jsclay/jsclay ( clay/clay )
113      !  3 : jwoxy/jspoc   ( O2/POC )
114      !  reaction rate is a function of solid=concentration in solid reactif in [mol/l]
115      !  and undersaturation in [mol/l].
116      !  Solid weight fractions should be in ie [mol/l])
117      !  second member and solution are in zundsat variable
118      !-------------------------------------------------------------------------
119
120      DO jk = 1, jpksed
121         DO ji = 1, jpoce
122            ! For Silicic Acid and clay
123            zundsat(ji,jk) = zsieq(ji) - pwcp(ji,jk,jwsil)
124         ENDDO
125      ENDDO
126     
127      ! Definition of reaction rates [rearat]=sans dim
128      ! For jk=1 no reaction (pure water without solid) for each solid compo
129      DO ji = 1, jpoce
130         zrearat1(ji,:) = 0.
131         zrearat2(ji,:) = 0.
132      ENDDO
133
134      ! left hand side of coefficient matrix
135      DO jk = 2, jpksed
136         DO ji = 1, jpoce
137            IF ( zundsat(ji,jk) > 0. ) THEN
138               zsolid1 = zvolc(ji,jk,jsopal) * solcp(ji,jk,jsopal)
139               zsatur = zundsat(ji,jk) / zsieq(ji)
140               zsatur2 = (1.0 + temp(ji) / 400.0 )**37
141               znusil = ( 0.225 * ( 1.0 + temp(ji) / 15.) + 0.775 * zsatur2 * zsatur**8.25 ) / zsieq(ji)
142               zbeta = 1.0 - reac_sil * znusil * dtsed2 * zundsat(ji,jk) + reac_sil * znusil * dtsed2 * zsolid1
143               zgamma = - reac_sil * znusil * dtsed2 * zundsat(ji,jk)
144               zundsat(ji,jk) = ( - zbeta + SQRT( zbeta**2 - 4.0 * zgamma ) ) / ( 2.0 * reac_sil * znusil * dtsed2 )
145               solcp(ji,jk,jsopal ) = zsolid1 / ( 1. + reac_sil * znusil * dtsed2 * zundsat(ji,jk) ) / zvolc(ji,jk,jsopal)
146               pwcp(ji,jk,jwsil) = zsieq(ji) - zundsat(ji,jk) 
147            ENDIF
148         ENDDO
149      ENDDO
150
151      !---------------------------------------------------------------
152      ! Performs CaCO3 particle deposition and redissolution (indice 9)
153      !--------------------------------------------------------------
154
155      ! computes co3por from the updated pwcp concentrations (note [co3por] = mol/l)
156
157      CALL sed_co3( kt )
158
159      ! *densSW(l)**2 converts aksps [mol2/kg sol2] into [mol2/l2] to get [undsat] in [mol/l]
160      DO jk = 1, jpksed
161         DO ji = 1, jpoce
162            zco3eq(ji)     = aksps(ji) * densSW(ji) * densSW(ji) / ( calcon2(ji) + rtrn )
163            zco3eq(ji)     = MAX( rtrn, zco3eq(ji) ) 
164            zundsat(ji,jk) = ( zco3eq(ji) - co3por(ji,jk) ) / zco3eq(ji)
165            saturco3(ji,jk) = zundsat(ji,jk)
166         ENDDO
167      ENDDO
168
169      DO jk = 2, jpksed
170         DO ji = 1, jpoce
171            IF ( zundsat(ji,jk) > 0. ) THEN
172               zsolid1 = zvolc(ji,jk,jscal) * solcp(ji,jk,jscal)
173               zbeta = 1.0 - reac_cal * dtsed2 * zundsat(ji,jk) + reac_cal * dtsed2 * zsolid1
174               zgamma = -reac_cal * dtsed2 * zundsat(ji,jk)
175               zundsat(ji,jk) = ( -zbeta + SQRT( zbeta**2 - 4.0 * zgamma ) ) / ( 2.0 * reac_cal * dtsed2 )
176               saturco3(ji,jk) = zundsat(ji,jk)
177               zreasat = reac_cal * dtsed2 * zundsat(ji,jk) * zsolid1 / ( 1. + reac_cal * dtsed2 * zundsat(ji,jk) )
178               solcp(ji,jk,jscal) = solcp(ji,jk,jscal) - zreasat / zvolc(ji,jk,jscal)
179               ! For DIC
180               pwcp(ji,jk,jwdic)  = pwcp(ji,jk,jwdic) + zreasat
181               ! For alkalinity
182               pwcp(ji,jk,jwalk)  = pwcp(ji,jk,jwalk) + 2.0 * zreasat
183            ENDIF
184         END DO
185      END DO
186
187      IF( ln_timing )  CALL timing_stop('sed_inorg')
188!     
189   END SUBROUTINE sed_inorg
190
191END MODULE sedinorg
Note: See TracBrowser for help on using the repository browser.