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/UKMO/dev_r12866_HPC-02_Daley_Tiling_trial_extra_halo/src/TOP/PISCES/SED – NEMO

source: NEMO/branches/UKMO/dev_r12866_HPC-02_Daley_Tiling_trial_extra_halo/src/TOP/PISCES/SED/sedinorg.F90 @ 12958

Last change on this file since 12958 was 12958, checked in by hadcv, 4 years ago

Merge in trunk changes to r12933
File size: 10.8 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 ) 
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         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            zsolid1 = zvolc(ji,jk,jsopal) * solcp(ji,jk,jsopal)
138            zsatur = MAX(0., zundsat(ji,jk) / zsieq(ji) )
139            zsatur2 = (1.0 + temp(ji) / 400.0 )**37
140            znusil = ( 0.225 * ( 1.0 + temp(ji) / 15.) + 0.775 * zsatur2 * zsatur**2.25 ) / zsieq(ji)
141            zrearat1(ji,jk)  = ( reac_sil * znusil * dtsed * zsolid1 ) / &
142               &                ( 1. + reac_sil * znusil * dtsed * zundsat(ji,jk) )
143         ENDDO
144      ENDDO
145
146      CALL sed_mat( jwsil, jpoce, jpksed, zrearat1, zrearat2, zundsat, dtsed )
147
148      ! New solid concentration values (jk=2 to jksed) for each couple
149      DO jk = 2, jpksed
150         DO ji = 1, jpoce
151            zreasat = zrearat1(ji,jk) * zundsat(ji,jk) / ( zvolc(ji,jk,jsopal) )
152            solcp(ji,jk,jsopal) = solcp(ji,jk,jsopal) - zreasat
153         ENDDO
154      ENDDO
155
156      ! New pore water concentrations   
157      DO jk = 1, jpksed
158         DO ji = 1, jpoce
159            ! Acid Silicic
160            pwcp(ji,jk,jwsil)  = zsieq(ji) - zundsat(ji,jk)
161         ENDDO
162      ENDDO
163
164      !---------------------------------------------------------------
165      ! Performs CaCO3 particle deposition and redissolution (indice 9)
166      !--------------------------------------------------------------
167
168      ! computes co3por from the updated pwcp concentrations (note [co3por] = mol/l)
169
170      CALL sed_co3( kt )
171
172      ! *densSW(l)**2 converts aksps [mol2/kg sol2] into [mol2/l2] to get [undsat] in [mol/l]
173      DO jk = 1, jpksed
174         DO ji = 1, jpoce
175            zco3eq(ji)     = aksps(ji) * densSW(ji) * densSW(ji) / ( calcon2(ji) + rtrn )
176            zco3eq(ji)     = MAX( rtrn, zco3eq(ji) ) 
177            zundsat(ji,jk) = MAX(0., zco3eq(ji) - co3por(ji,jk) )
178         ENDDO
179      ENDDO
180
181      DO jk = 2, jpksed
182         DO ji = 1, jpoce
183            zsolid1 = zvolc(ji,jk,jscal) * solcp(ji,jk,jscal)
184            zrearat1(ji,jk) = ( reac_cal * dtsed * zsolid1 / zco3eq(ji) ) / &
185                  &               ( 1. + reac_cal * dtsed * zundsat(ji,jk) / zco3eq(ji) )
186         END DO
187      END DO
188
189      ! solves tridiagonal system
190      CALL sed_mat( jwdic, jpoce, jpksed, zrearat1, zrearat2, zundsat, dtsed )
191
192      ! New solid concentration values (jk=2 to jksed) for cacO3
193      DO jk = 2, jpksed
194         DO ji = 1, jpoce
195            zreasat = zrearat1(ji,jk) * zundsat(ji,jk) / zvolc(ji,jk,jscal)
196            solcp(ji,jk,jscal) = solcp(ji,jk,jscal) - zreasat
197         ENDDO
198      ENDDO
199
200      ! New dissolved concentrations
201      DO jk = 1, jpksed
202         DO ji = 1, jpoce
203            zreasat = zrearat1(ji,jk) * zundsat(ji,jk)   
204            ! For DIC
205            pwcp(ji,jk,jwdic)  = pwcp(ji,jk,jwdic) + zreasat
206            ! For alkalinity
207            pwcp(ji,jk,jwalk)  = pwcp(ji,jk,jwalk) + 2.0 * zreasat 
208         ENDDO
209      ENDDO
210
211      !-------------------------------------------------
212      ! Beginning DIC, Alkalinity
213      !-------------------------------------------------
214     
215      DO jk = 1, jpksed
216         DO ji = 1, jpoce     
217            zundsat(ji,jk)   = pwcp(ji,jk,jwdic)
218            zrearat1(ji,jk)  = 0.
219         ENDDO
220      ENDDO
221
222      ! solves tridiagonal system
223      CALL sed_mat( jwdic, jpoce, jpksed, zrearat1, zrearat2, zundsat, dtsed )
224
225     ! New dissolved concentrations     
226      DO jk = 1, jpksed
227         DO ji = 1, jpoce                     
228            pwcp(ji,jk,jwdic) = zundsat(ji,jk)
229         ENDDO
230      ENDDO           
231
232      !-------------------------------------------------
233      ! Beginning DIC, Alkalinity
234      !-------------------------------------------------
235
236      DO jk = 1, jpksed
237         DO ji = 1, jpoce
238            zundsat(ji,jk) = pwcp(ji,jk,jwalk)
239            zrearat1(ji,jk) = 0.
240         ENDDO
241      ENDDO
242!
243!      ! solves tridiagonal system
244      CALL sed_mat( jwalk, jpoce, jpksed, zrearat1, zrearat2, zundsat, dtsed )
245!
246!      ! New dissolved concentrations     
247      DO jk = 1, jpksed
248         DO ji = 1, jpoce
249            pwcp(ji,jk,jwalk) = zundsat(ji,jk)
250         ENDDO
251      ENDDO
252     
253      !----------------------------------
254      !   Back to initial geometry
255      !-----------------------------
256     
257      !---------------------------------------------------------------------
258      !   1/ Compensation for ajustement of the bottom water concentrations
259      !      (see note n° 1 about *por(2))
260      !--------------------------------------------------------------------
261      DO jw = 1, jpwat
262         DO ji = 1, jpoce
263            pwcp(ji,1,jw) = pwcp(ji,1,jw) + &
264               &            pwcp(ji,2,jw) * dzdep(ji) * por(2) / dzkbot(ji)
265         END DO
266      ENDDO
267     
268      !-----------------------------------------------------------------------
269      !    2/ Det of new rainrg taking account of the new weight fraction obtained
270      !      in dz3d(2) after diffusion/reaction (react/diffu are also in dzdep!)
271      !      This new rain (rgntg rm) will be used in advection/burial routine
272      !------------------------------------------------------------------------
273      DO js = 1, jpsol
274         DO ji = 1, jpoce
275            rainrg(ji,js) = raintg(ji) * solcp(ji,2,js)
276            rainrm(ji,js) = rainrg(ji,js) / mol_wgt(js)
277         END DO
278      ENDDO
279
280      !  New raintg
281      raintg(:) = 0.
282      DO js = 1, jpsol
283         DO ji = 1, jpoce
284            raintg(ji) = raintg(ji) + rainrg(ji,js)
285         END DO
286      ENDDO
287     
288      !--------------------------------
289      !    3/ back to initial geometry
290      !--------------------------------
291      DO ji = 1, jpoce
292         dz3d  (ji,2) = dz(2)
293         volw3d(ji,2) = dz3d(ji,2) * por(2)
294         vols3d(ji,2) = dz3d(ji,2) * por1(2)
295      ENDDO
296
297      IF( ln_timing )  CALL timing_stop('sed_inorg')
298!     
299   END SUBROUTINE sed_inorg
300
301END MODULE sedinorg
Note: See TracBrowser for help on using the repository browser.