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.
sedsol.F90 in NEMO/branches/2021/dev_r14383_PISCES_NEWDEV_PISCO/src/TOP/PISCES/SED – NEMO

source: NEMO/branches/2021/dev_r14383_PISCES_NEWDEV_PISCO/src/TOP/PISCES/SED/sedsol.F90 @ 15351

Last change on this file since 15351 was 15351, checked in by cetlod, 3 years ago

dev_PISCO : minor bugfix

File size: 10.2 KB
Line 
1MODULE sedsol
2   !!======================================================================
3   !!              ***  MODULE  sedsol  ***
4   !!    Sediment : dissolution and reaction in pore water related
5   !!    related to organic matter
6   !!    Diffusion of solutes in pore water
7   !!=====================================================================
8   !! * Modules used
9   USE sed     ! sediment global variable
10   USE sed_oce
11   USE sedini
12   USE sedfunc
13   USE seddsr
14   USE sedjac
15   USE sedbtb
16   USE sedco3
17   USE sedmat
18   USE sedorg
19# if ! defined key_agrif
20   USE trosk
21#endif
22   USE lib_mpp         ! distribued memory computing library
23   USE lib_fortran
24
25   IMPLICIT NONE
26   PRIVATE
27
28   PUBLIC sed_sol
29
30   !! $Id: sedsol.F90 5215 2015-04-15 16:11:56Z nicolasmartin $
31CONTAINS
32   
33   SUBROUTINE sed_sol( kt ) 
34      !!----------------------------------------------------------------------
35      !!                   ***  ROUTINE sed_sol  ***
36      !!
37      !!  ** Purpose :  computes pore water diffusion and reactions
38      !!
39      !!  ** Methode :  Computation of the redox and dissolution reactions
40      !!                in the sediment.
41      !!                The main redox reactions are solved in sed_dsr whereas
42      !!                the secondary reactions are solved in sed_dsr_redoxb.
43      !!                Inorganic dissolution is solved in sed_inorg
44      !!                A strand spliting approach is being used here (see
45      !!                sed_dsr_redoxb for more information).
46      !!                Diffusive fluxes are computed in sed_diff
47      !!
48      !!   History :
49      !!        !  98-08 (E. Maier-Reimer, Christoph Heinze )  Original code
50      !!        !  04-10 (N. Emprin, M. Gehlen ) f90
51      !!        !  06-04 (C. Ethe)  Re-organization
52      !!        !  19-08 (O. Aumont) Debugging and improvement of the model.
53      !!                             The original method is replaced by a
54      !!                             Strand splitting method which deals
55      !!                             well with stiff reactions.
56      !!----------------------------------------------------------------------
57      !! Arguments
58      INTEGER, INTENT(in) ::   kt
59      ! --- local variables
60      INTEGER  :: ji, jk, js, jw, jn, neq   ! dummy looop indices
61      REAL(wp), DIMENSION( jpoce, jpvode * jpksed ) :: ZXIN, FVAL
62# if ! defined key_agrif     
63      REAL(wp), DIMENSION(jpoce,jpksed) :: preac
64      INTEGER :: JINDEX, ITOL, IJAC, MLJAC, IMAX
65      INTEGER :: MUJAC,LE1, LJAC, LWORK
66      INTEGER :: IDID, NMAXSTP, ROSM
67      REAL(wp) :: X, XEND
68      REAL(wp),DIMENSION(jpoce) :: H
69      INTEGER, DIMENSION(jpoce) :: accmask
70      REAL(wp), DIMENSION(jpvode * jpksed) :: RTOL, ATOL
71      REAL(wp), POINTER :: WORK(:)
72      INTEGER, DIMENSION(jpoce,3)   :: ISTAT
73      REAL(wp), DIMENSION(jpoce,2)  :: RSTAT
74#endif     
75      !!
76      !!----------------------------------------------------------------------
77
78      IF( ln_timing )  CALL timing_start('sed_sol')
79!
80      IF( kt == nitsed000 ) THEN
81         IF (lwp) THEN
82            WRITE(numsed,*) ' sed_sol : Organic/inorganic degradation related reactions and diffusion'
83            WRITE(numsed,*) ' '
84         ENDIF
85!         !
86      ENDIF
87
88      ! 1. Change of geometry
89      !    Increase of dz3d(2) thickness : dz3d(2) = dz3d(2)+dzdep
90      !    Warning : no change for dz(2)
91      !---------------------------------------------------------
92      dz3d(1:jpoce,2) = dz3d(1:jpoce,2) + dzdep(1:jpoce)
93
94      ! New values for volw3d(:,2) and vols3d(:,2)
95      ! Warning : no change neither for volw(2) nor  vols(2)
96      !------------------------------------------------------
97      volw3d(1:jpoce,2) = dz3d(1:jpoce,2) * por(2)
98      vols3d(1:jpoce,2) = dz3d(1:jpoce,2) * por1(2)
99
100      ! 2. Change of previous solid fractions (due to volum changes) for k=2
101      !---------------------------------------------------------------------
102      DO js = 1, jpsol
103         solcp(:,2,js) = solcp(:,2,js) * dz(2) / dz3d(:,2)
104      END DO
105
106      ! 3. New solid fractions (including solid rain fractions) for k=2
107      !------------------------------------------------------------------
108      DO js = 1, jpsol
109         DO ji = 1, jpoce
110            IF (dzdep(ji) .ne. 0) THEN
111               solcp(ji,2,js) = solcp(ji,2,js) + rainrg(ji,js) * dtsed / ( por1(2) * dz3d(ji,2) )
112               ! rainrm are temporary cancel
113            ENDIF
114         END DO
115      ENDDO
116
117      ! 4.  Adjustment of bottom water concen.(pwcp(1)):
118      !     We impose that pwcp(2) is constant. Including dzdep in dz3d(:,2) we assume
119      !     that dzdep has got a porosity of por(2). So pore water volum of jk=2 increase.
120      !     To keep pwcp(2) cste we must compensate this "increase" by a slight adjusment
121      !     of bottom water concentration.
122      !     This adjustment is compensate at the end of routine
123      !-------------------------------------------------------------
124      DO jw = 1, jpwat
125         pwcp(:,1,jw) = pwcp(:,1,jw) - pwcp(:,2,jw) * dzdep(:) * por(2) / dzkbot(:)
126      ENDDO
127
128      ! --------------------------------------------------
129      ! Computation of the diffusivities
130      ! --------------------------------------------------
131      DO js = 1, jpwat
132         DO jk = 1, jpksed
133            diff(:,jk,js) = ( diff1(js) + diff2(js) * temp(:) ) / ( 1.0 - 2.0 * log( por(jk) ) )
134         END DO
135      END DO
136
137      ! Impact of bioirrigation and adsorption on diffusion
138      ! ---------------------------------------------------
139      diff(:,:,jwsil) = diff(:,:,jwsil) * ( 1.0 + irrig(:,:) )
140      diff(:,:,jwoxy) = diff(:,:,jwoxy) * ( 1.0 + irrig(:,:) )
141      diff(:,:,jwdic) = diff(:,:,jwdic) * ( 1.0 + irrig(:,:) )
142      diff(:,:,jwno3) = diff(:,:,jwno3) * ( 1.0 + irrig(:,:) )
143      diff(:,:,jwpo4) = diff(:,:,jwpo4) * ( 1.0 + irrig(:,:) )
144      diff(:,:,jwalk) = diff(:,:,jwalk) * ( 1.0 + irrig(:,:) )
145      diff(:,:,jwh2s) = diff(:,:,jwh2s) * ( 1.0 + irrig(:,:) )
146      diff(:,:,jwso4) = diff(:,:,jwso4) * ( 1.0 + irrig(:,:) )
147      diff(:,:,jwlgw) = diff(:,:,jwlgw) * ( 1.0 + irrig(:,:) )
148      diff(:,:,jwnh4) = diff(:,:,jwnh4) * ( 1.0 + irrig(:,:) )
149      diff(:,:,jwfe2) = diff(:,:,jwfe2) * ( 1.0 + 0.1 * irrig(:,:) )
150
151      ! Conversion of volume units
152      !----------------------------
153      DO js = 1, jpsol
154         DO jk = 1, jpksed
155            volc(:,jk,js) = ( vols3d(:,jk) / mol_wgt(js) ) /  &
156            &                 ( volw3d(:,jk) * 1.e-3 )
157         ENDDO
158      ENDDO
159
160      ! Compute coefficients commonly used in diffusion
161      CALL sed_mat_coef
162      ! Apply bioturbation and compute the impact of the slow SMS on species
163      CALL sed_btb( kt )
164      ! Recompute pH after bioturbation and slow chemistry
165      CALL sed_co3( kt )
166
167# if ! defined key_agrif
168      ! The following part deals with the stiff ODEs
169      ! This is the expensive part of the code and should be carefully
170      ! chosen. We use the DVODE solver after many trials to find a cheap
171      ! way to solve the ODEs. This is not necessarily the most efficient
172      ! but this is the one that was not too much of a pain to code and that
173      ! was the most precise and quick.
174      ! The ones I tried : operator splitting (Strang), hybrid spectral methods
175      ! Brent, Powell's hybrid method, ...
176      ! ---------------------------------------------------------------------
177      NEQ  = jpvode * jpksed
178      XEND = dtsed 
179      RTOL = rosrtol
180      ATOL = rosatol
181      ITOL = 1
182      IJAC = 1
183      DO jn = 1, NEQ
184         js = jarr(jn,2)
185         IF (js == jwfe2) ATOL(jn) = rosatol / 100.0
186      END DO
187      MLJAC = jpvode
188      MUJAC = jpvode
189      LE1  = 2*MLJAC+MUJAC+1
190      LJAC = MLJAC+MUJAC+1
191      LWORK = NEQ*(LJAC+LE1+8) + 5
192      ALLOCATE(WORK(LWORK) )
193
194      X    = 0.0
195      H(:)    = dtsed
196      WORK  = 0.
197
198      ! Put all the species in one local array (nb of tracers * vertical
199      ! dimension
200      DO jn = 1, NEQ
201         jk = jarr(jn,1)
202         js = jarr(jn,2)
203         IF (js <= jpwat) THEN
204            zxin(:,jn) = pwcp(:,jk,js) * 1E6
205         ELSE
206            zxin(:,jn) = solcp(:,jk,js-jpwat) * 1E6
207         ENDIF
208      END DO
209
210      ! Set options for VODE : banded matrix. SParse option is much more
211      ! expensive except if one computes the sparse Jacobian explicitly
212      ! To speed up the computation, one way is to reduce ATOL and RTOL
213      ! which may be a good option at the beginning of the simulations
214      ! during the spin up
215      ! ----------------------------------------------------------------
216      CALL ROSK(NROSORDER, NEQ,X,zxin,XEND,H,RTOL,ATOL,ITOL, sed_jac,  &
217           &   MLJAC,MUJAC,WORK,LWORK,IDID,ISTAT,RSTAT)
218
219      accmask(:) = 0.0
220      CALL sed_func( NEQ, ZXIN, FVAL, ACCMASK )
221
222      DO jn = 1, NEQ
223         jk = jarr(jn,1)
224         js = jarr(jn,2)
225         IF (js <= jpwat) THEN
226            pwcp(:,jk,js) = zxin(:,jn) * 1E-6
227         ELSE
228            solcp(:,jk,js-jpwat) = zxin(:,jn) * 1E-6
229         ENDIF
230      END DO
231      rstepros(:) = ISTAT(:,3)
232
233      DEALLOCATE( WORK )
234#endif
235
236      preac(:,:) = 0.
237      CALL sed_mat_dsri( jwpo4, preac, pwcpa(:,:,jwpo4), dtsed, pwcp(:,:,jwpo4) )
238      CALL sed_mat_dsri( jwalk, preac, pwcpa(:,:,jwalk), dtsed, pwcp(:,:,jwalk) )
239
240      CALL sed_org( kt )
241
242      !----------------------------------
243      !   Back to initial geometry
244      !-----------------------------
245
246      !---------------------------------------------------------------------
247      !   1/ Compensation for ajustement of the bottom water concentrations
248      !      (see note n� 1 about *por(2))
249      !--------------------------------------------------------------------
250      DO jw = 1, jpwat
251         pwcp(:,1,jw) = pwcp(:,1,jw) + pwcp(:,2,jw) * dzdep(:) * por(2) / dzkbot(:)
252      ENDDO
253
254      ! 2. Change of previous solid fractions (due to volum changes) for k=2
255      !---------------------------------------------------------------------
256      DO js = 1, jpsol
257         solcp(:,2,js) = solcp(:,2,js) * dz3d(:,2) / dz(2)
258      END DO
259
260      !--------------------------------
261      !    3/ back to initial geometry
262      !--------------------------------
263      dz3d  (:,2) = dz(2)
264      volw3d(:,2) = dz3d(:,2) * por(2)
265      vols3d(:,2) = dz3d(:,2) * por1(2)
266
267      IF( ln_timing )  CALL timing_stop('sed_sol')
268!     
269   END SUBROUTINE sed_sol
270
271   SUBROUTINE JACFAC
272 
273   END SUBROUTINE JACFAC
274
275END MODULE sedsol
Note: See TracBrowser for help on using the repository browser.