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

Last change on this file since 14323 was 14323, checked in by aumont, 5 months ago

major update of the sediment model

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         ENDDO
166      ENDDO
167
168      DO jk = 2, jpksed
169         DO ji = 1, jpoce
170            IF ( zundsat(ji,jk) > 0. ) THEN
171               zsolid1 = zvolc(ji,jk,jscal) * solcp(ji,jk,jscal)
172               zbeta = 1.0 - reac_cal * dtsed2 * zundsat(ji,jk) + reac_cal * dtsed2 * zsolid1
173               zgamma = -reac_cal * dtsed2 * zundsat(ji,jk)
174               zundsat(ji,jk) = ( -zbeta + SQRT( zbeta**2 - 4.0 * zgamma ) ) / ( 2.0 * reac_sil * znusil * dtsed2 )
175               saturco3(ji,jk) = zundsat(ji,jk)
176               zreasat = reac_cal * dtsed2 * zundsat(ji,jk) * zsolid1 / ( 1. + reac_cal * dtsed2 * zundsat(ji,jk) )
177               solcp(ji,jk,jscal) = solcp(ji,jk,jscal) - zreasat / zvolc(ji,jk,jscal)
178               ! For DIC
179               pwcp(ji,jk,jwdic)  = pwcp(ji,jk,jwdic) + zreasat
180               ! For alkalinity
181               pwcp(ji,jk,jwalk)  = pwcp(ji,jk,jwalk) + 2.0 * zreasat
182            ENDIF
183         END DO
184      END DO
185
186      IF( ln_timing )  CALL timing_stop('sed_inorg')
187!     
188   END SUBROUTINE sed_inorg
189
190END MODULE sedinorg
Note: See TracBrowser for help on using the repository browser.