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.
Changeset 10288 for NEMO/branches/2018/dev_r9866_HPC_03_globcom/src/TOP/PISCES/SED/seddsr.F90 – NEMO

Ignore:
Timestamp:
2018-11-07T18:25:49+01:00 (5 years ago)
Author:
francesca
Message:

reduce global communications, see #2010

File:
1 edited

Legend:

Unmodified
Added
Removed

  • NEMO/branches/2018/dev_r9866_HPC_03_globcom/src/TOP/PISCES/SED/seddsr.F90

    r5215 r10288  
    11MODULE seddsr 
    2 #if defined key_sed 
    32   !!====================================================================== 
    43   !!              ***  MODULE  seddsr  *** 
    5    !!    Sediment : dissolution and reaction in pore water 
     4   !!    Sediment : dissolution and reaction in pore water related  
     5   !!    related to organic matter 
    66   !!===================================================================== 
    77   !! * Modules used 
    88   USE sed     ! sediment global variable 
    9    USE sedmat  ! linear system of equations 
    10    USE sedco3  ! carbonate ion and proton concentration  
     9   USE sed_oce 
     10   USE sedini 
     11   USE lib_mpp         ! distribued memory computing library 
     12   USE lib_fortran 
     13 
     14   IMPLICIT NONE 
     15   PRIVATE 
    1116 
    1217   PUBLIC sed_dsr 
    … …  
    1419   !! * Module variables 
    1520 
    16    REAL(wp), DIMENSION(:), ALLOCATABLE, PUBLIC :: cons_o2 
    17    REAL(wp), DIMENSION(:), ALLOCATABLE, PUBLIC :: cons_no3 
    18    REAL(wp), DIMENSION(:), ALLOCATABLE, PUBLIC :: sour_no3 
    19    REAL(wp), DIMENSION(:), ALLOCATABLE, PUBLIC :: sour_c13 
    20    REAL(wp), DIMENSION(:), ALLOCATABLE, PUBLIC ::  dens_mol_wgt  ! molecular density  
     21   REAL(wp) :: zadsnh4 
     22   REAL(wp), DIMENSION(jpsol), PUBLIC      :: dens_mol_wgt  ! molecular density  
     23   REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: zvolc    ! temp. variables 
     24 
    2125 
    2226   !! $Id$ 
    2327CONTAINS 
    2428    
    25    SUBROUTINE sed_dsr( kt )  
     29   SUBROUTINE sed_dsr( kt, knt )  
    2630      !!---------------------------------------------------------------------- 
    2731      !!                   ***  ROUTINE sed_dsr  *** 
    … …  
    2933      !!  ** Purpose :  computes pore water dissolution and reaction 
    3034      !! 
    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  
     35      !!  ** Methode :  Computation of the redox reactions in sediment. 
     36      !!                The main redox reactions are solved in sed_dsr whereas 
     37      !!                the secondary reactions are solved in sed_dsr_redoxb. 
     38      !!                A strand spliting approach is being used here (see  
     39      !!                sed_dsr_redoxb for more information).  
    4040      !! 
    4141      !!   History : 
    … …  
    4343      !!        !  04-10 (N. Emprin, M. Gehlen ) f90 
    4444      !!        !  06-04 (C. Ethe)  Re-organization 
     45      !!        !  19-08 (O. Aumont) Debugging and improvement of the model. 
     46      !!                             The original method is replaced by a  
     47      !!                              Strand splitting method which deals  
     48      !!                              well with stiff reactions. 
    4549      !!---------------------------------------------------------------------- 
    4650      !! Arguments 
    47       INTEGER, INTENT(in) ::   kt       ! number of iteration 
     51      INTEGER, INTENT(in) ::   kt, knt       ! number of iteration 
    4852      ! --- local variables 
    49       INTEGER :: ji, jk, js, jw   ! dummy looop indices 
    50       INTEGER :: nv               ! number of variables in linear tridiagonal eq 
    51  
    52       REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: zrearat    ! reaction rate in pore water 
    53       REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: zundsat    ! undersaturation ; indice jpwatp1 is for calcite    
    54       REAL(wp), DIMENSION(:    ), ALLOCATABLE :: zmo2_0, zmo2_1  ! temp. array for mass balance calculation 
    55       REAL(wp), DIMENSION(:    ), ALLOCATABLE :: zmno3_0, zmno3_1, zmno3_2 
    56       REAL(wp), DIMENSION(:    ), ALLOCATABLE :: zmc13_0, zmc13_1, zmc13_2, zmc13_3 
    57       REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: zvolc    ! temp. variables 
     53      INTEGER :: ji, jk, js, jw, jn   ! dummy looop indices 
     54 
     55      REAL(wp), DIMENSION(jpoce,jpksed) :: zrearat1, zrearat2, zrearat3    ! reaction rate in pore water 
     56      REAL(wp), DIMENSION(jpoce,jpksed) :: zundsat    ! undersaturation ; indice jpwatp1 is for calcite    
     57      REAL(wp), DIMENSION(jpoce,jpksed) :: zkpoc, zkpos, zkpor, zlimo2, zlimno3, zlimso4, zlimfeo    ! undersaturation ; indice jpwatp1 is for calcite    
     58      REAL(wp), DIMENSION(jpoce)        :: zsumtot 
    5859      REAL(wp)  ::  zsolid1, zsolid2, zsolid3, zvolw, zreasat 
    59  
     60      REAL(wp)  ::  zsatur, zsatur2, znusil, zkpoca, zkpocb, zkpocc 
     61      REAL(wp)  ::  zratio, zgamma, zbeta, zlimtmp, zundsat2 
    6062      !! 
    6163      !!---------------------------------------------------------------------- 
    6264 
    63       IF( kt == nitsed000 ) THEN 
    64          WRITE(numsed,*) ' sed_dsr : Dissolution reaction ' 
    65          WRITE(numsed,*) ' ' 
    66          !  
    67          ALLOCATE( dens_mol_wgt((jpoce) )  
    68          dens_mol_wgt(1:jpsol) = dens / mol_wgt(1:jpsol) 
    69          !  
    70          ALLOCATE( cons_o2 (jpoce) ) ;   ALLOCATE( cons_no3(jpoce) )  
    71          ALLOCATE( sour_no3(jpoce) ) ;   ALLOCATE( sour_c13(jpoce) )   
     65      IF( ln_timing )  CALL timing_start('sed_dsr') 
     66! 
     67      IF( kt == nitsed000 .AND. knt == 1 ) THEN 
     68         IF (lwp) THEN 
     69            WRITE(numsed,*) ' sed_dsr : Dissolution reaction ' 
     70            WRITE(numsed,*) ' ' 
     71         ENDIF 
    7272      ENDIF 
    7373 
    74       ! Initialization of data for mass balance calculation 
    75       !--------------------------------------------------- 
    76  
    77       tokbot(:,:)  = 0. 
    78       cons_o2 (:)  = 0.  
    79       cons_no3(:)  = 0.   
    80       sour_no3(:)  = 0.   
    81       sour_c13(:)  = 0.        
    82     
    83       ! Initializations 
    84       !---------------------- 
    85       ALLOCATE( zmo2_0 (jpoce) ) ;  ALLOCATE( zmo2_1 (jpoce) )  
    86       ALLOCATE( zmno3_0(jpoce) ) ;  ALLOCATE( zmno3_1(jpoce) )  ;  ALLOCATE( zmno3_2(jpoce) )  
    87       ALLOCATE( zmc13_0(jpoce) ) ;  ALLOCATE( zmc13_1(jpoce) )  ;  ALLOCATE( zmc13_2(jpoce) ) ;  ALLOCATE( zmc13_3(jpoce) )  
    88  
    89       zmo2_0 (:) = 0.  ; zmo2_1 (:) = 0. 
    90       zmno3_0(:) = 0.  ; zmno3_1(:) = 0.  ;  zmno3_2(:) = 0. 
    91       zmc13_0(:) = 0.  ; zmc13_1(:) = 0.  ;  zmc13_2(:) = 0.  ; zmc13_3(:) = 0. 
     74     ! Initializations 
     75     !---------------------- 
    9276       
    93       ALLOCATE( zrearat(jpoce,jpksed,3) ) ;  ALLOCATE( zundsat(jpoce,jpksed,3) )       
    94       zrearat(:,:,:)   = 0.    ;   zundsat(:,:,:) = 0.  
    95  
    96  
    97       ALLOCATE( zvolc(jpoce,jpksed,jpsol) )  
    98       zvolc(:,:,:)   = 0. 
    99  
    100       !-------------------------------------------------------------------- 
    101       ! Temporary accomodation to take account of  particule rain deposition 
    102       !--------------------------------------------------------------------- 
    103        
    104        
    105       ! 1. Change of geometry 
    106       !    Increase of dz3d(2) thickness : dz3d(2) = dz3d(2)+dzdep 
    107       !    Warning : no change for dz(2) 
    108       !--------------------------------------------------------- 
    109       dz3d(1:jpoce,2) = dz3d(1:jpoce,2) + dzdep(1:jpoce) 
    110  
    111        
    112       ! New values for volw3d(:,2) and vols3d(:,2) 
    113       ! Warning : no change neither for volw(2) nor  vols(2) 
    114       !------------------------------------------------------ 
    115       volw3d(1:jpoce,2) = dz3d(1:jpoce,2) * por(2) 
    116       vols3d(1:jpoce,2) = dz3d(1:jpoce,2) * por1(2) 
     77      zrearat1(:,:)   = 0.    ;   zundsat(:,:) = 0. ; zkpoc(:,:) = 0. 
     78      zlimo2 (:,:)    = 0.    ;   zlimno3(:,:) = 0. ; zrearat2(:,:) = 0. 
     79      zlimso4(:,:)    = 0.    ;   zkpor(:,:)   = 0. ; zrearat3(:,:) = 0. 
     80      zkpos  (:,:)    = 0. 
     81      zsumtot(:)      = rtrn 
     82   
     83      ALLOCATE( zvolc(jpoce, jpksed, jpsol) ) 
     84      zvolc(:,:,:)    = 0. 
     85      zadsnh4 = 1.0 / ( 1.0 + adsnh4 ) 
     86 
     87      ! Inhibition terms for the different redox equations 
     88      ! -------------------------------------------------- 
     89      DO jk = 1, jpksed 
     90         DO ji = 1, jpoce 
     91            zkpoc(ji,jk) = reac_pocl  
     92            zkpos(ji,jk) = reac_pocs 
     93            zkpor(ji,jk) = reac_pocr 
     94         END DO 
     95      END DO 
    11796 
    11897      ! Conversion of volume units 
    … …  
    12099      DO js = 1, jpsol 
    121100         DO jk = 1, jpksed 
    122             DO ji = 1, jpoce     
     101            DO ji = 1, jpoce 
    123102               zvolc(ji,jk,js) = ( vols3d(ji,jk) * dens_mol_wgt(js) ) /  & 
    124                   &              ( volw3d(ji,jk) * 1.e-3 )      
     103                  &              ( volw3d(ji,jk) * 1.e-3 ) 
    125104            ENDDO 
    126105         ENDDO 
    127106      ENDDO 
    128107 
    129       ! 2. Change of previous solid fractions (due to volum changes) for k=2 
    130       !--------------------------------------------------------------------- 
    131  
    132       DO js = 1, jpsol 
    133          DO ji = 1, jpoce 
    134             solcp(ji,2,js) = solcp(ji,2,js) * dz(2) / dz3d(ji,2) 
    135          ENDDO 
    136       END DO 
    137  
    138       ! 3. New solid fractions (including solid rain fractions) for k=2   
    139       !------------------------------------------------------------------    
    140       DO js = 1, jpsol 
    141          DO ji = 1, jpoce 
    142             solcp(ji,2,js) = solcp(ji,2,js) + & 
    143             &           ( rainrg(ji,js) / raintg(ji) ) * ( dzdep(ji) / dz3d(ji,2) ) 
    144             ! rainrm are temporary cancel 
    145             rainrm(ji,js) = 0. 
    146          END DO 
    147       ENDDO 
    148  
    149       ! 4.  Adjustment of bottom water concen.(pwcp(1)):  
    150       !     We impose that pwcp(2) is constant. Including dzdep in dz3d(:,2) we assume  
    151       !     that dzdep has got a porosity of por(2). So pore water volum of jk=2 increase. 
    152       !     To keep pwcp(2) cste we must compensate this "increase" by a slight adjusment 
    153       !     of bottom water concentration. 
    154       !     This adjustment is compensate at the end of routine 
    155       !------------------------------------------------------------- 
    156       DO jw = 1, jpwat 
    157          DO ji = 1, jpoce 
    158             pwcp(ji,1,jw) = pwcp(ji,1,jw) - & 
    159                &            pwcp(ji,2,jw) * dzdep(ji) * por(2) / dzkbot(ji) 
    160          END DO 
    161       ENDDO 
    162  
    163   
    164108      !---------------------------------------------------------- 
    165       ! 5.  Beginning of  Pore Water diffusion and solid reaction 
     109      ! 5.  Beginning of solid reaction 
    166110      !--------------------------------------------------------- 
    167        
    168       !----------------------------------------------------------------------------- 
    169       ! For jk=2,jpksed, and for couple  
    170       !  1 : jwsil/jsopal  ( SI/Opal ) 
    171       !  2 : jsclay/jsclay ( clay/clay )  
    172       !  3 : jwoxy/jspoc   ( O2/POC ) 
    173       !  reaction rate is a function of solid=concentration in solid reactif in [mol/l]  
    174       !  and undersaturation in [mol/l]. 
    175       !  Solid weight fractions should be in ie [mol/l]) 
    176       !  second member and solution are in zundsat variable 
    177       !------------------------------------------------------------------------- 
    178  
    179       !number of variables 
    180       nv  = 3 
    181       
    182       DO jk = 1, jpksed 
    183          DO ji = 1, jpoce 
    184             ! For Silicic Acid and clay 
    185             zundsat(ji,jk,1) = sat_sil   - pwcp(ji,jk,jwsil) 
    186             zundsat(ji,jk,2) = sat_clay 
    187             ! For O2 
    188             zundsat(ji,jk,3) = pwcp(ji,jk,jwoxy) / so2ut  
    189          ENDDO 
    190       ENDDO 
    191        
    192111       
    193112      ! Definition of reaction rates [rearat]=sans dim  
    194113      ! For jk=1 no reaction (pure water without solid) for each solid compo 
    195       DO ji = 1, jpoce 
    196          zrearat(ji,1,:) = 0. 
    197       ENDDO 
    198  
     114      zrearat1(:,:) = 0. 
     115      zrearat2(:,:) = 0. 
     116      zrearat3(:,:) = 0. 
     117 
     118      zundsat(:,:) = pwcp(:,:,jwoxy) 
     119 
     120      DO jk = 2, jpksed 
     121         DO ji = 1, jpoce 
     122            zlimo2(ji,jk) = 1.0 / ( zundsat(ji,jk) + xksedo2 ) 
     123            zsolid1 = zvolc(ji,jk,jspoc)  * solcp(ji,jk,jspoc) 
     124            zsolid2 = zvolc(ji,jk,jspos)  * solcp(ji,jk,jspos) 
     125            zsolid3 = zvolc(ji,jk,jspor)  * solcp(ji,jk,jspor) 
     126            zkpoca  = zkpoc(ji,jk) * zlimo2(ji,jk) 
     127            zkpocb  = zkpos(ji,jk) * zlimo2(ji,jk) 
     128            zkpocc  = zkpor(ji,jk) * zlimo2(ji,jk) 
     129            zrearat1(ji,jk)  = ( zkpoc(ji,jk) * dtsed2 * zsolid1 ) / & 
     130            &                 ( 1. + zkpoca * zundsat(ji,jk ) * dtsed2 ) 
     131            zrearat2(ji,jk)  = ( zkpos(ji,jk) * dtsed2 * zsolid2 ) / & 
     132            &                 ( 1. + zkpocb * zundsat(ji,jk ) * dtsed2 ) 
     133            zrearat3(ji,jk)  = ( zkpor(ji,jk) * dtsed2 * zsolid3 ) / & 
     134            &                 ( 1. + zkpocc * zundsat(ji,jk ) * dtsed2 ) 
     135         ENDDO 
     136      ENDDO 
    199137 
    200138      ! left hand side of coefficient matrix 
    201       DO jk = 2, jpksed 
    202          DO ji = 1, jpoce 
    203             zsolid1 = zvolc(ji,jk,jsopal) * solcp(ji,jk,jsopal) 
    204             zsolid2 = zvolc(ji,jk,jsclay) * solcp(ji,jk,jsclay) 
    205             zsolid3 = zvolc(ji,jk,jspoc)  * solcp(ji,jk,jspoc) 
    206  
    207             zrearat(ji,jk,1)  = ( reac_sil * dtsed * zsolid1 ) / & 
    208                &                ( 1. + reac_sil * dtsed * zundsat(ji,jk,1 ) ) 
    209             zrearat(ji,jk,2)  = ( reac_clay * dtsed * zsolid2 ) / & 
    210                &                ( 1. + reac_clay * dtsed * zundsat(ji,jk,2 ) ) 
    211             zrearat(ji,jk,3)  = ( reac_poc  * dtsed * zsolid3 ) / & 
    212                &                ( 1. + reac_poc  * dtsed * zundsat(ji,jk,3 ) ) 
    213          ENDDO 
    214       ENDDO 
    215  
    216  
    217       CALL sed_mat( nv, jpoce, jpksed, zrearat, zundsat ) 
    218  
     139!      DO jn = 1, 5 
     140      DO jk = 2, jpksed 
     141         DO ji = 1, jpoce 
     142jflag1:     DO jn = 1, 10 
     143               zsolid1 = zvolc(ji,jk,jspoc)  * solcp(ji,jk,jspoc) 
     144               zsolid2 = zvolc(ji,jk,jspos)  * solcp(ji,jk,jspos) 
     145               zsolid3 = zvolc(ji,jk,jspor)  * solcp(ji,jk,jspor) 
     146               zbeta   = xksedo2 - pwcp(ji,jk,jwoxy) + so2ut * ( zrearat1(ji,jk)    & 
     147               &         + zrearat2(ji,jk) + zrearat3(ji,jk) ) 
     148               zgamma = - xksedo2 * pwcp(ji,jk,jwoxy) 
     149               zundsat2 = zundsat(ji,jk) 
     150               zundsat(ji,jk) = ( - zbeta + SQRT( zbeta**2 - 4.0 * zgamma ) ) / 2.0 
     151               zlimo2(ji,jk) = 1.0 / ( zundsat(ji,jk) + xksedo2 ) 
     152               zkpoca  = zkpoc(ji,jk) * zlimo2(ji,jk) 
     153               zkpocb  = zkpos(ji,jk) * zlimo2(ji,jk) 
     154               zkpocc  = zkpor(ji,jk) * zlimo2(ji,jk) 
     155               zrearat1(ji,jk)  = ( zkpoc(ji,jk) * dtsed2 * zsolid1 ) / & 
     156               &                 ( 1. + zkpoca * zundsat(ji,jk ) * dtsed2 ) 
     157               zrearat2(ji,jk)  = ( zkpos(ji,jk) * dtsed2 * zsolid2 ) / & 
     158               &                 ( 1. + zkpocb * zundsat(ji,jk ) * dtsed2 ) 
     159               zrearat3(ji,jk)  = ( zkpor(ji,jk) * dtsed2 * zsolid3 ) / & 
     160               &                 ( 1. + zkpocc * zundsat(ji,jk ) * dtsed2 ) 
     161               IF ( ABS( (zundsat(ji,jk)-zundsat2)/(zundsat2+rtrn)) < 1E-8 ) THEN 
     162                  EXIT jflag1 
     163               ENDIF 
     164            END DO jflag1 
     165         END DO 
     166      END DO 
    219167 
    220168      ! New solid concentration values (jk=2 to jksed) for each couple  
    221       DO js = 1, nv 
    222          DO jk = 2, jpksed 
    223             DO ji = 1, jpoce 
    224                zreasat = zrearat(ji,jk,js) * zundsat(ji,jk,js) / zvolc(ji,jk,js) 
    225                solcp(ji,jk,js) = solcp(ji,jk,js) - zreasat 
    226             ENDDO 
    227          ENDDO 
    228       ENDDO 
    229       ! mass of O2/NO3 before POC remin. for mass balance check  
    230       ! det. of o2 consomation/NO3 production Mc13 
    231       DO jk = 1, jpksed 
    232          DO ji = 1, jpoce 
    233             zvolw = volw3d(ji,jk) * 1.e-3 
    234             zmo2_0 (ji)  = zmo2_0 (ji) + pwcp(ji,jk,jwoxy) * zvolw 
    235             zmno3_0(ji)  = zmno3_0(ji) + pwcp(ji,jk,jwno3) * zvolw 
    236             zmc13_0(ji)  = zmc13_0(ji) + pwcp(ji,jk,jwc13) * zvolw 
     169      DO jk = 2, jpksed 
     170         DO ji = 1, jpoce 
     171            zreasat = zrearat1(ji,jk) * zlimo2(ji,jk) * zundsat(ji,jk) / zvolc(ji,jk,jspoc) 
     172            solcp(ji,jk,jspoc) = solcp(ji,jk,jspoc) - zreasat 
     173            zreasat = zrearat2(ji,jk) * zlimo2(ji,jk) * zundsat(ji,jk) / zvolc(ji,jk,jspos) 
     174            solcp(ji,jk,jspos) = solcp(ji,jk,jspos) - zreasat 
     175            zreasat = zrearat3(ji,jk) * zlimo2(ji,jk) * zundsat(ji,jk) / zvolc(ji,jk,jspor) 
     176            solcp(ji,jk,jspor) = solcp(ji,jk,jspor) - zreasat 
    237177         ENDDO 
    238178      ENDDO 
    239179 
    240180      ! New pore water concentrations     
    241       DO jk = 1, jpksed 
     181      DO jk = 2, jpksed 
    242182         DO ji = 1, jpoce 
    243183            ! Acid Silicic  
    244             pwcp(ji,jk,jwsil)  = sat_sil - zundsat(ji,jk,1)             
    245             ! For O2 (in mol/l) 
    246             pwcp(ji,jk,jwoxy)  = zundsat(ji,jk,3) * so2ut  
    247             zreasat = zrearat(ji,jk,3) * zundsat(ji,jk,3)    ! oxygen          
     184            pwcp(ji,jk,jwoxy)  = zundsat(ji,jk) 
     185            zreasat = ( zrearat1(ji,jk) + zrearat2(ji,jk) + zrearat3(ji,jk) ) * zlimo2(ji,jk) * zundsat(ji,jk)    ! oxygen          
    248186            ! For DIC 
    249187            pwcp(ji,jk,jwdic)  = pwcp(ji,jk,jwdic) + zreasat 
    250             ! For nitrates 
    251             pwcp(ji,jk,jwno3)  = pwcp(ji,jk,jwno3) + zreasat * srno3             
     188            zsumtot(ji) = zsumtot(ji) + zreasat / dtsed2 * volw3d(ji,jk) * 1.e-3 * 86400. * 365. * 1E3 
    252189            ! For Phosphate (in mol/l) 
    253             pwcp(ji,jk,jwpo4)  = pwcp(ji,jk,jwpo4) + zreasat * spo4r             
     190            pwcp(ji,jk,jwpo4)  = pwcp(ji,jk,jwpo4) + zreasat * spo4r 
     191            ! For iron (in mol/l) 
     192            pwcp(ji,jk,jwfe2)  = pwcp(ji,jk,jwfe2) + fecratio(ji) * zreasat 
    254193            ! For alkalinity 
    255             pwcp(ji,jk,jwalk)  = pwcp(ji,jk,jwalk) - zreasat * ( srno3 + 2.* spo4r )            
    256             ! For DIC13 
    257             pwcp(ji,jk,jwc13)  = pwcp(ji,jk,jwc13) + zreasat * rc13P * pdb 
    258          ENDDO 
    259       ENDDO 
    260  
    261  
    262       ! Mass of O2 for mass balance check and det. of o2 consomation 
    263       DO jk = 1, jpksed 
    264          DO ji = 1, jpoce 
    265             zvolw = volw3d(ji,jk) * 1.e-3 
    266             zmo2_1 (ji) = zmo2_1 (ji) + pwcp(ji,jk,jwoxy) * zvolw 
    267             zmno3_1(ji) = zmno3_1(ji) + pwcp(ji,jk,jwno3) * zvolw 
    268             zmc13_1(ji) = zmc13_1(ji) + pwcp(ji,jk,jwc13) * zvolw 
    269          ENDDO 
    270       ENDDO 
    271  
    272       DO ji = 1, jpoce 
    273          cons_o2 (ji) = zmo2_0 (ji) - zmo2_1 (ji) 
    274          sour_no3(ji) = zmno3_1(ji) - zmno3_0(ji)      
    275          sour_c13(ji) = zmc13_1(ji) - zmc13_0(ji)  
    276       ENDDO 
    277   
     194            pwcp(ji,jk,jwalk)  = pwcp(ji,jk,jwalk) + zreasat * ( srno3 * zadsnh4 - 2.* spo4r ) 
     195            ! Ammonium 
     196            pwcp(ji,jk,jwnh4)  = pwcp(ji,jk,jwnh4) + zreasat * srno3 * zadsnh4 
     197            ! Ligands 
     198            sedligand(ji,jk)   = sedligand(ji,jk) + ratligc * zreasat - reac_ligc * sedligand(ji,jk) 
     199         ENDDO 
     200      ENDDO 
    278201 
    279202      !-------------------------------------------------------------------- 
    … …  
    282205      !-------------------------------------------------------------------- 
    283206 
    284       nv = 1 
    285       DO jk = 1, jpksed 
    286          DO ji = 1, jpoce 
    287             zundsat(ji,jk,1) = pwcp(ji,jk,jwno3) / srDnit 
    288          ENDDO 
    289       ENDDO 
    290       DO jk = 2, jpksed 
    291          DO ji = 1, jpoce 
    292             IF( pwcp(ji,jk,jwoxy) < sthrO2 ) THEN 
     207      zrearat1(:,:) = 0. 
     208      zrearat2(:,:) = 0. 
     209      zrearat3(:,:) = 0. 
     210 
     211      zundsat(:,:) = pwcp(:,:,jwno3) 
     212 
     213      DO jk = 2, jpksed 
     214         DO ji = 1, jpoce 
     215            zlimno3(ji,jk) = ( 1.0 - pwcp(ji,jk,jwoxy) * zlimo2(ji,jk) ) / ( zundsat(ji,jk) + xksedno3 ) 
     216            zsolid1 = zvolc(ji,jk,jspoc) * solcp(ji,jk,jspoc) 
     217            zsolid2 = zvolc(ji,jk,jspos) * solcp(ji,jk,jspos) 
     218            zsolid3 = zvolc(ji,jk,jspor) * solcp(ji,jk,jspor) 
     219            zkpoca = zkpoc(ji,jk) * zlimno3(ji,jk) 
     220            zkpocb = zkpos(ji,jk) * zlimno3(ji,jk) 
     221            zkpocc = zkpor(ji,jk) * zlimno3(ji,jk) 
     222            zrearat1(ji,jk)  = ( zkpoc(ji,jk) * dtsed2 * zsolid1 ) / & 
     223            &                 ( 1. + zkpoca * zundsat(ji,jk ) * dtsed2 ) 
     224            zrearat2(ji,jk)  = ( zkpos(ji,jk) * dtsed2 * zsolid2 ) / & 
     225            &                 ( 1. + zkpocb * zundsat(ji,jk ) * dtsed2 ) 
     226            zrearat3(ji,jk)  = ( zkpor(ji,jk) * dtsed2 * zsolid3 ) / & 
     227            &                 ( 1. + zkpocc * zundsat(ji,jk ) * dtsed2 ) 
     228        END DO 
     229      END DO 
     230 
     231!      DO jn = 1, 5 
     232      DO jk = 2, jpksed 
     233         DO ji = 1, jpoce 
     234jflag2:    DO jn = 1, 10 
     235               zlimtmp = ( 1.0 - pwcp(ji,jk,jwoxy) * zlimo2(ji,jk) ) 
    293236               zsolid1 = zvolc(ji,jk,jspoc) * solcp(ji,jk,jspoc) 
    294                zrearat(ji,jk,1) = ( reac_no3 * dtsed * zsolid1 ) / & 
    295                   &                    ( 1. + reac_no3 * dtsed * zundsat(ji,jk,1 ) ) 
    296             ELSE 
    297                zrearat(ji,jk,1) = 0. 
    298             ENDIF 
    299          END DO 
    300       END DO 
    301  
    302  
    303       ! solves tridiagonal system 
    304       CALL sed_mat( nv, jpoce, jpksed, zrearat, zundsat ) 
     237               zsolid2 = zvolc(ji,jk,jspos) * solcp(ji,jk,jspos) 
     238               zsolid3 = zvolc(ji,jk,jspor) * solcp(ji,jk,jspor) 
     239               zbeta   = xksedno3 - pwcp(ji,jk,jwno3) + srDnit * ( zrearat1(ji,jk)    & 
     240               &         + zrearat2(ji,jk) + zrearat3(ji,jk) ) * zlimtmp 
     241               zgamma = - xksedno3 * pwcp(ji,jk,jwno3) 
     242               zundsat2 = zundsat(ji,jk) 
     243               zundsat(ji,jk) = ( - zbeta + SQRT( zbeta**2 - 4.0 * zgamma ) ) / 2.0 
     244               zlimno3(ji,jk) = ( 1.0 - pwcp(ji,jk,jwoxy) * zlimo2(ji,jk) ) / ( zundsat(ji,jk) + xksedno3 ) 
     245               zkpoca  = zkpoc(ji,jk) * zlimno3(ji,jk) 
     246               zkpocb  = zkpos(ji,jk) * zlimno3(ji,jk) 
     247               zkpocc  = zkpor(ji,jk) * zlimno3(ji,jk) 
     248               zrearat1(ji,jk)  = ( zkpoc(ji,jk) * dtsed2 * zsolid1 ) / & 
     249               &                 ( 1. + zkpoca * zundsat(ji,jk ) * dtsed2 ) 
     250               zrearat2(ji,jk)  = ( zkpos(ji,jk) * dtsed2 * zsolid2 ) / & 
     251               &                 ( 1. + zkpocb * zundsat(ji,jk ) * dtsed2 ) 
     252               zrearat3(ji,jk)  = ( zkpor(ji,jk) * dtsed2 * zsolid3 ) / & 
     253               &                 ( 1. + zkpocc * zundsat(ji,jk ) * dtsed2 ) 
     254               IF ( ABS( (zundsat(ji,jk)-zundsat2)/(zundsat2+rtrn)) < 1E-8 ) THEN 
     255                  EXIT jflag2 
     256               ENDIF 
     257            END DO jflag2 
     258         END DO 
     259      END DO 
    305260 
    306261 
    … …  
    308263      DO jk = 2, jpksed 
    309264         DO ji = 1, jpoce 
    310             zreasat = zrearat(ji,jk,1) * zundsat(ji,jk,1) / zvolc(ji,jk,jspoc) 
     265            zreasat = zrearat1(ji,jk) * zlimno3(ji,jk) * zundsat(ji,jk) / zvolc(ji,jk,jspoc) 
    311266            solcp(ji,jk,jspoc) = solcp(ji,jk,jspoc) - zreasat 
     267            zreasat = zrearat2(ji,jk) * zlimno3(ji,jk) * zundsat(ji,jk) / zvolc(ji,jk,jspos) 
     268            solcp(ji,jk,jspos) = solcp(ji,jk,jspos) - zreasat 
     269            zreasat = zrearat3(ji,jk) * zlimno3(ji,jk) * zundsat(ji,jk) / zvolc(ji,jk,jspor) 
     270            solcp(ji,jk,jspor) = solcp(ji,jk,jspor) - zreasat 
    312271         ENDDO 
    313272      ENDDO 
    314273 
    315274      ! New dissolved concentrations 
    316       DO jk = 1, jpksed 
    317          DO ji = 1, jpoce 
    318             zreasat = zrearat(ji,jk,1) * zundsat(ji,jk,1)     
     275      DO jk = 2, jpksed 
     276         DO ji = 1, jpoce 
    319277            ! For nitrates 
    320             pwcp(ji,jk,jwno3)  =  zundsat(ji,jk,1) * srDnit 
     278            pwcp(ji,jk,jwno3)  =  zundsat(ji,jk) 
     279            zreasat = ( zrearat1(ji,jk) + zrearat2(ji,jk) + zrearat3(ji,jk) ) * zlimno3(ji,jk) * zundsat(ji,jk) 
    321280            ! For DIC 
    322281            pwcp(ji,jk,jwdic)  = pwcp(ji,jk,jwdic) + zreasat 
     282            zsumtot(ji) = zsumtot(ji) + zreasat / dtsed2 * volw3d(ji,jk) * 1.e-3 * 86400. * 365. * 1E3 
    323283            ! For Phosphate (in mol/l) 
    324284            pwcp(ji,jk,jwpo4)  = pwcp(ji,jk,jwpo4) + zreasat * spo4r             
     285            ! Ligands 
     286            sedligand(ji,jk)   = sedligand(ji,jk) + ratligc * zreasat 
     287            ! For iron (in mol/l) 
     288            pwcp(ji,jk,jwfe2)  = pwcp(ji,jk,jwfe2) + fecratio(ji) * zreasat 
    325289            ! For alkalinity 
    326             pwcp(ji,jk,jwalk)  = pwcp(ji,jk,jwalk) + zreasat * ( srDnit - 2.* spo4r )            
    327             ! For DIC13 
    328             pwcp(ji,jk,jwc13)  = pwcp(ji,jk,jwc13) + zreasat * rc13P * pdb 
    329          ENDDO 
    330       ENDDO 
    331  
    332  
    333       ! Mass of O2 for mass balance check and det. of o2 consomation 
    334       DO jk = 1, jpksed 
    335          DO ji = 1, jpoce 
    336             zvolw = volw3d(ji,jk) * 1.e-3 
    337             zmno3_2(ji) = zmno3_2(ji) + pwcp(ji,jk ,jwno3) * zvolw 
    338             zmc13_2(ji) = zmc13_2(ji) + pwcp(ji,jk ,jwc13) * zvolw    
    339          ENDDO 
    340       ENDDO 
    341  
    342       DO ji = 1, jpoce 
    343          cons_no3(ji) = zmno3_1(ji) - zmno3_2(ji)   
    344          sour_c13(ji) = sour_c13(ji) + zmc13_2(ji) - zmc13_1(ji)      
    345       ENDDO 
    346  
    347  
    348       !--------------------------- 
    349       ! Solves PO4 diffusion  
    350       !---------------------------- 
    351  
    352       nv = 1 
    353       DO jk = 1, jpksed 
    354          DO ji = 1, jpoce 
    355             zundsat(ji,jk,1) = pwcp(ji,jk,jwpo4) 
    356             zrearat(ji,jk,1) = 0. 
    357          ENDDO 
    358       ENDDO 
    359  
    360  
    361       ! solves tridiagonal system 
    362       CALL sed_mat( nv, jpoce, jpksed, zrearat, zundsat ) 
    363  
    364  
    365       ! New undsaturation values and dissolved concentrations 
    366       DO jk = 1, jpksed 
    367          DO ji = 1, jpoce 
    368             pwcp(ji,jk,jwpo4) = zundsat(ji,jk,1) 
    369          ENDDO 
    370       ENDDO 
    371  
    372  
    373       !--------------------------------------------------------------- 
    374       ! Performs CaCO3 particle deposition and redissolution (indice 9) 
    375       !-------------------------------------------------------------- 
    376  
    377       ! computes co3por from the updated pwcp concentrations (note [co3por] = mol/l) 
    378  
    379       CALL sed_co3( kt ) 
    380  
    381  
    382       nv = 1 
    383       ! *densSW(l)**2 converts aksps [mol2/kg sol2] into [mol2/l2] to get [undsat] in [mol/l] 
    384       DO jk = 1, jpksed 
    385          DO ji = 1, jpoce 
    386             zundsat(ji,jk,1) = aksps(ji) * densSW(ji) * densSW(ji) / calcon2(ji) & 
    387                &                     - co3por(ji,jk) 
    388             ! positive values of undersaturation 
    389             zundsat(ji,jk,1) = MAX( 0., zundsat(ji,jk,1) )             
    390          ENDDO 
    391       ENDDO 
    392  
    393       DO jk = 2, jpksed 
    394          DO ji = 1, jpoce 
    395             zsolid1 = zvolc(ji,jk,jscal) * solcp(ji,jk,jscal) 
    396             zrearat(ji,jk,1) = ( reac_cal * dtsed * zsolid1 ) / & 
    397                   &               ( 1. + reac_cal * dtsed * zundsat(ji,jk,1) ) 
    398          END DO 
    399       END DO 
    400  
    401  
    402       ! solves tridiagonal system 
    403       CALL sed_mat( nv, jpoce, jpksed, zrearat, zundsat ) 
    404  
    405  
    406       ! New solid concentration values (jk=2 to jksed) for cacO3 
    407       DO jk = 2, jpksed 
    408          DO ji = 1, jpoce 
    409             zreasat = zrearat(ji,jk,1) * zundsat(ji,jk,1) / zvolc(ji,jk,jscal) 
    410             solcp(ji,jk,jscal) = solcp(ji,jk,jscal) - zreasat 
    411          ENDDO 
    412       ENDDO 
     290            pwcp(ji,jk,jwalk)  = pwcp(ji,jk,jwalk) + zreasat * ( srDnit + srno3 * zadsnh4 - 2.* spo4r )            
     291            ! Ammonium 
     292            pwcp(ji,jk,jwnh4)  = pwcp(ji,jk,jwnh4) + zreasat * srno3 * zadsnh4 
     293         ENDDO 
     294      ENDDO 
     295 
     296      !-------------------------------------------------------------------- 
     297      ! Begining POC iron reduction 
     298      ! (indice n�5 for couple POFe(OH)3 ie solcp(:,:,jspoc)/pwcp(:,:,jsfeo)) 
     299      !-------------------------------------------------------------------- 
     300 
     301      zrearat1(:,:) = 0. 
     302      zrearat2(:,:) = 0. 
     303      zrearat3(:,:) = 0. 
     304 
     305      zundsat(:,:) = solcp(:,:,jsfeo) 
     306 
     307      DO jk = 2, jpksed 
     308         DO ji = 1, jpoce 
     309            zlimfeo(ji,jk) = ( 1.0 - pwcp(ji,jk,jwoxy) * zlimo2(ji,jk) ) * ( 1.0 - pwcp(ji,jk,jwno3)    & 
     310            &                / ( pwcp(ji,jk,jwno3) + xksedno3 ) ) / ( zundsat(ji,jk) + xksedfeo ) 
     311            zsolid1 = zvolc(ji,jk,jspoc) * solcp(ji,jk,jspoc) 
     312            zsolid2 = zvolc(ji,jk,jspos) * solcp(ji,jk,jspos) 
     313            zsolid3 = zvolc(ji,jk,jspor) * solcp(ji,jk,jspor) 
     314            zkpoca = zkpoc(ji,jk) * zlimfeo(ji,jk) 
     315            zkpocb = zkpos(ji,jk) * zlimfeo(ji,jk) 
     316            zkpocc = zkpor(ji,jk) * zlimfeo(ji,jk) 
     317            zrearat1(ji,jk) = ( zkpoc(ji,jk) * dtsed2 * zsolid1 ) / & 
     318            &                    ( 1. + zkpoca * zundsat(ji,jk) * dtsed2 ) 
     319            zrearat2(ji,jk) = ( zkpos(ji,jk) * dtsed2 * zsolid2 ) / & 
     320            &                    ( 1. + zkpocb * zundsat(ji,jk) * dtsed2 ) 
     321            zrearat3(ji,jk) = ( zkpor(ji,jk) * dtsed2 * zsolid3 ) / & 
     322            &                    ( 1. + zkpocc * zundsat(ji,jk) * dtsed2 ) 
     323         END DO 
     324      END DO 
     325 
     326!      DO jn = 1, 5 
     327      DO jk = 2, jpksed 
     328         DO ji = 1, jpoce 
     329jflag3:     DO jn = 1, 10 
     330               zlimtmp = ( 1.0 - pwcp(ji,jk,jwoxy) * zlimo2(ji,jk) ) * ( 1.0 - pwcp(ji,jk,jwno3)    & 
     331               &                / ( pwcp(ji,jk,jwno3) + xksedno3 ) ) 
     332               zsolid1 = zvolc(ji,jk,jspoc) * solcp(ji,jk,jspoc) 
     333               zsolid2 = zvolc(ji,jk,jspos) * solcp(ji,jk,jspos) 
     334               zsolid3 = zvolc(ji,jk,jspor) * solcp(ji,jk,jspor) 
     335               zreasat = ( zrearat1(ji,jk) + zrearat2(ji,jk) + zrearat3(ji,jk) ) / zvolc(ji,jk,jsfeo) 
     336               zbeta   = xksedfeo - solcp(ji,jk,jsfeo) + 4.0 * zreasat * zlimtmp 
     337               zgamma  = -xksedfeo * solcp(ji,jk,jsfeo) 
     338               zundsat2 = zundsat(ji,jk) 
     339               zundsat(ji,jk) = ( - zbeta + SQRT( zbeta**2 - 4.0 * zgamma ) ) / 2.0 
     340               zlimfeo(ji,jk) = ( 1.0 - pwcp(ji,jk,jwoxy) * zlimo2(ji,jk) ) * ( 1.0 - pwcp(ji,jk,jwno3)    & 
     341               &                / ( pwcp(ji,jk,jwno3) + xksedno3 ) ) / ( zundsat(ji,jk) + xksedfeo ) 
     342               zkpoca  = zkpoc(ji,jk) * zlimfeo(ji,jk) 
     343               zkpocb  = zkpos(ji,jk) * zlimfeo(ji,jk) 
     344               zkpocc  = zkpor(ji,jk) * zlimfeo(ji,jk) 
     345               zrearat1(ji,jk) = ( zkpoc(ji,jk) * dtsed2 * zsolid1 ) / & 
     346               &                    ( 1. + zkpoca * zundsat(ji,jk) * dtsed2 ) 
     347               zrearat2(ji,jk) = ( zkpos(ji,jk) * dtsed2 * zsolid2 ) / & 
     348               &                    ( 1. + zkpocb * zundsat(ji,jk) * dtsed2 ) 
     349               zrearat3(ji,jk) = ( zkpor(ji,jk) * dtsed2 * zsolid3 ) / & 
     350               &                    ( 1. + zkpocc * zundsat(ji,jk) * dtsed2 ) 
     351               IF ( ABS( (zundsat(ji,jk)-zundsat2)/( MAX(0.,zundsat2)+rtrn)) < 1E-8 ) THEN 
     352                  EXIT jflag3 
     353               ENDIF 
     354            END DO jflag3 
     355         END DO 
     356      END DO 
     357 
     358 
     359         ! New solid concentration values (jk=2 to jksed) for each couple  
     360      DO jk = 2, jpksed 
     361         DO ji = 1, jpoce 
     362            zreasat = zrearat1(ji,jk) * zlimfeo(ji,jk) * zundsat(ji,jk) / zvolc(ji,jk,jspoc) 
     363            solcp(ji,jk,jspoc) = solcp(ji,jk,jspoc) - zreasat 
     364            zreasat = zrearat2(ji,jk) * zlimfeo(ji,jk) * zundsat(ji,jk) / zvolc(ji,jk,jspos) 
     365            solcp(ji,jk,jspos) = solcp(ji,jk,jspos) - zreasat 
     366            zreasat = zrearat3(ji,jk) * zlimfeo(ji,jk) * zundsat(ji,jk) / zvolc(ji,jk,jspor) 
     367            solcp(ji,jk,jspor) = solcp(ji,jk,jspor) - zreasat 
     368         END DO 
     369      END DO 
    413370 
    414371      ! New dissolved concentrations 
    415       DO jk = 1, jpksed 
    416          DO ji = 1, jpoce 
    417             zreasat = zrearat(ji,jk,1) * zundsat(ji,jk,1)     
     372      DO jk = 2, jpksed 
     373         DO ji = 1, jpoce 
     374            zreasat = ( zrearat1(ji,jk) + zrearat2(ji,jk) + zrearat3(ji,jk) ) * zlimfeo(ji,jk) * zundsat(ji,jk) 
     375            ! For FEOH 
     376            solcp(ji,jk,jsfeo) = zundsat(ji,jk) 
    418377            ! For DIC 
    419378            pwcp(ji,jk,jwdic)  = pwcp(ji,jk,jwdic) + zreasat 
     379            zsumtot(ji) = zsumtot(ji) + zreasat / dtsed2 * volw3d(ji,jk) * 1.e-3 * 86400. * 365. * 1E3 
     380            ! For Phosphate (in mol/l) 
     381            pwcp(ji,jk,jwpo4)  = pwcp(ji,jk,jwpo4) + zreasat * ( spo4r + 4.0 * redfep ) 
     382            ! Ligands 
     383            sedligand(ji,jk)   = sedligand(ji,jk) + ratligc * zreasat 
     384            ! For iron (in mol/l) 
     385            pwcp(ji,jk,jwfe2)  = pwcp(ji,jk,jwfe2) + fecratio(ji) * zreasat 
    420386            ! For alkalinity 
    421             pwcp(ji,jk,jwalk)  = pwcp(ji,jk,jwalk) + 2.* zreasat  
    422             ! For DIC13 
    423             pwcp(ji,jk,jwc13)  = pwcp(ji,jk,jwc13) + zreasat * rc13Ca * pdb 
    424          ENDDO 
    425       ENDDO 
    426  
    427       DO jk = 1, jpksed 
    428          DO ji = 1, jpoce 
    429             zmc13_3(ji) = zmc13_3(ji) + pwcp(ji,jk,jwc13) * volw3d(ji,jk) * 1.e-3                   
    430          ENDDO 
    431       ENDDO 
    432  
    433       DO ji = 1, jpoce      
    434          sour_c13(ji) = sour_c13(ji) + zmc13_3(ji) - zmc13_2(ji)    
    435       ENDDO 
    436        
    437       !------------------------------------------------- 
    438       ! Beginning DIC, Alkalinity and DIC13 diffusion 
    439       !------------------------------------------------- 
    440        
    441       nv = 3 
    442       DO jk = 1, jpksed 
    443          DO ji = 1, jpoce       
    444             zundsat(ji,jk,1) = pwcp(ji,jk,jwdic) 
    445             zundsat(ji,jk,2) = pwcp(ji,jk,jwalk) 
    446             zundsat(ji,jk,3) = pwcp(ji,jk,jwc13) 
    447        
    448             zrearat(ji,jk,1) = 0. 
    449             zrearat(ji,jk,2) = 0. 
    450             zrearat(ji,jk,3) = 0. 
    451        
    452          ENDDO 
    453       ENDDO 
    454  
    455  
    456       ! solves tridiagonal system 
    457       CALL sed_mat( nv, jpoce, jpksed, zrearat, zundsat ) 
    458  
    459  
    460       ! New dissolved concentrations       
    461       DO jk = 1, jpksed 
    462          DO ji = 1, jpoce                       
    463             pwcp(ji,jk,jwdic) = zundsat(ji,jk,1) 
    464             pwcp(ji,jk,jwalk) = zundsat(ji,jk,2) 
    465             pwcp(ji,jk,jwc13) = zundsat(ji,jk,3) 
    466          ENDDO 
    467       ENDDO             
    468        
    469       !---------------------------------- 
    470       !   Back to initial geometry 
    471       !----------------------------- 
    472        
    473       !--------------------------------------------------------------------- 
    474       !   1/ Compensation for ajustement of the bottom water concentrations 
    475       !      (see note n° 1 about *por(2)) 
     387            pwcp(ji,jk,jwalk)  = pwcp(ji,jk,jwalk) + zreasat * ( srno3 * zadsnh4 - 2.* spo4r ) + 8.0 * zreasat 
     388            ! Ammonium 
     389            pwcp(ji,jk,jwnh4)  = pwcp(ji,jk,jwnh4) + zreasat * srno3 * zadsnh4 
     390            pwcp(ji,jk,jwfe2)  = pwcp(ji,jk,jwfe2) + zreasat * 4.0 
     391         ENDDO 
     392      ENDDO 
     393 
    476394      !-------------------------------------------------------------------- 
    477       DO jw = 1, jpwat 
    478          DO ji = 1, jpoce 
    479             pwcp(ji,1,jw) = pwcp(ji,1,jw) + & 
    480                &            pwcp(ji,2,jw) * dzdep(ji) * por(2) / dzkbot(ji) 
    481          END DO 
    482       ENDDO 
    483        
    484       !----------------------------------------------------------------------- 
    485       !    2/ Det of new rainrg taking account of the new weight fraction obtained  
    486       !      in dz3d(2) after diffusion/reaction (react/diffu are also in dzdep!) 
    487       !      This new rain (rgntg rm) will be used in advection/burial routine 
    488       !------------------------------------------------------------------------ 
    489       DO js = 1, jpsol 
    490          DO ji = 1, jpoce 
    491             rainrg(ji,js) = raintg(ji) * solcp(ji,2,js) 
    492             rainrm(ji,js) = rainrg(ji,js) / mol_wgt(js) 
    493          END DO 
    494       ENDDO 
    495        
    496       !  New raintg 
    497       raintg(:) = 0. 
    498       DO js = 1, jpsol 
    499          DO ji = 1, jpoce 
    500             raintg(ji) = raintg(ji) + rainrg(ji,js) 
    501          END DO 
    502       ENDDO 
    503        
    504       !-------------------------------- 
    505       !    3/ back to initial geometry 
    506       !-------------------------------- 
     395      ! Begining POC denitrification and NO3- diffusion 
     396      ! (indice n�5 for couple POC/NO3- ie solcp(:,:,jspoc)/pwcp(:,:,jwno3)) 
     397      !-------------------------------------------------------------------- 
     398 
     399      zrearat1(:,:) = 0. 
     400      zrearat2(:,:) = 0. 
     401      zrearat3(:,:) = 0. 
     402 
     403      zundsat(:,:) = pwcp(:,:,jwso4) 
     404 
     405      DO jk = 2, jpksed 
     406         DO ji = 1, jpoce 
     407            zlimso4(ji,jk) = ( 1.0 - pwcp(ji,jk,jwoxy) * zlimo2(ji,jk) ) * ( 1.0 - pwcp(ji,jk,jwno3)    & 
     408            &                / ( pwcp(ji,jk,jwno3) + xksedno3 ) ) * ( 1. - solcp(ji,jk,jsfeo)  & 
     409            &                / ( solcp(ji,jk,jsfeo) + xksedfeo ) ) / ( zundsat(ji,jk) + xksedso4 ) 
     410            zsolid1 = zvolc(ji,jk,jspoc) * solcp(ji,jk,jspoc) 
     411            zsolid2 = zvolc(ji,jk,jspos) * solcp(ji,jk,jspos) 
     412            zsolid3 = zvolc(ji,jk,jspor) * solcp(ji,jk,jspor) 
     413            zkpoca = zkpoc(ji,jk) * zlimso4(ji,jk) 
     414            zkpocb = zkpos(ji,jk) * zlimso4(ji,jk) 
     415            zkpocc = zkpor(ji,jk) * zlimso4(ji,jk) 
     416            zrearat1(ji,jk)  = ( zkpoc(ji,jk) * dtsed2 * zsolid1 ) / & 
     417            &                 ( 1. + zkpoca * zundsat(ji,jk ) * dtsed2 ) 
     418            zrearat2(ji,jk)  = ( zkpos(ji,jk) * dtsed2 * zsolid2 ) / & 
     419            &                 ( 1. + zkpocb * zundsat(ji,jk ) * dtsed2 ) 
     420            zrearat3(ji,jk)  = ( zkpor(ji,jk) * dtsed2 * zsolid3 ) / & 
     421            &                 ( 1. + zkpocc * zundsat(ji,jk ) * dtsed2 ) 
     422        END DO 
     423      END DO 
     424! 
     425!      DO jn = 1, 5  
     426      DO jk = 2, jpksed 
     427         DO ji = 1, jpoce 
     428jflag4:     DO jn = 1, 10 
     429               zlimtmp = ( 1.0 - pwcp(ji,jk,jwoxy) * zlimo2(ji,jk) ) * ( 1.0 - pwcp(ji,jk,jwno3)    & 
     430               &         / ( pwcp(ji,jk,jwno3) + xksedno3 ) ) * ( 1. - solcp(ji,jk,jsfeo)  & 
     431               &         / ( solcp(ji,jk,jsfeo) + xksedfeo ) )  
     432               zsolid1 = zvolc(ji,jk,jspoc) * solcp(ji,jk,jspoc) 
     433               zsolid2 = zvolc(ji,jk,jspos) * solcp(ji,jk,jspos) 
     434               zsolid3 = zvolc(ji,jk,jspor) * solcp(ji,jk,jspor) 
     435               zreasat = ( zrearat1(ji,jk) + zrearat2(ji,jk) + zrearat3(ji,jk) )  
     436               zbeta   = xksedso4 - pwcp(ji,jk,jwso4) + 0.5 * zreasat * zlimtmp 
     437               zgamma = - xksedso4 * pwcp(ji,jk,jwso4) 
     438               zundsat2 = zundsat(ji,jk) 
     439               zundsat(ji,jk) = ( - zbeta + SQRT( zbeta**2 - 4.0 * zgamma ) ) / 2.0 
     440               zlimso4(ji,jk) = ( 1.0 - pwcp(ji,jk,jwoxy) * zlimo2(ji,jk) ) * ( 1.0 - pwcp(ji,jk,jwno3)    & 
     441               &                / ( pwcp(ji,jk,jwno3) + xksedno3 ) ) * ( 1. - solcp(ji,jk,jsfeo)  & 
     442               &                / ( solcp(ji,jk,jsfeo) + xksedfeo ) ) / ( zundsat(ji,jk) + xksedso4 ) 
     443               zkpoca  = zkpoc(ji,jk) * zlimso4(ji,jk) 
     444               zkpocb  = zkpos(ji,jk) * zlimso4(ji,jk) 
     445               zkpocc  = zkpor(ji,jk) * zlimso4(ji,jk) 
     446               zrearat1(ji,jk)  = ( zkpoc(ji,jk) * dtsed2 * zsolid1 ) / & 
     447               &                 ( 1. + zkpoca * zundsat(ji,jk ) * dtsed2 ) 
     448               zrearat2(ji,jk)  = ( zkpos(ji,jk) * dtsed2 * zsolid2 ) / & 
     449               &                 ( 1. + zkpocb * zundsat(ji,jk ) * dtsed2 ) 
     450               zrearat3(ji,jk)  = ( zkpor(ji,jk) * dtsed2 * zsolid3 ) / & 
     451               &                 ( 1. + zkpocc * zundsat(ji,jk ) * dtsed2 ) 
     452               IF ( ABS( (zundsat(ji,jk)-zundsat2)/(zundsat2+rtrn)) < 1E-8 ) THEN 
     453                  EXIT jflag4 
     454               ENDIF 
     455            END DO jflag4 
     456         END DO 
     457      END DO 
     458 
     459     ! New solid concentration values (jk=2 to jksed) for each couple  
     460      DO jk = 2, jpksed 
     461         DO ji = 1, jpoce 
     462            zreasat = zrearat1(ji,jk) * zlimso4(ji,jk) * zundsat(ji,jk) / zvolc(ji,jk,jspoc) 
     463            solcp(ji,jk,jspoc) = solcp(ji,jk,jspoc) - zreasat 
     464            zreasat = zrearat2(ji,jk) * zlimso4(ji,jk) * zundsat(ji,jk) / zvolc(ji,jk,jspos) 
     465            solcp(ji,jk,jspos) = solcp(ji,jk,jspos) - zreasat 
     466            zreasat = zrearat3(ji,jk) * zlimso4(ji,jk) * zundsat(ji,jk) / zvolc(ji,jk,jspor) 
     467            solcp(ji,jk,jspor) = solcp(ji,jk,jspor) - zreasat 
     468         ENDDO 
     469      ENDDO 
     470! 
     471      ! New dissolved concentrations 
     472      DO jk = 2, jpksed 
     473         DO ji = 1, jpoce 
     474            ! For sulfur 
     475            pwcp(ji,jk,jwh2s)  = pwcp(ji,jk,jwh2s) - ( zundsat(ji,jk) - pwcp(ji,jk,jwso4) )  
     476            pwcp(ji,jk,jwso4)  =  zundsat(ji,jk) 
     477            zreasat = ( zrearat1(ji,jk) + zrearat2(ji,jk) + zrearat3(ji,jk) ) * zlimso4(ji,jk) * zundsat(ji,jk) 
     478            ! For DIC 
     479            pwcp(ji,jk,jwdic)  = pwcp(ji,jk,jwdic) + zreasat 
     480            zsumtot(ji) = zsumtot(ji) + zreasat / dtsed2 * volw3d(ji,jk) * 1.e-3 * 86400. * 365. * 1E3 
     481            ! For Phosphate (in mol/l) 
     482            pwcp(ji,jk,jwpo4)  = pwcp(ji,jk,jwpo4) + zreasat * spo4r 
     483            ! Ligands 
     484            sedligand(ji,jk)   = sedligand(ji,jk) + ratligc * zreasat 
     485            ! For iron (in mol/l) 
     486            pwcp(ji,jk,jwfe2)  = pwcp(ji,jk,jwfe2) + fecratio(ji) * zreasat 
     487            ! For alkalinity 
     488            pwcp(ji,jk,jwalk)  = pwcp(ji,jk,jwalk) + zreasat * ( srno3 * zadsnh4 - 2.* spo4r ) + zreasat 
     489            ! Ammonium 
     490            pwcp(ji,jk,jwnh4)  = pwcp(ji,jk,jwnh4) + zreasat * srno3 * zadsnh4 
     491         ENDDO 
     492      ENDDO 
     493 
     494      ! Oxydation of the reduced products. Here only ammonium and ODU is accounted for 
     495      ! There are two options here: A simple time splitting scheme and a modified  
     496     ! Patankar scheme 
     497      ! ------------------------------------------------------------------------------ 
     498 
     499      call sed_dsr_redoxb 
     500 
     501      ! --------------------------------------------------------------  
     502      !    4/ Computation of the bioturbation coefficient 
     503      !       This parameterization is taken from Archer et al. (2002) 
     504      ! -------------------------------------------------------------- 
     505 
    507506      DO ji = 1, jpoce 
    508          dz3d  (ji,2) = dz(2) 
    509          volw3d(ji,2) = dz3d(ji,2) * por(2) 
    510          vols3d(ji,2) = dz3d(ji,2) * por1(2) 
    511       ENDDO 
    512        
    513       !---------------------------------------------------------------------- 
    514       !    4/ Saving new amount of material in dzkbot for mass balance check 
    515       !       tokbot in [mol] (implicit *1cm*1cm for spacial dim) 
    516       !---------------------------------------------------------------------- 
    517       DO jw = 1, jpwat  
    518          DO ji = 1, jpoce 
    519             tokbot(ji,jw) = pwcp(ji,1,jw) * 1.e-3 * dzkbot(ji) 
    520          END DO 
    521       ENDDO 
    522  
    523       DEALLOCATE( zmo2_0  ) ;  DEALLOCATE( zmno3_1 )  ;  DEALLOCATE( zmno3_2 )  
    524       DEALLOCATE( zmc13_0 ) ;  DEALLOCATE( zmc13_1 )  ;  DEALLOCATE( zmc13_2 ) ;  DEALLOCATE( zmc13_3 )  
    525        
    526       DEALLOCATE( zrearat ) ;  DEALLOCATE( zundsat )  ;  DEALLOCATE( zvolc )    
    527        
     507         db(ji,:) = dbiot * zsumtot(ji) * pwcp(ji,1,jwoxy) / (pwcp(ji,1,jwoxy) + 20.E-6) 
     508      END DO 
     509 
     510      ! ------------------------------------------------------ 
     511      !    Vertical variations of the bioturbation coefficient 
     512      ! ------------------------------------------------------ 
     513      IF (ln_btbz) THEN 
     514         DO ji = 1, jpoce 
     515            db(ji,:) = db(ji,:) * exp( -(profsedw(:) / dbtbzsc)**2 ) / (365.0 * 86400.0) 
     516         END DO 
     517      ELSE 
     518         DO jk = 1, jpksed 
     519            IF (profsedw(jk) > dbtbzsc) THEN 
     520                db(:,jk) = 0.0 
     521            ENDIF 
     522         END DO 
     523      ENDIF 
     524 
     525      IF (ln_irrig) THEN 
     526         DO jk = 1, jpksed 
     527            DO ji = 1, jpoce 
     528               irrig(ji,jk) = ( 7.63752 - 7.4465 * exp( -0.89603 * zsumtot(ji) ) ) * pwcp(ji,1,jwoxy)  & 
     529               &             / (pwcp(ji,1,jwoxy) + 20.E-6)  
     530               irrig(ji,jk) = irrig(ji,jk) * exp( -(profsedw(jk) / xirrzsc) ) 
     531            END DO 
     532         END DO 
     533      ELSE 
     534         irrig(:,:) = 0.0 
     535      ENDIF 
     536 
     537      DEALLOCATE( zvolc ) 
     538 
     539      IF( ln_timing )  CALL timing_stop('sed_dsr') 
     540!       
    528541   END SUBROUTINE sed_dsr 
    529 #else 
    530    !!====================================================================== 
    531    !! MODULE seddsr  :   Dummy module  
    532    !!====================================================================== 
    533    !! $Id$ 
    534 CONTAINS 
    535    SUBROUTINE sed_dsr ( kt ) 
    536      INTEGER, INTENT(in) :: kt 
    537      WRITE(*,*) 'sed_dsr: You should not have seen this print! error?', kt 
    538   END SUBROUTINE sed_dsr 
    539 #endif 
    540     
     542 
     543   SUBROUTINE sed_dsr_redoxb 
     544      !!---------------------------------------------------------------------- 
     545      !!                   ***  ROUTINE sed_dsr_redox  *** 
     546      !!  
     547      !!  ** Purpose :  computes secondary redox reactions 
     548      !! 
     549      !!  ** Methode :  It uses Strand splitter algorithm proposed by  
     550      !!                Nguyen et al. (2009) and modified by Wang et al. (2018) 
     551      !!                Basically, each equation is solved analytically when  
     552      !!                feasible, otherwise numerically at t+1/2. Then  
     553      !!                the last equation is solved at t+1. The other equations 
     554      !!                are then solved at t+1 starting in the reverse order. 
     555      !!                Ideally, it's better to start from the fastest reaction 
     556      !!                to the slowest and then reverse the order to finish up 
     557      !!                with the fastest one. But random order works well also. 
     558      !!                The scheme is second order, positive and mass  
     559      !!                conserving. It works well for stiff systems. 
     560      !!  
     561      !!   History : 
     562      !!        !  18-08 (O. Aumont)  Original code 
     563      !!---------------------------------------------------------------------- 
     564      !! Arguments 
     565      ! --- local variables 
     566      INTEGER   ::  ji, jk, jn   ! dummy looop indices 
     567 
     568      REAL, DIMENSION(6)  :: zsedtrn, zsedtra 
     569      REAL(wp)  ::  zalpha, zbeta, zgamma, zdelta, zepsi, zsedfer 
     570      !! 
     571      !!---------------------------------------------------------------------- 
     572 
     573      IF( ln_timing )  CALL timing_start('sed_dsr_redoxb') 
     574 
     575      DO ji = 1, jpoce 
     576         DO jk = 2, jpksed 
     577            zsedtrn(1)  = pwcp(ji,jk,jwoxy) 
     578            zsedtrn(2)  = MAX(0., pwcp(ji,jk,jwh2s) ) 
     579            zsedtrn(3)  = pwcp(ji,jk,jwnh4) 
     580            zsedtrn(4)  = MAX(0., pwcp(ji,jk,jwfe2) - sedligand(ji,jk) ) 
     581            zsedfer     = MIN(0., pwcp(ji,jk,jwfe2) - sedligand(ji,jk) ) 
     582            zsedtrn(5)  = solcp(ji,jk,jsfeo) * zvolc(ji,jk,jsfeo) 
     583            zsedtrn(6)  = solcp(ji,jk,jsfes) * zvolc(ji,jk,jsfes) 
     584            zsedtra(:)  = zsedtrn(:)  
     585 
     586            ! First pass of the scheme. At the end, it is 1st order  
     587            ! ----------------------------------------------------- 
     588            ! Fe + O2 
     589            zalpha = zsedtra(1) - 0.25 * zsedtra(4) 
     590            zbeta  = zsedtra(4) + zsedtra(5)  
     591            zgamma = pwcp(ji,jk,jwalk) - 2.0 * zsedtra(4) 
     592            zdelta = pwcp(ji,jk,jwpo4) - redfep * zsedtra(4) 
     593            zsedtra(4) = ( zsedtra(4) * zalpha ) / ( 0.25 * zsedtra(4) * ( exp( reac_fe2 * zalpha * dtsed2 / 2. ) - 1.0 )  & 
     594            &            + zalpha * exp( reac_fe2 * zalpha * dtsed2 / 2. ) + rtrn ) 
     595            zsedtra(1) = zalpha + 0.25 * zsedtra(4)  
     596            zsedtra(5) = zbeta  - zsedtra(4) 
     597            pwcp(ji,jk,jwalk) = zgamma + 2.0 * zsedtra(4) 
     598            pwcp(ji,jk,jwpo4) = zdelta + redfep * zsedtra(4) 
     599            ! H2S + O2 
     600            zalpha = zsedtra(1) - 2.0 * zsedtra(2) 
     601            zbeta  = pwcp(ji,jk,jwso4) + zsedtra(2) 
     602            zgamma = pwcp(ji,jk,jwalk) - 2.0 * zsedtra(2) 
     603            zsedtra(2) = ( zsedtra(2) * zalpha ) / ( 2.0 * zsedtra(2) * ( exp( reac_h2s * zalpha * dtsed2 / 2. ) - 1.0 )  & 
     604            &            + zalpha * exp( reac_h2s * zalpha * dtsed2 / 2. ) + rtrn ) 
     605            zsedtra(1) = zalpha + 2.0 * zsedtra(2) 
     606            pwcp(ji,jk,jwalk) = zgamma + 2.0 * zsedtra(2) 
     607            pwcp(ji,jk,jwso4) = zbeta - zsedtra(2) 
     608            ! NH4 + O2 
     609            zalpha = zsedtra(1) - 2.0 * zsedtra(3) 
     610            zgamma = pwcp(ji,jk,jwalk) - 2.0 * zsedtra(3) 
     611            zsedtra(3) = ( zsedtra(3) * zalpha ) / ( 2.0 * zsedtra(3) * ( exp( reac_nh4 * zadsnh4 * zalpha * dtsed2 / 2. ) - 1.0 )  & 
     612            &            + zalpha * exp( reac_nh4 * zadsnh4 * zalpha * dtsed2 /2. ) + rtrn ) 
     613            zsedtra(1) = zalpha + 2.0 * zsedtra(3) 
     614            pwcp(ji,jk,jwalk) = zgamma + 2.0 * zsedtra(3) 
     615            ! FeS - O2 
     616            zalpha = zsedtra(1) - 2.0 * zsedtra(6) 
     617            zbeta  = zsedtra(4) + zsedtra(6) 
     618            zgamma = pwcp(ji,jk,jwso4) + zsedtra(6) 
     619            zsedtra(6) = ( zsedtra(6) * zalpha ) / ( 2.0 * zsedtra(6) * ( exp( reac_feso * zalpha * dtsed2 / 2. ) - 1.0 )  & 
     620            &            + zalpha * exp( reac_feso * zalpha * dtsed2 /2. ) + rtrn ) 
     621            zsedtra(1) = zalpha + 2.0 * zsedtra(6) 
     622            zsedtra(4) = zbeta  - zsedtra(6) 
     623            pwcp(ji,jk,jwso4) = zgamma - zsedtra(6) 
     624!            ! Fe - H2S 
     625            zalpha = zsedtra(2) - zsedtra(4) 
     626            zbeta  = zsedtra(4) + zsedtra(6) 
     627            zgamma = pwcp(ji,jk,jwalk) - 2.0 * zsedtra(4) 
     628            zsedtra(4) = ( zsedtra(4) * zalpha ) / ( zsedtra(4) * ( exp( reac_fes * zalpha * dtsed2 / 2. ) - 1.0 )  & 
     629            &            + zalpha * exp( reac_fes * zalpha * dtsed2 /2. ) + rtrn ) 
     630            zsedtra(2) = zalpha + zsedtra(4) 
     631            zsedtra(6) = zbeta  - zsedtra(4) 
     632            pwcp(ji,jk,jwalk) = zgamma + 2.0 * zsedtra(4) 
     633            ! FEOH + H2S 
     634            zalpha = zsedtra(5) - 2.0 * zsedtra(2) 
     635            zbeta  = zsedtra(5) + zsedtra(4) 
     636            zgamma = pwcp(ji,jk,jwalk) - 2.0 * zsedtra(4) 
     637            zdelta = pwcp(ji,jk,jwso4) + zsedtra(2) 
     638            zepsi  = pwcp(ji,jk,jwpo4) + redfep * zsedtra(5) 
     639            zsedtra(2) = ( zsedtra(2) * zalpha ) / ( 2.0 * zsedtra(2) * ( exp( reac_feh2s * zalpha * dtsed2 ) - 1.0 )  & 
     640            &            + zalpha * exp( reac_feh2s * zalpha * dtsed2 ) + rtrn ) 
     641            zsedtra(5) = zalpha + 2.0 * zsedtra(2) 
     642            zsedtra(4) = zbeta  - zsedtra(5) 
     643            pwcp(ji,jk,jwso4) = zdelta - zsedtra(2) 
     644            pwcp(ji,jk,jwalk) = zgamma + 2.0 * zsedtra(4) 
     645            pwcp(ji,jk,jwpo4) = zepsi - redfep * zsedtra(5) 
     646            ! Fe - H2S 
     647            zalpha = zsedtra(2) - zsedtra(4) 
     648            zbeta  = zsedtra(4) + zsedtra(6) 
     649            zgamma = pwcp(ji,jk,jwalk) - 2.0 * zsedtra(4) 
     650            zsedtra(4) = ( zsedtra(4) * zalpha ) / ( zsedtra(4) * ( exp( reac_fes * zalpha * dtsed2 / 2. ) - 1.0 )  & 
     651            &            + zalpha * exp( reac_fes * zalpha * dtsed2 /2. ) + rtrn ) 
     652            zsedtra(2) = zalpha + zsedtra(4) 
     653            zsedtra(6) = zbeta  - zsedtra(4) 
     654            pwcp(ji,jk,jwalk) = zgamma + 2.0 * zsedtra(4) 
     655            ! FeS - O2 
     656            zalpha = zsedtra(1) - 2.0 * zsedtra(6) 
     657            zbeta  = zsedtra(4) + zsedtra(6) 
     658            zgamma = pwcp(ji,jk,jwso4) + zsedtra(6) 
     659            zsedtra(6) = ( zsedtra(6) * zalpha ) / ( 2.0 * zsedtra(6) * ( exp( reac_feso * zalpha * dtsed2 / 2. ) - 1.0 )  & 
     660            &            + zalpha * exp( reac_feso * zalpha * dtsed2 /2. ) + rtrn ) 
     661            zsedtra(1) = zalpha + 2.0 * zsedtra(6) 
     662            zsedtra(4) = zbeta  - zsedtra(6) 
     663            pwcp(ji,jk,jwso4) = zgamma - zsedtra(6) 
     664            ! NH4 + O2 
     665            zalpha = zsedtra(1) - 2.0 * zsedtra(3) 
     666            zgamma = pwcp(ji,jk,jwalk) - 2.0 * zsedtra(3) 
     667            zsedtra(3) = ( zsedtra(3) * zalpha ) / ( 2.0 * zsedtra(3) * ( exp( reac_nh4 * zadsnh4 * zalpha * dtsed2 / 2. ) - 1.0 )  & 
     668            &            + zalpha * exp( reac_nh4 * zadsnh4 * zalpha * dtsed2 /2. ) + rtrn ) 
     669            zsedtra(1) = zalpha + 2.0 * zsedtra(3) 
     670            pwcp(ji,jk,jwalk) = zgamma + 2.0 * zsedtra(3) 
     671            ! H2S + O2 
     672            zalpha = zsedtra(1) - 2.0 * zsedtra(2) 
     673            zbeta  = pwcp(ji,jk,jwso4) + zsedtra(2) 
     674            zgamma = pwcp(ji,jk,jwalk) - 2.0 * zsedtra(2) 
     675            zsedtra(2) = ( zsedtra(2) * zalpha ) / ( 2.0 * zsedtra(2) * ( exp( reac_h2s * zalpha * dtsed2 / 2. ) - 1.0 )  & 
     676            &            + zalpha * exp( reac_h2s * zalpha * dtsed2 / 2. ) + rtrn ) 
     677            zsedtra(1) = zalpha + 2.0 * zsedtra(2) 
     678            pwcp(ji,jk,jwso4) = zbeta - zsedtra(2) 
     679            pwcp(ji,jk,jwalk) = zgamma + 2.0 * zsedtra(2) 
     680            ! Fe + O2 
     681            zalpha = zsedtra(1) - 0.25 * zsedtra(4) 
     682            zbeta  = zsedtra(4) + zsedtra(5) 
     683            zgamma = pwcp(ji,jk,jwalk) - 2.0 * zsedtra(4) 
     684            zdelta = pwcp(ji,jk,jwpo4) - redfep * zsedtra(4) 
     685            zsedtra(4) = ( zsedtra(4) * zalpha ) / ( 0.25 * zsedtra(4) * ( exp( reac_fe2 * zalpha * dtsed2 / 2. ) - 1.0 )  & 
     686            &            + zalpha * exp( reac_fe2 * zalpha * dtsed2 / 2. ) + rtrn ) 
     687            zsedtra(1) = zalpha + 0.25 * zsedtra(4) 
     688            zsedtra(5) = zbeta  - zsedtra(4) 
     689            pwcp(ji,jk,jwpo4) = zdelta + redfep * zsedtra(4) 
     690            pwcp(ji,jk,jwalk) = zgamma + 2.0 * zsedtra(4) 
     691            pwcp(ji,jk,jwoxy)  = zsedtra(1) 
     692            pwcp(ji,jk,jwh2s)  = zsedtra(2) 
     693            pwcp(ji,jk,jwnh4)  = zsedtra(3) 
     694            pwcp(ji,jk,jwfe2)  = zsedtra(4) + sedligand(ji,jk) + zsedfer 
     695            pwcp(ji,jk,jwno3)  = pwcp(ji,jk,jwno3)  + ( zsedtrn(3) - pwcp(ji,jk,jwnh4) ) 
     696            solcp(ji,jk,jsfeo) = zsedtra(5) / zvolc(ji,jk,jsfeo) 
     697            solcp(ji,jk,jsfes) = zsedtra(6) / zvolc(ji,jk,jsfes) 
     698         END DO 
     699      END DO 
     700 
     701      IF( ln_timing )  CALL timing_stop('sed_dsr_redoxb') 
     702 
     703  END SUBROUTINE sed_dsr_redoxb 
     704 
    541705END MODULE seddsr 
Note: See TracChangeset for help on using the changeset viewer.