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/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/src/TOP/PISCES/SED – NEMO

source: NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/src/TOP/PISCES/SED/sedinorg.F90 @ 10345

Last change on this file since 10345 was 10345, checked in by smasson, 5 years ago

dev_r10164_HPC09_ESIWACE_PREP_MERGE: merge with trunk@10344, see #2133

File size: 10.8 KB
RevLine 
[10225]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 ) 
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       ! 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
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 ) 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.    ;   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) * dz(jk)
89            zsolcpcl = zsolcpcl + solcp(ji,jk,jsclay) * dz(jk)
90         END DO
91         zsieq(ji) = sieqs(ji) * MAX(0.25, 1.0 - (0.045 * zsolcpcl / zsolcpsi )**0.58 )
92         zsieq(ji) = MAX( rtrn, sieqs(ji) )
93      END DO
94
95      DO js = 1, jpsol
96         DO jk = 1, jpksed
97            DO ji = 1, jpoce   
98               zvolc(ji,jk,js) = ( vols3d(ji,jk) * dens_mol_wgt(js) ) /  &
99                  &              ( volw3d(ji,jk) * 1.e-3  )     
100            ENDDO
101         ENDDO
102      ENDDO
103
104      !----------------------------------------------------------
105      ! 5.  Beginning of  Pore Water diffusion and solid reaction
106      !---------------------------------------------------------
107     
108      !-----------------------------------------------------------------------------
109      ! For jk=2,jpksed, and for couple
110      !  1 : jwsil/jsopal  ( SI/Opal )
111      !  2 : jsclay/jsclay ( clay/clay )
112      !  3 : jwoxy/jspoc   ( O2/POC )
113      !  reaction rate is a function of solid=concentration in solid reactif in [mol/l]
114      !  and undersaturation in [mol/l].
115      !  Solid weight fractions should be in ie [mol/l])
116      !  second member and solution are in zundsat variable
117      !-------------------------------------------------------------------------
118
119      DO jk = 1, jpksed
120         DO ji = 1, jpoce
121            ! For Silicic Acid and clay
122            zundsat(ji,jk) = zsieq(ji) - pwcp(ji,jk,jwsil)
123         ENDDO
124      ENDDO
125     
126      ! Definition of reaction rates [rearat]=sans dim
127      ! For jk=1 no reaction (pure water without solid) for each solid compo
128      DO ji = 1, jpoce
129         zrearat1(ji,:) = 0.
130         zrearat2(ji,:) = 0.
131      ENDDO
132
133      ! left hand side of coefficient matrix
134      DO jk = 2, jpksed
135         DO ji = 1, jpoce
136            zsolid1 = zvolc(ji,jk,jsopal) * solcp(ji,jk,jsopal)
137            zsatur = MAX(0., zundsat(ji,jk) / zsieq(ji) )
138            zsatur2 = (1.0 + temp(ji) / 400.0 )**37
139            znusil = ( 0.225 * ( 1.0 + temp(ji) / 15.) + 0.775 * zsatur2 * zsatur**2.25 ) / zsieq(ji)
140            zrearat1(ji,jk)  = ( reac_sil * znusil * dtsed * zsolid1 ) / &
141               &                ( 1. + reac_sil * znusil * dtsed * zundsat(ji,jk) )
142         ENDDO
143      ENDDO
144
145      CALL sed_mat( jwsil, jpoce, jpksed, zrearat1, zrearat2, zundsat, dtsed )
146
147      ! New solid concentration values (jk=2 to jksed) for each couple
148      DO jk = 2, jpksed
149         DO ji = 1, jpoce
150            zreasat = zrearat1(ji,jk) * zundsat(ji,jk) / ( zvolc(ji,jk,jsopal) )
151            solcp(ji,jk,jsopal) = solcp(ji,jk,jsopal) - zreasat
152         ENDDO
153      ENDDO
154
155      ! New pore water concentrations   
156      DO jk = 1, jpksed
157         DO ji = 1, jpoce
158            ! Acid Silicic
159            pwcp(ji,jk,jwsil)  = zsieq(ji) - zundsat(ji,jk)
160         ENDDO
161      ENDDO
162
163      !---------------------------------------------------------------
164      ! Performs CaCO3 particle deposition and redissolution (indice 9)
165      !--------------------------------------------------------------
166
167      ! computes co3por from the updated pwcp concentrations (note [co3por] = mol/l)
168
169      CALL sed_co3( kt )
170
171      ! *densSW(l)**2 converts aksps [mol2/kg sol2] into [mol2/l2] to get [undsat] in [mol/l]
172      DO jk = 1, jpksed
173         DO ji = 1, jpoce
174            zco3eq(ji)     = aksps(ji) * densSW(ji) * densSW(ji) / ( calcon2(ji) + rtrn )
175            zco3eq(ji)     = MAX( rtrn, zco3eq(ji) ) 
176            zundsat(ji,jk) = MAX(0., zco3eq(ji) - co3por(ji,jk) )
177         ENDDO
178      ENDDO
179
180      DO jk = 2, jpksed
181         DO ji = 1, jpoce
182            zsolid1 = zvolc(ji,jk,jscal) * solcp(ji,jk,jscal)
183            zrearat1(ji,jk) = ( reac_cal * dtsed * zsolid1 / zco3eq(ji) ) / &
184                  &               ( 1. + reac_cal * dtsed * zundsat(ji,jk) / zco3eq(ji) )
185         END DO
186      END DO
187
188      ! solves tridiagonal system
189      CALL sed_mat( jwdic, jpoce, jpksed, zrearat1, zrearat2, zundsat, dtsed )
190
191      ! New solid concentration values (jk=2 to jksed) for cacO3
192      DO jk = 2, jpksed
193         DO ji = 1, jpoce
194            zreasat = zrearat1(ji,jk) * zundsat(ji,jk) / zvolc(ji,jk,jscal)
195            solcp(ji,jk,jscal) = solcp(ji,jk,jscal) - zreasat
196         ENDDO
197      ENDDO
198
199      ! New dissolved concentrations
200      DO jk = 1, jpksed
201         DO ji = 1, jpoce
202            zreasat = zrearat1(ji,jk) * zundsat(ji,jk)   
203            ! For DIC
204            pwcp(ji,jk,jwdic)  = pwcp(ji,jk,jwdic) + zreasat
205            ! For alkalinity
206            pwcp(ji,jk,jwalk)  = pwcp(ji,jk,jwalk) + 2.0 * zreasat 
207         ENDDO
208      ENDDO
209
210      !-------------------------------------------------
211      ! Beginning DIC, Alkalinity
212      !-------------------------------------------------
213     
214      DO jk = 1, jpksed
215         DO ji = 1, jpoce     
216            zundsat(ji,jk)   = pwcp(ji,jk,jwdic)
217            zrearat1(ji,jk)  = 0.
218         ENDDO
219      ENDDO
220
221      ! solves tridiagonal system
222      CALL sed_mat( jwdic, jpoce, jpksed, zrearat1, zrearat2, zundsat, dtsed )
223
224     ! New dissolved concentrations     
225      DO jk = 1, jpksed
226         DO ji = 1, jpoce                     
227            pwcp(ji,jk,jwdic) = zundsat(ji,jk)
228         ENDDO
229      ENDDO           
230
231      !-------------------------------------------------
232      ! Beginning DIC, Alkalinity
233      !-------------------------------------------------
234
235      DO jk = 1, jpksed
236         DO ji = 1, jpoce
237            zundsat(ji,jk) = pwcp(ji,jk,jwalk)
238            zrearat1(ji,jk) = 0.
239         ENDDO
240      ENDDO
241!
242!      ! solves tridiagonal system
243      CALL sed_mat( jwalk, jpoce, jpksed, zrearat1, zrearat2, zundsat, dtsed )
244!
245!      ! New dissolved concentrations     
246      DO jk = 1, jpksed
247         DO ji = 1, jpoce
248            pwcp(ji,jk,jwalk) = zundsat(ji,jk)
249         ENDDO
250      ENDDO
251     
252      !----------------------------------
253      !   Back to initial geometry
254      !-----------------------------
255     
256      !---------------------------------------------------------------------
257      !   1/ Compensation for ajustement of the bottom water concentrations
258      !      (see note n° 1 about *por(2))
259      !--------------------------------------------------------------------
260      DO jw = 1, jpwat
261         DO ji = 1, jpoce
262            pwcp(ji,1,jw) = pwcp(ji,1,jw) + &
263               &            pwcp(ji,2,jw) * dzdep(ji) * por(2) / dzkbot(ji)
264         END DO
265      ENDDO
266     
267      !-----------------------------------------------------------------------
268      !    2/ Det of new rainrg taking account of the new weight fraction obtained
269      !      in dz3d(2) after diffusion/reaction (react/diffu are also in dzdep!)
270      !      This new rain (rgntg rm) will be used in advection/burial routine
271      !------------------------------------------------------------------------
272      DO js = 1, jpsol
273         DO ji = 1, jpoce
274            rainrg(ji,js) = raintg(ji) * solcp(ji,2,js)
275            rainrm(ji,js) = rainrg(ji,js) / mol_wgt(js)
276         END DO
277      ENDDO
278
279      !  New raintg
280      raintg(:) = 0.
281      DO js = 1, jpsol
282         DO ji = 1, jpoce
283            raintg(ji) = raintg(ji) + rainrg(ji,js)
284         END DO
285      ENDDO
286     
287      !--------------------------------
288      !    3/ back to initial geometry
289      !--------------------------------
290      DO ji = 1, jpoce
291         dz3d  (ji,2) = dz(2)
292         volw3d(ji,2) = dz3d(ji,2) * por(2)
293         vols3d(ji,2) = dz3d(ji,2) * por1(2)
294      ENDDO
295
296      IF( ln_timing )  CALL timing_stop('sed_inorg')
297!     
298   END SUBROUTINE sed_inorg
299
300END MODULE sedinorg
Note: See TracBrowser for help on using the repository browser.