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 14323 – NEMO

Changeset 14323


Ignore:
Timestamp:
2021-01-21T00:33:21+01:00 (3 years ago)
Author:
aumont
Message:

major update of the sediment model

Location:
NEMO/branches/2019/dev_r11708_aumont_PISCES_QUOTA
Files:
1 added
1 deleted
12 edited

Legend:

Unmodified
Added
Removed

  • NEMO/branches/2019/dev_r11708_aumont_PISCES_QUOTA/cfgs/SHARED/field_def_nemo-pisces.xml

    r14276 r14323  
    137137       <field id="SedpH"           long_name="PH"                                       unit="1"          /> 
    138138       <field id="SedCO3por"       long_name="Bicarbonates"                             unit="mol/m3" /> 
    139        <field id="Sedligand"       long_name="Ligands Concentration"                    unit="mol/m3" /> 
     139       <field id="Sedligand"       long_name="Ligands"                                  unit="mol/m3" /> 
     140       <field id="SaturCO3"        long_name="CO3 Saturation"                           unit="-" /> 
    140141     </field_group> 
    141142 
    … …  
    153154       <field id="FlxFe2"      long_name="Fe2+ sediment flux"                      unit="mol/cm2/s"     /> 
    154155       <field id="dzdep"       long_name="Sedimentation rate"                      unit="cm/s"     /> 
    155        <field id="sflxclay"    long_name="Clay sedimentation rate"                 unit="g/m2/s"     /> 
    156        <field id="sflxcal"     long_name="Calcite sedimentation rate"              unit="mol/m2/s"     /> 
    157        <field id="sflxbsi"     long_name="BSi Sedimentation rate"                  unit="mol/m2/s"     /> 
    158        <field id="sflxpoc"     long_name="POC Sedimentation rate"                  unit="mol/m2/s"     /> 
    159        <field id="sflxfeo"     long_name="Fe(OH)3 Sedimentation rate"              unit="mol/m2/s"     /> 
    160        <field id="FlxBSi"      long_name="BSi burial flux"                         unit="g/cm2/s"     /> 
    161        <field id="FlxClay"     long_name="Clay burial flux"                        unit="g/cm2/s"     /> 
    162        <field id="FlxPOC"      long_name="POC burial flux"                         unit="g/cm2/s"     /> 
    163        <field id="FlxCaCO3"    long_name="Calcite burial flux"                     unit="g/cm2/s"     /> 
    164        <field id="FlxPOS"      long_name="POS burial flux"                         unit="g/cm2/s"     /> 
    165        <field id="FlxPOR"      long_name="POR burial flux"                         unit="g/cm2/s"     /> 
    166        <field id="FlxFeO"      long_name="FeO burial flux"                         unit="g/cm2/s"     /> 
    167        <field id="FlxFeS"      long_name="FeS burial flux"                         unit="g/cm2/s"     /> 
     156       <field id="sflxclay"    long_name="Clay sedimentation rate"                 unit="g/cm2/s"     /> 
     157       <field id="sflxbsi"     long_name="BSi sedimentation rate"                  unit="g/cm2/s"     /> 
     158       <field id="sflxpoc"     long_name="POC sedimentation rate"                  unit="g/cm2/s"     /> 
     159       <field id="sflxcal"     long_name="Calcite sedimentation rate"              unit="g/cm2/s"     /> 
     160       <field id="FlxClay"     long_name="Clay burial rate"                        unit="g/cm2/s"     /> 
     161       <field id="FlxCaCO3"    long_name="Calcite burial rate"                     unit="g/cm2/s"     /> 
     162       <field id="FlxBSi"      long_name="BSi burial rate"                         unit="g/cm2/s"     /> 
     163       <field id="FlxPOC"      long_name="POC burial rate"                         unit="g/cm2/s"     /> 
     164       <field id="FlxFeO"      long_name="Fe(OH)3 burial rate"                     unit="g/cm2/s"     /> 
     165       <field id="FlxFeS"      long_name="FeS burial rate"                         unit="g/cm2/s"     /> 
     166       <field id="FlxPOR"      long_name="POR burial rate"                         unit="g/cm2/s"     /> 
     167       <field id="FlxPOS"      long_name="POS burial rate"                         unit="g/cm2/s"     /> 
     168       <field id="Flxtot"      long_name="total burial flux"                       unit="g/cm2/s"     /> 
     169 
    168170</field_group> 
    169171 
    … …  
    215217       <field id="SIZEP"       long_name="Mean relative size of picophyto."        unit="-"          grid_ref="grid_T_3D" /> 
    216218       <field id="SIZED"       long_name="Mean relative size of diatoms"           unit="-"          grid_ref="grid_T_3D" /> 
    217        <field id="RASSD"       long_name="Mean relative size of diatoms"           unit="-"          grid_ref="grid_T_3D" /> 
    218        <field id="RASSN"       long_name="Mean relative size of diatoms"           unit="-"          grid_ref="grid_T_3D" /> 
    219        <field id="RASSP"       long_name="Mean relative size of diatoms"           unit="-"          grid_ref="grid_T_3D" /> 
    220219       <field id="Fe3"         long_name="Iron III concentration"                  unit="nmol/m3"    grid_ref="grid_T_3D" /> 
    221220       <field id="FeL1"        long_name="Complexed Iron concentration with L1"    unit="nmol/m3"    grid_ref="grid_T_3D" /> 

  • NEMO/branches/2019/dev_r11708_aumont_PISCES_QUOTA/cfgs/SHARED/namelist_sediment_ref

    r14306 r14323  
    1212&nam_run      !  Characteristics of the simulation 
    1313!----------------------------------------------------------------------- 
    14   nrseddt   = 1          ! Nb of iterations for fast species 
     14  nrseddt   = 2          ! Nb of iterations for fast species 
    1515  ln_sed_2way = .false.  ! 2 way coupling with pisces 
    1616/ 
    … …  
    5858   seddiag3d(1)   = 'SedpH      ' , 'pH                                     ',  '-      ' 
    5959   seddiag3d(2)   = 'SedCO3por  ' , 'Dissolved CO3 concentration            ',  'mol/L' 
    60    seddiag3d(3)   = 'Sedligand'   , 'ligand concentration                   ',  'mol/L' 
     60   seddiag3d(3)   = 'Sedligand  ' , 'ligand concentration                   ',  'mol/L' 
     61   seddiag3d(4)   = 'SaturCO3   ' , 'CO3 saturation                         ',  '-' 
    6162   seddiag2d(1)   = 'FlxSi      ' , 'Silicate flux                          ',  'mol/cm2/s' 
    6263   seddiag2d(2)   = 'FlxO2      ' , 'Dissolved Oxygen Flux                  ',  'mol/L' 
    … …  
    99100   rcorgl   =  10.     ! Reactivity for labile POC [an-1] 
    100101   rcorgs   =  0.1     ! Reactivity for semi-refractory POC [an-1] 
    101    rcorgr   =  1.E-4   ! Reactivity for refractory POC [an-1] 
     102   rcorgr   =  5.E-4   ! Reactivity for refractory POC [an-1] 
    102103   rcnh4    =  1E7     ! Reactivity for O2/NH4 [l.mol-1.an-1] 
    103104   rch2s    =  2E8     ! Reactivity for O2/H2S [l.mol-1.an-1] 
    … …  
    106107   rcfes    =  1E6     ! Reactivity for FE2+/H2S [l.mol-1.an-1] 
    107108   rcfeso   =  2E7     ! Reactivity for FES/O2 [l.mol-1.an-1] 
    108    xksedo2  =  1.E-6   ! Half-saturation constant for oxic remin [mol/l] 
     109   xksedo2  =  4.E-6   ! Half-saturation constant for oxic remin [mol/l] 
    109110   xksedno3 =  10.E-6  ! Half-saturation constant for denitrification [mol/l]  
    110111   xksedfeo =  0.007   ! Half-saturation constant for iron remin [%] 
    … …  
    130131   cn_sedrst_outdir = "."          !  directory to which to write output sediment restarts 
    131132/ 
    132 !----------------------------------------------------------------------- 
    133 &nam_output     !   parameters for outputing the sediment module 
    134 !----------------------------------------------------------------------- 
    135    ldefsedpis_avg = .true.            !  write averaged output variables 
    136    cn_sedwri_out  = "output_sed.nc"   !  name of the input restart file name of the sediment module 
    137    nrpfsedpis_avg = 0                 ! Frequency of the averaged outputs 
    138    nwrtsedpis_avg = 24                ! Frequency of the averaged outputs 
    139    ntssedpis_avg  =  1                ! ??? 
    140 / 

  • NEMO/branches/2019/dev_r11708_aumont_PISCES_QUOTA/src/TOP/PISCES/SED/par_sed.F90

    r14276 r14323  
    6161   INTEGER, PARAMETER ::  & 
    6262      jptrased   = jpsol + jpwat , & 
    63       jpdia3dsed = 3             , & 
     63      jpdia3dsed = 4             , & 
    6464      jpdia2dsed = 20  
    6565 

  • NEMO/branches/2019/dev_r11708_aumont_PISCES_QUOTA/src/TOP/PISCES/SED/sed.F90

    r10425 r14323  
    7373   REAL(wp), PUBLIC, DIMENSION(:,:  ), ALLOCATABLE ::  fromsed    !: 
    7474   REAL(wp), PUBLIC, DIMENSION(:,:  ), ALLOCATABLE ::  tosed      !: 
    75    REAL(wp), PUBLIC, DIMENSION(:,:  ), ALLOCATABLE ::  rloss      !:  
    76    REAL(wp), PUBLIC, DIMENSION(:,:  ), ALLOCATABLE ::  tokbot         
    77    ! 
    7875   REAL(wp), PUBLIC, DIMENSION(:    ), ALLOCATABLE ::  temp       !: temperature 
    7976   REAL(wp), PUBLIC, DIMENSION(:    ), ALLOCATABLE ::  salt       !: salinity 
    80    REAL(wp), PUBLIC, DIMENSION(:    ), ALLOCATABLE ::  press      !: pressure 
    8177   REAL(wp), PUBLIC, DIMENSION(:    ), ALLOCATABLE ::  raintg     !: total massic flux rained in each cell (sum of sol. comp.) 
    8278   REAL(wp), PUBLIC, DIMENSION(:    ), ALLOCATABLE ::  fecratio   !: Fe/C ratio in falling particles to the sediments 
    8379   REAL(wp), PUBLIC, DIMENSION(:    ), ALLOCATABLE ::  dzdep      !: total thickness of solid material rained [cm] in each cell 
    84    REAL(wp), PUBLIC, DIMENSION(:    ), ALLOCATABLE ::  zkbot      !: total thickness of solid material rained [cm] in each cell 
     80   REAL(wp), PUBLIC, DIMENSION(:    ), ALLOCATABLE ::  zkbot      !: Total water column depth 
    8581   REAL(wp), PUBLIC, DIMENSION(:    ), ALLOCATABLE ::  wacc       !: total thickness of solid material rained [cm] in each cell 
    8682   ! 
    … …  
    9187   REAL(wp), PUBLIC, DIMENSION(:,:  ), ALLOCATABLE ::  vols3d     !:  ??? 
    9288 
    93  
    9489   !! Chemistry 
    95    REAL(wp), PUBLIC, DIMENSION(:    ), ALLOCATABLE ::  densSW  
    96    REAL(wp), PUBLIC, DIMENSION(:    ), ALLOCATABLE ::  borats  
     90   REAL(wp), PUBLIC, DIMENSION(:    ), ALLOCATABLE ::  densSW 
     91   REAL(wp), PUBLIC, DIMENSION(:    ), ALLOCATABLE ::  borats 
    9792   REAL(wp), PUBLIC, DIMENSION(:    ), ALLOCATABLE ::  calcon2 
    98    REAL(wp), PUBLIC, DIMENSION(:    ), ALLOCATABLE ::  akbs   
    99    REAL(wp), PUBLIC, DIMENSION(:    ), ALLOCATABLE ::  ak1s  
    100    REAL(wp), PUBLIC, DIMENSION(:    ), ALLOCATABLE ::  ak2s    
    101    REAL(wp), PUBLIC, DIMENSION(:    ), ALLOCATABLE ::  akws   
    102    REAL(wp), PUBLIC, DIMENSION(:    ), ALLOCATABLE ::  ak12s   
    103    REAL(wp), PUBLIC, DIMENSION(:    ), ALLOCATABLE ::  ak1ps  
    104    REAL(wp), PUBLIC, DIMENSION(:    ), ALLOCATABLE ::  ak2ps   
    105    REAL(wp), PUBLIC, DIMENSION(:    ), ALLOCATABLE ::  ak3ps  
    106    REAL(wp), PUBLIC, DIMENSION(:    ), ALLOCATABLE ::  ak12ps  
     93   REAL(wp), PUBLIC, DIMENSION(:    ), ALLOCATABLE ::  akbs 
     94   REAL(wp), PUBLIC, DIMENSION(:    ), ALLOCATABLE ::  ak1s 
     95   REAL(wp), PUBLIC, DIMENSION(:    ), ALLOCATABLE ::  ak2s 
     96   REAL(wp), PUBLIC, DIMENSION(:    ), ALLOCATABLE ::  akws 
     97   REAL(wp), PUBLIC, DIMENSION(:    ), ALLOCATABLE ::  ak12s 
     98   REAL(wp), PUBLIC, DIMENSION(:    ), ALLOCATABLE ::  ak1ps 
     99   REAL(wp), PUBLIC, DIMENSION(:    ), ALLOCATABLE ::  ak2ps 
     100   REAL(wp), PUBLIC, DIMENSION(:    ), ALLOCATABLE ::  ak3ps 
     101   REAL(wp), PUBLIC, DIMENSION(:    ), ALLOCATABLE ::  ak12ps 
    107102   REAL(wp), PUBLIC, DIMENSION(:    ), ALLOCATABLE ::  ak123ps 
    108    REAL(wp), PUBLIC, DIMENSION(:    ), ALLOCATABLE ::  aksis  
    109    REAL(wp), PUBLIC, DIMENSION(:    ), ALLOCATABLE ::  aksps  
     103   REAL(wp), PUBLIC, DIMENSION(:    ), ALLOCATABLE ::  aksis 
     104   REAL(wp), PUBLIC, DIMENSION(:    ), ALLOCATABLE ::  aksps 
    110105   REAL(wp), PUBLIC, DIMENSION(:    ), ALLOCATABLE ::  sieqs 
    111106   REAL(wp), PUBLIC, DIMENSION(:    ), ALLOCATABLE ::  aks3s 
    … …  
    130125   REAL(wp), PUBLIC, DIMENSION(:,:  ), ALLOCATABLE ::  db            !: bioturbation ceofficient 
    131126   REAL(wp), PUBLIC, DIMENSION(:,:  ), ALLOCATABLE ::  irrig        !: bioturbation ceofficient 
    132    REAL(wp), PUBLIC, DIMENSION(:    ), ALLOCATABLE ::  rdtsed        !:  sediment model time-step 
    133127   REAL(wp), PUBLIC, DIMENSION(:,:  ), ALLOCATABLE :: sedligand 
     128   REAL(wp), PUBLIC, DIMENSION(:,:  ), ALLOCATABLE :: saturco3 
    134129   REAL(wp)  ::   dens               !: density of solid material 
    135130   !! Inputs / Outputs 
    … …  
    157152      !!------------------------------------------------------------------- 
    158153      ! 
    159       ALLOCATE( trc_data(jpi,jpj,jpdta)                                   ,   & 
    160          &      epkbot(jpi,jpj), gdepbot(jpi,jpj)        ,   & 
    161          &      dz(jpksed)  , por(jpksed) , por1(jpksed)                  ,   & 
    162          &      volw(jpksed), vols(jpksed), rdtsed(jpksed)  ,   & 
    163          &      trcsedi  (jpi,jpj,jpksed,jptrased)                        ,   & 
    164          &      flxsedi3d(jpi,jpj,jpksed,jpdia3dsed)                      ,   & 
    165          &      flxsedi2d(jpi,jpj,jpdia2dsed)                             ,   & 
    166          &      mol_wgt(jpsol),                                           STAT=sed_alloc ) 
     154      ALLOCATE( trc_data(jpi,jpj,jpdta), trcsedi  (jpi,jpj,jpksed,jptrased) ,  & 
     155         &      epkbot(jpi,jpj), gdepbot(jpi,jpj), dz(jpksed), por(jpksed)  ,  & 
     156         &      por1(jpksed), volw(jpksed), vols(jpksed), mol_wgt(jpsol)    ,  & 
     157         &      flxsedi2d(jpi,jpj,jpdia2dsed), flxsedi3d(jpi,jpj,jpksed,jpdia3dsed),  STAT=sed_alloc ) 
    167158 
    168159      IF( sed_alloc /= 0 )   CALL ctl_stop( 'STOP', 'sed_alloc: failed to allocate arrays' ) 

  • NEMO/branches/2019/dev_r11708_aumont_PISCES_QUOTA/src/TOP/PISCES/SED/sedadv.F90

    r10425 r14323  
    8888      fromsed(:,:) = 0. 
    8989      tosed  (:,:) = 0.  
    90       rloss  (:,:) = 0. 
    9190      ikwneg = 1 
    9291      nztime = jpksed 
    … …  
    254253 
    255254         ELSE IF( raintg(ji) < eps ) THEN ! rain  = 0 
    256 !! Nadia    rloss(:,:) = rainrm(:,:)   bug ??????            
    257  
    258             rloss(ji,1:jpsol) = rainrm(ji,1:jpsol) 
    259255 
    260256            zfull (2) = zfilled(2) 
    … …  
    300296            IF( ikwneg == 2 ) THEN ! advection is reversed in the first sediment layer 
    301297 
    302                zkwnup = rdtsed(ikwneg) * raintg(ji) / dz(ikwneg) 
     298               zkwnup = dtsed * raintg(ji) / ( denssol * por1(ikwneg) * dz(ikwneg) ) 
    303299               zkwnlo = ABS( zwb(ji,ikwneg) ) / dz(ikwneg) 
    304300               zfull (ikwneg+1) = zfilled(ikwneg+1) - zkwnlo * dvolsm(ikwneg+1) 
    … …  
    347343               ! Heinze  fromsed(ji,jsclay) = zkwnlo * 1. * denssol * por1(jpksed) / mol_wgt(jsclay) 
    348344               fromsed(ji,jsclay) = zkwnlo * 1.* por1clay 
    349             ELSE   ! 2 < ikwneg(ji) <= jpksedm1 
     345            ELSE IF( ikwneg > 2 .AND. ikwneg < jpksed-1 ) THEN 
    350346 
    351347               zkwnup = ABS( zwb(ji,ikwneg-1) ) * por1(ikwneg-1) / ( dz(ikwneg) * por1(ikwneg) ) 

  • NEMO/branches/2019/dev_r11708_aumont_PISCES_QUOTA/src/TOP/PISCES/SED/seddiff.F90

    r10225 r14323  
    4343      INTEGER, INTENT(in) ::   kt, knt       ! number of iteration 
    4444      ! --- local variables 
    45       INTEGER :: ji, jk, js   ! dummy looop indices 
     45      INTEGER :: ji, jk, jw   ! dummy looop indices 
    4646 
    4747      REAL(wp), DIMENSION(jpoce,jpksed) :: zrearat1, zrearat2   ! reaction rate in pore water 
    … …  
    6060     ! Initializations 
    6161     !---------------------- 
    62       zrearat1(:,:)   = 0. 
     62      zrearat1(:,:) = 0. 
    6363      zrearat2(:,:) = 0. 
    6464 
    65       !--------------------------- 
    66       ! Solves PO4 diffusion  
    67       !---------------------------- 
     65      ! -------------------------------------- 
     66      ! Solves solute diffusion in pore waters 
     67      ! -------------------------------------- 
    6868 
    69       ! solves tridiagonal system 
    70       CALL sed_mat( jwpo4, jpoce, jpksed, zrearat1, zrearat2, pwcp(:,:,jwpo4), dtsed2 / 2.0 ) 
    71  
    72       !--------------------------- 
    73       ! Solves NH4 diffusion  
    74       !---------------------------- 
    75  
    76       ! solves tridiagonal system 
    77       CALL sed_mat( jwnh4, jpoce, jpksed, zrearat1, zrearat2, pwcp(:,:,jwnh4), dtsed2 / 2.0 ) 
    78  
    79       !--------------------------- 
    80       ! Solves Fe2+ diffusion  
    81       !---------------------------- 
    82  
    83       ! solves tridiagonal system 
    84       CALL sed_mat( jwfe2, jpoce, jpksed, zrearat1, zrearat2, pwcp(:,:,jwfe2), dtsed2 / 2.0 ) 
    85  
    86       !--------------------------- 
    87       ! Solves H2S diffusion  
    88       !---------------------------- 
    89  
    90       ! solves tridiagonal system 
    91       CALL sed_mat( jwh2s, jpoce, jpksed, zrearat1, zrearat2, pwcp(:,:,jwh2s), dtsed2 / 2.0  ) 
    92  
    93       !--------------------------- 
    94       ! Solves SO4 diffusion  
    95       !---------------------------- 
    96  
    97       ! solves tridiagonal system 
    98       CALL sed_mat( jwso4, jpoce, jpksed, zrearat1, zrearat2, pwcp(:,:,jwso4), dtsed2 / 2.0 ) 
    99  
    100       !--------------------------- 
    101       ! Solves O2 diffusion 
    102       !---------------------------- 
    103  
    104       ! solves tridiagonal system 
    105       CALL sed_mat( jwoxy, jpoce, jpksed, zrearat1, zrearat2, pwcp(:,:,jwoxy), dtsed2 / 2.0 ) 
    106  
    107       !--------------------------- 
    108       ! Solves NO3 diffusion 
    109       !---------------------------- 
    110  
    111       ! solves tridiagonal system 
    112       CALL sed_mat( jwno3, jpoce, jpksed, zrearat1, zrearat2, pwcp(:,:,jwno3), dtsed2 / 2.0 ) 
     69      ! Solves tridiagonal system 
     70      DO jw = 1, jpwat 
     71         CALL sed_mat( jw, jpoce, jpksed, zrearat1, zrearat2, pwcp(:,:,jw), dtsed2 / 2.0 ) 
     72      END DO 
    11373 
    11474      CALL sed_mat( jwdic, jpoce, jpksed, zrearat1, zrearat2, sedligand(:,:), dtsed2 / 2.0 ) 

  • NEMO/branches/2019/dev_r11708_aumont_PISCES_QUOTA/src/TOP/PISCES/SED/seddsr.F90

    r10362 r14323  
    5555      REAL(wp), DIMENSION(jpoce,jpksed) :: zrearat1, zrearat2, zrearat3    ! reaction rate in pore water 
    5656      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    
     57      REAL(wp), DIMENSION(jpoce,jpksed) :: zlimo2, zlimno3, zlimso4, zlimfeo    ! undersaturation ; indice jpwatp1 is for calcite    
    5858      REAL(wp), DIMENSION(jpoce)        :: zsumtot 
    5959      REAL(wp)  ::  zsolid1, zsolid2, zsolid3, zvolw, zreasat 
    60       REAL(wp)  ::  zsatur, zsatur2, znusil, zkpoca, zkpocb, zkpocc 
    61       REAL(wp)  ::  zratio, zgamma, zbeta, zlimtmp, zundsat2 
     60      REAL(wp)  ::  zgamma, zbeta, zlimtmp 
    6261      !! 
    6362      !!---------------------------------------------------------------------- 
    … …  
    7574     !---------------------- 
    7675       
    77       zrearat1(:,:)   = 0.    ;   zundsat(:,:) = 0. ; zkpoc(:,:) = 0. 
    78       zlimo2 (:,:)    = 0.    ;   zlimno3(:,:) = 0. ; zrearat2(:,:) = 0. 
    79       zlimso4(:,:)    = 0.    ;   zkpor(:,:)   = 0. ; zrearat3(:,:) = 0. 
    80       zkpos  (:,:)    = 0. 
     76      zrearat1(:,:)   = 0.    ;   zundsat(:,:)  = 0. 
     77      zlimo2 (:,:)    = 0.    ;   zlimno3(:,:)  = 0. ; zrearat2(:,:) = 0. 
     78      zlimso4(:,:)    = 0.    ;   zrearat3(:,:) = 0. 
    8179      zsumtot(:)      = rtrn 
    8280   
    … …  
    8482      zvolc(:,:,:)    = 0. 
    8583      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 
    9684 
    9785      ! Conversion of volume units 
    … …  
    116104      zrearat3(:,:) = 0. 
    117105 
    118       zundsat(:,:) = pwcp(:,:,jwoxy) 
     106      zundsat(:,:) = MAX( pwcp(:,:,jwoxy) - rtrn, 0. ) 
    119107 
    120108      DO jk = 2, jpksed 
    121109         DO ji = 1, jpoce 
    122             zlimo2(ji,jk) = 1.0 / ( zundsat(ji,jk) + xksedo2 ) 
    123110            zsolid1 = zvolc(ji,jk,jspoc)  * solcp(ji,jk,jspoc) 
    124111            zsolid2 = zvolc(ji,jk,jspos)  * solcp(ji,jk,jspos) 
    125112            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 
    137  
    138       ! left hand side of coefficient matrix 
    139 !      DO jn = 1, 5 
    140       DO jk = 2, jpksed 
    141          DO ji = 1, jpoce 
    142 jflag1:     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 
    167  
    168       ! New solid concentration values (jk=2 to jksed) for each couple  
    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 
    177          ENDDO 
    178       ENDDO 
    179  
    180       ! New pore water concentrations     
    181       DO jk = 2, jpksed 
    182          DO ji = 1, jpoce 
    183             ! Acid Silicic  
    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          
     113            zrearat1(ji,jk)  = ( reac_pocl * dtsed2 * zsolid1 ) / ( 1. + reac_pocl * dtsed2 ) 
     114            zrearat2(ji,jk)  = ( reac_pocs * dtsed2 * zsolid2 ) / ( 1. + reac_pocs * dtsed2 ) 
     115            zrearat3(ji,jk)  = ( reac_pocr * dtsed2 * zsolid3 ) / ( 1. + reac_pocr * dtsed2 ) 
     116            solcp(ji,jk,jspoc) = solcp(ji,jk,jspoc) - zrearat1(ji,jk) / zvolc(ji,jk,jspoc) 
     117            solcp(ji,jk,jspos) = solcp(ji,jk,jspos) - zrearat2(ji,jk) / zvolc(ji,jk,jspos) 
     118            solcp(ji,jk,jspor) = solcp(ji,jk,jspor) - zrearat3(ji,jk) / zvolc(ji,jk,jspor) 
     119            zreasat = zrearat1(ji,jk) + zrearat2(ji,jk) + zrearat3(ji,jk) 
     120 
    186121            ! For DIC 
    187122            pwcp(ji,jk,jwdic)  = pwcp(ji,jk,jwdic) + zreasat 
    188123            zsumtot(ji) = zsumtot(ji) + zreasat / dtsed2 * volw3d(ji,jk) * 1.e-3 * 86400. * 365. * 1E3 
     124 
     125            ! For alkalinity 
     126            pwcp(ji,jk,jwalk)  = pwcp(ji,jk,jwalk) + zreasat * ( srno3 - 2.* spo4r ) 
     127 
    189128            ! For Phosphate (in mol/l) 
    190129            pwcp(ji,jk,jwpo4)  = pwcp(ji,jk,jwpo4) + zreasat * spo4r 
     130 
    191131            ! For iron (in mol/l) 
    192132            pwcp(ji,jk,jwfe2)  = pwcp(ji,jk,jwfe2) + fecratio(ji) * zreasat 
    193             ! For alkalinity 
    194             pwcp(ji,jk,jwalk)  = pwcp(ji,jk,jwalk) + zreasat * ( srno3 * zadsnh4 - 2.* spo4r ) 
     133 
    195134            ! Ammonium 
    196135            pwcp(ji,jk,jwnh4)  = pwcp(ji,jk,jwnh4) + zreasat * srno3 * zadsnh4 
     136 
    197137            ! Ligands 
    198             sedligand(ji,jk)   = sedligand(ji,jk) + ratligc * zreasat - reac_ligc * sedligand(ji,jk) 
     138            sedligand(ji,jk)   = sedligand(ji,jk) + ratligc * zreasat 
     139 
     140         ENDDO 
     141      ENDDO 
     142 
     143      ! left hand side of coefficient matrix 
     144      DO jk = 2, jpksed 
     145         DO ji = 1, jpoce 
     146            zreasat = zrearat1(ji,jk) + zrearat2(ji,jk) + zrearat3(ji,jk) 
     147            zbeta   = xksedo2 - pwcp(ji,jk,jwoxy) + so2ut * zreasat 
     148            zgamma = - xksedo2 * pwcp(ji,jk,jwoxy) 
     149            zundsat(ji,jk) = ( - zbeta + SQRT( zbeta**2 - 4.0 * zgamma ) ) / 2.0 
     150            zlimo2(ji,jk) = zundsat(ji,jk) / ( zundsat(ji,jk) + xksedo2 ) 
     151 
     152            ! Oxygen 
     153            pwcp(ji,jk,jwoxy)  = zundsat(ji,jk) + MIN( pwcp(ji,jk,jwoxy), rtrn ) 
     154            zreasat = zreasat * zlimo2(ji,jk) 
     155 
     156            ! Ligands 
     157            sedligand(ji,jk)   = sedligand(ji,jk) - reac_ligc * sedligand(ji,jk) 
    199158         ENDDO 
    200159      ENDDO 
    … …  
    205164      !-------------------------------------------------------------------- 
    206165 
    207       zrearat1(:,:) = 0. 
    208       zrearat2(:,:) = 0. 
    209       zrearat3(:,:) = 0. 
    210  
    211       zundsat(:,:) = pwcp(:,:,jwno3) 
     166      zundsat(:,:) = MAX(0., pwcp(:,:,jwno3) - rtrn) 
    212167 
    213168      DO jk = 2, jpksed 
    214169         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 
    234 jflag2:    DO jn = 1, 10 
    235                zlimtmp = ( 1.0 - pwcp(ji,jk,jwoxy) * zlimo2(ji,jk) ) 
    236                zsolid1 = zvolc(ji,jk,jspoc) * solcp(ji,jk,jspoc) 
    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 
    260  
    261  
    262       ! New solid concentration values (jk=2 to jksed) for each couple  
    263       DO jk = 2, jpksed 
    264          DO ji = 1, jpoce 
    265             zreasat = zrearat1(ji,jk) * zlimno3(ji,jk) * zundsat(ji,jk) / zvolc(ji,jk,jspoc) 
    266             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 
    271          ENDDO 
    272       ENDDO 
    273  
    274       ! New dissolved concentrations 
    275       DO jk = 2, jpksed 
    276          DO ji = 1, jpoce 
     170            zlimtmp = ( 1.0 - zlimo2(ji,jk) ) 
     171            zreasat = zrearat1(ji,jk) + zrearat2(ji,jk) + zrearat3(ji,jk) 
     172            zbeta   = xksedno3 - pwcp(ji,jk,jwno3) + srDnit * zreasat * zlimtmp 
     173            zgamma  = - xksedno3 * pwcp(ji,jk,jwno3) 
     174            zundsat(ji,jk) = ( - zbeta + SQRT( zbeta**2 - 4.0 * zgamma ) ) / 2.0 
     175            zlimno3(ji,jk) = zlimtmp * zundsat(ji,jk) / ( zundsat(ji,jk) + xksedno3 ) 
     176 
    277177            ! For nitrates 
    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) 
    280             ! For DIC 
    281             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 
    283             ! For Phosphate (in mol/l) 
    284             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 
     178            pwcp(ji,jk,jwno3) = zundsat(ji,jk) + MIN(pwcp(ji,jk,jwno3), rtrn) 
     179            zreasat = zreasat * zlimno3(ji,jk) 
     180 
    289181            ! For alkalinity 
    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 
     182            pwcp(ji,jk,jwalk) = pwcp(ji,jk,jwalk) + zreasat * srDnit 
    293183         ENDDO 
    294184      ENDDO 
    … …  
    299189      !-------------------------------------------------------------------- 
    300190 
    301       zrearat1(:,:) = 0. 
    302       zrearat2(:,:) = 0. 
    303       zrearat3(:,:) = 0. 
    304  
    305       zundsat(:,:) = solcp(:,:,jsfeo) 
     191      zundsat(:,:) = MAX(0., solcp(:,:,jsfeo) - rtrn) 
    306192 
    307193      DO jk = 2, jpksed 
    308194         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 
    329 jflag3:     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 
    370  
    371       ! New dissolved concentrations 
    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) 
     195            zlimtmp = 1.0 - zlimo2(ji,jk) - zlimno3(ji,jk) 
     196            zreasat = ( zrearat1(ji,jk) + zrearat2(ji,jk) + zrearat3(ji,jk) ) / zvolc(ji,jk,jsfeo) 
     197            zbeta   = xksedfeo - solcp(ji,jk,jsfeo) + 4.0 * zreasat * zlimtmp 
     198            zgamma  = -xksedfeo * solcp(ji,jk,jsfeo) 
     199            zundsat(ji,jk) = ( - zbeta + SQRT( zbeta**2 - 4.0 * zgamma ) ) / 2.0 
     200            zlimfeo(ji,jk) = zlimtmp * zundsat(ji,jk) / ( zundsat(ji,jk) + xksedfeo ) 
     201            zreasat = ( zrearat1(ji,jk) + zrearat2(ji,jk) + zrearat3(ji,jk) ) * zlimfeo(ji,jk) 
     202 
    375203            ! For FEOH 
    376             solcp(ji,jk,jsfeo) = zundsat(ji,jk) 
    377             ! For DIC 
    378             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 
     204            solcp(ji,jk,jsfeo) = zundsat(ji,jk) + MIN(solcp(ji,jk,jsfeo), rtrn) 
     205 
    380206            ! 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 
     207            pwcp(ji,jk,jwpo4)  = pwcp(ji,jk,jwpo4) + zreasat * 4.0 * redfep 
     208 
    386209            ! For alkalinity 
    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 
     210            pwcp(ji,jk,jwalk)  = pwcp(ji,jk,jwalk) + 8.0 * zreasat 
     211 
     212            ! Iron 
    390213            pwcp(ji,jk,jwfe2)  = pwcp(ji,jk,jwfe2) + zreasat * 4.0 
    391214         ENDDO 
    392215      ENDDO 
    393216 
    394       !-------------------------------------------------------------------- 
    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) 
     217!      !-------------------------------------------------------------------- 
     218!      ! Begining POC denitrification and NO3- diffusion 
     219!      ! (indice n�5 for couple POC/NO3- ie solcp(:,:,jspoc)/pwcp(:,:,jwno3)) 
     220!      !-------------------------------------------------------------------- 
     221! 
     222      zundsat(:,:) = MAX(0., pwcp(:,:,jwso4) - rtrn) 
    404223 
    405224      DO jk = 2, jpksed 
    406225         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 
    428 jflag4:     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 
     226            zlimtmp = 1.0 - zlimo2(ji,jk) - zlimno3(ji,jk) - zlimfeo(ji,jk) 
     227            zreasat = zrearat1(ji,jk) + zrearat2(ji,jk) + zrearat3(ji,jk) 
     228            zbeta   = xksedso4 - pwcp(ji,jk,jwso4) + 0.5 * zreasat * zlimtmp 
     229            zgamma = - xksedso4 * pwcp(ji,jk,jwso4) 
     230            zundsat(ji,jk) = ( - zbeta + SQRT( zbeta**2 - 4.0 * zgamma ) ) / 2.0 
     231            zlimso4(ji,jk) = zlimtmp * zundsat(ji,jk) / ( zundsat(ji,jk) + xksedso4 ) 
     232 
    474233            ! 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 
     234            pwcp(ji,jk,jwso4)  =  zundsat(ji,jk) + MIN(pwcp(ji,jk,jwso4), rtrn) 
     235            zreasat = zreasat * zlimso4(ji,jk) 
     236            pwcp(ji,jk,jwh2s)  = pwcp(ji,jk,jwh2s) + 0.5 * zreasat 
     237 
    487238            ! 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 
     239            pwcp(ji,jk,jwalk)  = pwcp(ji,jk,jwalk) + zreasat 
    491240         ENDDO 
    492241      ENDDO 
    … …  
    496245     ! Patankar scheme 
    497246      ! ------------------------------------------------------------------------------ 
    498  
    499247      call sed_dsr_redoxb 
    500248 
    … …  
    564312      !! Arguments 
    565313      ! --- local variables 
    566       INTEGER   ::  ji, jk, jn   ! dummy looop indices 
     314      INTEGER   ::  ji, jk, jn, jw   ! dummy looop indices 
    567315 
    568316      REAL, DIMENSION(6)  :: zsedtrn, zsedtra 
    … …  
    582330            zsedtrn(5)  = solcp(ji,jk,jsfeo) * zvolc(ji,jk,jsfeo) 
    583331            zsedtrn(6)  = solcp(ji,jk,jsfes) * zvolc(ji,jk,jsfes) 
    584             zsedtra(:)  = zsedtrn(:)  
     332            zsedtra(:)  = MAX(0., zsedtrn(:) - rtrn ) 
    585333 
    586334            ! First pass of the scheme. At the end, it is 1st order  
    … …  
    592340            zdelta = pwcp(ji,jk,jwpo4) - redfep * zsedtra(4) 
    593341            IF ( zalpha == 0. ) THEN 
    594                zsedtra(4) = zsedtra(4) / ( 1.0 + zsedtra(4) * reac_fe2 * dtsed2 / 2.0 ) 
    595             ELSE 
    596                zsedtra(4) = ( zsedtra(4) * zalpha ) / ( 0.25 * zsedtra(4) * ( exp( reac_fe2 * zalpha * dtsed2 / 2. ) - 1.0 )  & 
     342               zsedtra(4) = zsedtra(4) / ( 1.0 + 0.25 * zsedtra(4) * reac_fe2 * dtsed2 / 2.0 ) 
     343            ELSE 
     344               zsedtra(4) = ( zsedtra(4) * zalpha ) / ( 0.25 * zsedtra(4)    & 
     345               &            * ( exp( reac_fe2 * zalpha * dtsed2 / 2. ) - 1.0 )  & 
    597346               &            + zalpha * exp( reac_fe2 * zalpha * dtsed2 / 2. ) ) 
    598347            ENDIF 
    599             zsedtra(1) = zalpha + 0.25 * zsedtra(4)  
     348            zsedtra(1) = zalpha + 0.25 * zsedtra(4) 
    600349            zsedtra(5) = zbeta  - zsedtra(4) 
    601350            pwcp(ji,jk,jwalk) = zgamma + 2.0 * zsedtra(4) 
    … …  
    606355            zgamma = pwcp(ji,jk,jwalk) - 2.0 * zsedtra(2) 
    607356            IF ( zalpha == 0. ) THEN 
    608                zsedtra(2) = zsedtra(2) / ( 1.0 + zsedtra(2) * reac_h2s * dtsed2 / 2.0 ) 
     357               zsedtra(2) = zsedtra(2) / ( 1.0 + 2.0 * zsedtra(2) * reac_h2s * dtsed2 / 2.0 ) 
    609358            ELSE 
    610359               zsedtra(2) = ( zsedtra(2) * zalpha ) / ( 2.0 * zsedtra(2) * ( exp( reac_h2s * zalpha * dtsed2 / 2. ) - 1.0 )  & 
    … …  
    615364            pwcp(ji,jk,jwso4) = zbeta - zsedtra(2) 
    616365            ! NH4 + O2 
    617             zalpha = zsedtra(1) - 2.0 * zsedtra(3) 
    618             zgamma = pwcp(ji,jk,jwalk) - 2.0 * zsedtra(3) 
    619             IF ( zalpha == 0. ) THEN 
    620                zsedtra(3) = zsedtra(3) / ( 1.0 + zsedtra(3) * reac_nh4 * zadsnh4 * dtsed2 / 2.0 ) 
    621             ELSE 
    622                zsedtra(3) = ( zsedtra(3) * zalpha ) / ( 2.0 * zsedtra(3) * ( exp( reac_nh4 * zadsnh4 * zalpha * dtsed2 / 2. ) - 1.0 )  & 
    623                &            + zalpha * exp( reac_nh4 * zadsnh4 * zalpha * dtsed2 /2. ) ) 
    624             ENDIF 
    625             zsedtra(1) = zalpha + 2.0 * zsedtra(3) 
    626             pwcp(ji,jk,jwalk) = zgamma + 2.0 * zsedtra(3) 
     366            zalpha = zsedtra(1) - 2.0 * zsedtra(3) / zadsnh4  
     367            zgamma = pwcp(ji,jk,jwalk) - 2.0 * zsedtra(3) /zadsnh4  
     368            IF ( zalpha == 0. ) THEN 
     369               zsedtra(3) = zsedtra(3) / ( 1.0 + 2.0 * zsedtra(3) * reac_nh4 / zadsnh4 * dtsed2 / 2.0 ) 
     370            ELSE 
     371               zsedtra(3) = ( zsedtra(3) * zalpha * zadsnh4 ) / ( 2.0 * zsedtra(3) * ( exp( reac_nh4 * zadsnh4 * zalpha * dtsed2 / 2. ) - 1.0 )  & 
     372               &            + zalpha * zadsnh4 * exp( reac_nh4 * zadsnh4 * zalpha * dtsed2 /2. ) ) 
     373            ENDIF 
     374            zsedtra(1) = zalpha + 2.0 * zsedtra(3) / zadsnh4  
     375            pwcp(ji,jk,jwalk) = zgamma + 2.0 * zsedtra(3) / zadsnh4  
    627376            ! FeS - O2 
    628377            zalpha = zsedtra(1) - 2.0 * zsedtra(6) 
    … …  
    630379            zgamma = pwcp(ji,jk,jwso4) + zsedtra(6) 
    631380            IF ( zalpha == 0. ) THEN 
    632                zsedtra(6) = zsedtra(6) / ( 1.0 + zsedtra(6) * reac_feso * dtsed2 / 2.0 ) 
     381               zsedtra(6) = zsedtra(6) / ( 1.0 + 2.0 * zsedtra(6) * reac_feso * dtsed2 / 2.0 ) 
    633382            ELSE 
    634383               zsedtra(6) = ( zsedtra(6) * zalpha ) / ( 2.0 * zsedtra(6) * ( exp( reac_feso * zalpha * dtsed2 / 2. ) - 1.0 )  & 
    … …  
    658407            zepsi  = pwcp(ji,jk,jwpo4) + redfep * zsedtra(5) 
    659408            IF ( zalpha == 0. ) THEN 
    660                zsedtra(2) = zsedtra(2) / ( 1.0 + zsedtra(2) * reac_feh2s * dtsed2 ) 
     409               zsedtra(2) = zsedtra(2) / ( 1.0 + 2.0 * zsedtra(2) * reac_feh2s * dtsed2 ) 
    661410            ELSE 
    662411               zsedtra(2) = ( zsedtra(2) * zalpha ) / ( 2.0 * zsedtra(2) * ( exp( reac_feh2s * zalpha * dtsed2 ) - 1.0 )  & 
    … …  
    668417            pwcp(ji,jk,jwalk) = zgamma + 2.0 * zsedtra(4) 
    669418            pwcp(ji,jk,jwpo4) = zepsi - redfep * zsedtra(5) 
    670             ! Fe - H2S 
     419           ! Fe - H2S 
    671420            zalpha = zsedtra(2) - zsedtra(4) 
    672421            zbeta  = zsedtra(4) + zsedtra(6) 
    … …  
    686435            zgamma = pwcp(ji,jk,jwso4) + zsedtra(6) 
    687436            IF (zalpha == 0.) THEN 
    688                zsedtra(6) = zsedtra(6) / ( 1.0 + zsedtra(6) * reac_feso * dtsed2 / 2. ) 
     437               zsedtra(6) = zsedtra(6) / ( 1.0 + 2.0 * zsedtra(6) * reac_feso * dtsed2 / 2. ) 
    689438            ELSE 
    690439               zsedtra(6) = ( zsedtra(6) * zalpha ) / ( 2.0 * zsedtra(6) * ( exp( reac_feso * zalpha * dtsed2 / 2. ) - 1.0 )  & 
    … …  
    694443            zsedtra(4) = zbeta  - zsedtra(6) 
    695444            pwcp(ji,jk,jwso4) = zgamma - zsedtra(6) 
     445 
    696446            ! NH4 + O2 
    697             zalpha = zsedtra(1) - 2.0 * zsedtra(3) 
    698             zgamma = pwcp(ji,jk,jwalk) - 2.0 * zsedtra(3) 
    699             IF (zalpha == 0.) THEN 
    700                zsedtra(3) = zsedtra(3) / ( 1.0 + zsedtra(3) * reac_nh4 * zadsnh4 * dtsed2 / 2.0)  
    701             ELSE 
    702                zsedtra(3) = ( zsedtra(3) * zalpha ) / ( 2.0 * zsedtra(3) * ( exp( reac_nh4 * zadsnh4 * zalpha * dtsed2 / 2. ) - 1.0 )  & 
    703                &            + zalpha * exp( reac_nh4 * zadsnh4 * zalpha * dtsed2 /2. ) ) 
    704             ENDIF 
    705             zsedtra(1) = zalpha + 2.0 * zsedtra(3) 
    706             pwcp(ji,jk,jwalk) = zgamma + 2.0 * zsedtra(3) 
     447            zalpha = zsedtra(1) - 2.0 * zsedtra(3) / zadsnh4  
     448            zgamma = pwcp(ji,jk,jwalk) - 2.0 * zsedtra(3) / zadsnh4  
     449            IF ( zalpha == 0. ) THEN 
     450               zsedtra(3) = zsedtra(3) / ( 1.0 + 2.0 * zsedtra(3) * reac_nh4 / zadsnh4 * dtsed2 / 2.0 ) 
     451            ELSE 
     452               zsedtra(3) = ( zsedtra(3) * zalpha * zadsnh4 ) / ( 2.0 * zsedtra(3) * ( exp( reac_nh4 * zadsnh4 * zalpha * dtsed2 / 2. ) - 1.0 )  & 
     453               &            + zalpha * zadsnh4 * exp( reac_nh4 * zadsnh4 * zalpha * dtsed2 /2. ) ) 
     454            ENDIF 
     455            zsedtra(1) = zalpha + 2.0 * zsedtra(3) / zadsnh4  
     456            pwcp(ji,jk,jwalk) = zgamma + 2.0 * zsedtra(3) / zadsnh4  
    707457            ! H2S + O2 
    708458            zalpha = zsedtra(1) - 2.0 * zsedtra(2) 
    … …  
    710460            zgamma = pwcp(ji,jk,jwalk) - 2.0 * zsedtra(2) 
    711461            IF ( zalpha == 0. ) THEN 
    712                zsedtra(2) = zsedtra(2) / ( 1.0 + zsedtra(2) * reac_h2s * dtsed2 / 2.0 ) 
     462               zsedtra(2) = zsedtra(2) / ( 1.0 + 2.0 * zsedtra(2) * reac_h2s * dtsed2 / 2.0 ) 
    713463            ELSE 
    714464               zsedtra(2) = ( zsedtra(2) * zalpha ) / ( 2.0 * zsedtra(2) * ( exp( reac_h2s * zalpha * dtsed2 / 2. ) - 1.0 )  & 
    … …  
    718468            pwcp(ji,jk,jwso4) = zbeta - zsedtra(2) 
    719469            pwcp(ji,jk,jwalk) = zgamma + 2.0 * zsedtra(2) 
     470 
    720471            ! Fe + O2 
    721472            zalpha = zsedtra(1) - 0.25 * zsedtra(4) 
    … …  
    724475            zdelta = pwcp(ji,jk,jwpo4) - redfep * zsedtra(4) 
    725476            IF ( zalpha == 0. ) THEN 
    726                zsedtra(4) = zsedtra(4) / ( 1.0 + zsedtra(4) * reac_fe2 * dtsed2 / 2.0 ) 
    727             ELSE 
    728                zsedtra(4) = ( zsedtra(4) * zalpha ) / ( 0.25 * zsedtra(4) * ( exp( reac_fe2 * zalpha * dtsed2 / 2. ) - 1.0 )  & 
     477               zsedtra(4) = zsedtra(4) / ( 1.0 + 0.25 * zsedtra(4) * reac_fe2 * dtsed2 / 2.0 ) 
     478            ELSE 
     479               zsedtra(4) = ( zsedtra(4) * zalpha ) / ( 0.25 * zsedtra(4)     & 
     480               &            * ( exp( reac_fe2  * zalpha * dtsed2 / 2. ) - 1.0 )  & 
    729481               &            + zalpha * exp( reac_fe2 * zalpha * dtsed2 / 2. ) ) 
    730482            ENDIF 
    … …  
    733485            pwcp(ji,jk,jwpo4) = zdelta + redfep * zsedtra(4) 
    734486            pwcp(ji,jk,jwalk) = zgamma + 2.0 * zsedtra(4) 
     487 
     488            ! Update the concentrations after the secondary reactions  
     489            zsedtra(:) = zsedtra(:) + MIN( zsedtrn(:), rtrn ) 
    735490            pwcp(ji,jk,jwoxy)  = zsedtra(1) 
    736491            pwcp(ji,jk,jwh2s)  = zsedtra(2) 

  • NEMO/branches/2019/dev_r11708_aumont_PISCES_QUOTA/src/TOP/PISCES/SED/seddta.F90

    r14306 r14323  
    4949 
    5050      REAL(wp), DIMENSION(jpoce) :: zdtap, zdtag 
    51       REAL(wp), DIMENSION(jpi,jpj) :: zwsbio4, zwsbio3 
     51      REAL(wp), DIMENSION(jpi,jpj) :: zwsbio4, zwsbio3, zddust 
    5252      REAL(wp) :: zf0, zf1, zf2, zkapp, zratio, zdep 
    5353 
    … …  
    7676!         conv2   = 1.0e+3 / ( 1.0e+4 * rsecday * 30. ) 
    7777         conv2 = 1.0e+3 /  1.0e+4  
    78          rdtsed(2:jpksed) = dtsed / ( denssol * por1(2:jpksed) ) 
    7978      ENDIF 
    8079 
    … …  
    8281      zdtap(:)    = 0.  
    8382      zdtag(:)    = 0.   
     83      zddust(:,:) = 0.0 
    8484 
    8585      ! reading variables 
    … …  
    157157!        zf1    = ( 0.02 * (reac_poc - reac_poc * 0.) + zkapp - reac_poc ) / ( reac_poc / 100. - reac_poc ) 
    158158!        zf1    = MIN(0.98, MAX(0., zf1 ) ) 
    159          zf1    = 0.7 
     159         zf1    = 0.48 
    160160         zf2    = 0.02 
    161161         zf0    = 1.0 - zf2 - zf1 
    … …  
    180180      rainrm_dta(1:jpoce,jsfeo)  = rainrm_dta(1:jpoce,jsclay) * mol_wgt(jsclay) / mol_wgt(jsfeo) * 0.035 * 0.5 * 0.333 
    181181      rainrm_dta(1:jpoce,jsclay) = rainrm_dta(1:jpoce,jsclay) * (1.0 - 0.035 * 0.5 * 0.333 ) 
     182      CALL unpack_arr ( jpoce, zddust(1:jpi,1:jpj), iarroce(1:jpoce), wacc(1:jpoce) ) 
     183      zddust(:,:) = dust(:,:) + zddust(:,:) / ( rsecday * 365.0 ) * por1(2) * denssol / conv2 
    182184 
    183185!    rainrm_dta(1:jpoce,jsclay) = 1.0E-4 * conv2 / mol_wgt(jsclay) 
    … …  
    213215 
    214216      ! computation of dzdep = total thickness of solid material rained [cm] in each cell 
    215       dzdep(1:jpoce) = raintg(1:jpoce) * rdtsed(2)  
     217      dzdep(1:jpoce) = raintg(1:jpoce) * dtsed / ( denssol * por1(2) ) 
    216218 
    217219      IF( lk_iomput ) THEN 
    218           IF( iom_use("sflxclay" ) ) CALL iom_put( "sflxclay", dust(:,:) * 1E3 / 1.E4 ) 
     220          IF( iom_use("sflxclay" ) ) CALL iom_put( "sflxclay", zddust(:,:) * 1E3 / 1.E4 ) 
    219221          IF( iom_use("sflxcal" ) )  CALL iom_put( "sflxcal", trc_data(:,:,14) / 1.E4 ) 
    220222          IF( iom_use("sflxbsi" ) )  CALL iom_put( "sflxbsi", trc_data(:,:,11) / 1.E4 ) 

  • NEMO/branches/2019/dev_r11708_aumont_PISCES_QUOTA/src/TOP/PISCES/SED/sedini.F90

    r14306 r14323  
    161161      ALLOCATE( diff(jpoce,jpksed,jpwat ) )  ;  ALLOCATE( irrig(jpoce, jpksed) ) 
    162162      ALLOCATE( wacc(jpoce) )                ;  ALLOCATE( fecratio(jpoce) ) 
    163       ALLOCATE( press(jpoce) )               ;  ALLOCATE( densSW(jpoce) )  
    164163      ALLOCATE( hipor(jpoce,jpksed) )        ;  ALLOCATE( co3por(jpoce,jpksed) ) 
    165164      ALLOCATE( dz3d(jpoce,jpksed) )         ;  ALLOCATE( volw3d(jpoce,jpksed) )       ;  ALLOCATE( vols3d(jpoce,jpksed) ) 
    166       ALLOCATE( sedligand(jpoce, jpksed) ) 
     165      ALLOCATE( sedligand(jpoce, jpksed) )   ;  ALLOCATE( saturco3(jpoce,jpksed) )  ;  ALLOCATE( densSW(jpoce) ) 
    167166 
    168167      ! Initialization of global variables 
    169       pwcp  (:,:,:) = 0.   ;  pwcp0 (:,:,:) = 0.  ; pwcp_dta  (:,:) = 0.   
    170       solcp (:,:,:) = 0.   ;  solcp0(:,:,:) = 0.  ; rainrm_dta(:,:) = 0. 
    171       rainrm(:,:  ) = 0.   ;  rainrg(:,:  ) = 0.  ; raintg    (:  ) = 0.  
    172       dzdep (:    ) = 0.   ;  iarroce(:   ) = 0   ; dzkbot    (:  ) = 0. 
    173       temp  (:    ) = 0.   ;  salt   (:   ) = 0.  ; zkbot     (:  ) = 0. 
    174       press (:    ) = 0.   ;  densSW (:   ) = 0.  ; db        (:,:) = 0.  
    175       hipor (:,:  ) = 0.   ;  co3por (:,: ) = 0.  ; irrig     (:,:) = 0.  
    176       dz3d  (:,:  ) = 0.   ;  volw3d (:,: ) = 0.  ; vols3d    (:,:) = 0.  
     168      pwcp  (:,:,:) = 0.   ;  pwcp0  (:,:,:) = 0.  ; pwcp_dta  (:,:) = 0.   
     169      solcp (:,:,:) = 0.   ;  solcp0 (:,:,:) = 0.  ; rainrm_dta(:,:) = 0. 
     170      rainrm(:,:  ) = 0.   ;  rainrg (:,:  ) = 0.  ; raintg    (:  ) = 0.  
     171      dzdep (:    ) = 0.   ;  iarroce(:   )  = 0   ; dzkbot    (:  ) = 0. 
     172      temp  (:    ) = 0.   ;  salt   (:   )  = 0.  ; zkbot     (:  ) = 0. 
     173      db    (:,:  ) = 0.   ;  densSW (:   ) = 0.  
     174      hipor (:,:  ) = 0.   ;  co3por (:,: )  = 0.  ; irrig     (:,:) = 0.  
     175      dz3d  (:,:  ) = 0.   ;  volw3d (:,: )  = 0.  ; vols3d    (:,:) = 0.  
    177176      fecratio(:)   = 1E-5  
    178177      sedligand(:,:) = 0.6E-9 
    179178 
    180179      ! Chemical variables       
    181       ALLOCATE( akbs  (jpoce) )  ;  ALLOCATE( ak1s   (jpoce) )  ;  ALLOCATE( ak2s  (jpoce) ) ;  ALLOCATE( akws  (jpoce) )      
    182       ALLOCATE( ak1ps (jpoce) )  ;  ALLOCATE( ak2ps  (jpoce) )  ;  ALLOCATE( ak3ps (jpoce) ) ;  ALLOCATE( aksis (jpoce) )     
    183       ALLOCATE( aksps (jpoce) )  ;  ALLOCATE( ak12s  (jpoce) )  ;  ALLOCATE( ak12ps(jpoce) ) ;  ALLOCATE( ak123ps(jpoce) )     
    184       ALLOCATE( borats(jpoce) )  ;  ALLOCATE( calcon2(jpoce) )  ;  ALLOCATE( sieqs (jpoce) )  
     180      ALLOCATE( akbs  (jpoce) )  ;  ALLOCATE( ak1s   (jpoce) )  ;  ALLOCATE( ak2s  (jpoce) ) ;  ALLOCATE( akws  (jpoce) ) 
     181      ALLOCATE( ak1ps (jpoce) )  ;  ALLOCATE( ak2ps  (jpoce) )  ;  ALLOCATE( ak3ps (jpoce) ) ;  ALLOCATE( aksis (jpoce) ) 
     182      ALLOCATE( aksps (jpoce) )  ;  ALLOCATE( ak12s  (jpoce) )  ;  ALLOCATE( ak12ps(jpoce) ) ;  ALLOCATE( ak123ps(jpoce) ) 
     183      ALLOCATE( borats(jpoce) )  ;  ALLOCATE( calcon2(jpoce) )  ;  ALLOCATE( sieqs (jpoce) ) 
    185184      ALLOCATE( aks3s(jpoce) )   ;  ALLOCATE( akf3s(jpoce) )    ;  ALLOCATE( sulfats(jpoce) ) 
    186185      ALLOCATE( fluorids(jpoce) ) 
    … …  
    193192 
    194193      ! Mass balance calculation   
    195       ALLOCATE( fromsed(jpoce, jpsol) ) ; ALLOCATE( tosed(jpoce, jpsol) ) ;  ALLOCATE( rloss(jpoce, jpsol) ) 
    196       ALLOCATE( tokbot (jpoce, jpwat) )  
    197  
    198       fromsed(:,:) = 0.    ;   tosed(:,:) = 0. ;  rloss(:,:) = 0.  ;   tokbot(:,:) = 0.  
     194      ALLOCATE( fromsed(jpoce, jpsol) ) ; ALLOCATE( tosed(jpoce, jpsol) ) 
     195 
     196      fromsed(:,:) = 0.    ;   tosed(:,:) = 0. 
    199197 
    200198      ! Initialization of sediment geometry 

  • NEMO/branches/2019/dev_r11708_aumont_PISCES_QUOTA/src/TOP/PISCES/SED/sedinorg.F90

    r12360 r14323  
    2323CONTAINS 
    2424    
    25    SUBROUTINE sed_inorg( kt )  
     25   SUBROUTINE sed_inorg( kt, knt )  
    2626      !!---------------------------------------------------------------------- 
    2727      !!                   ***  ROUTINE sed_inorg  *** 
    … …  
    4646      !!---------------------------------------------------------------------- 
    4747      !! Arguments 
    48       INTEGER, INTENT(in) ::   kt       ! number of iteration 
     48      INTEGER, INTENT(in) ::   kt, knt       ! number of iteration 
    4949      ! --- local variables 
    5050      INTEGER :: ji, jk, js, jw         ! dummy looop indices 
    … …  
    5454      REAL(wp), DIMENSION(jpoce,jpksed,jpsol) :: zvolc    ! temp. variables 
    5555      REAL(wp), DIMENSION(jpoce) :: zsieq 
    56       REAL(wp)  ::  zsolid1, zvolw, zreasat 
     56      REAL(wp)  ::  zsolid1, zvolw, zreasat, zbeta, zgamma 
    5757      REAL(wp)  ::  zsatur, zsatur2, znusil, zsolcpcl, zsolcpsi 
    5858      !! 
    … …  
    6161      IF( ln_timing )  CALL timing_start('sed_inorg') 
    6262! 
    63       IF( kt == nitsed000 ) THEN 
     63      IF( kt == nitsed000 .AND. knt == 1 ) THEN 
    6464         IF (lwp) THEN 
    6565            WRITE(numsed,*) ' sed_inorg : Dissolution reaction ' 
    … …  
    7373       
    7474      zrearat1(:,:) = 0.    ;   zundsat(:,:)  = 0.  
    75       zrearat2(:,:) = 0.    ;   zrearat2(:,:) = 0. 
     75      zrearat2(:,:) = 0. 
    7676      zco3eq(:)     = rtrn 
    7777      zvolc(:,:,:)  = 0. 
    … …  
    8686         zsolcpsi = 0.0 
    8787         DO jk = 1, jpksed 
    88             zsolcpsi = zsolcpsi + solcp(ji,jk,jsopal) * dz(jk) 
    89             zsolcpcl = zsolcpcl + solcp(ji,jk,jsclay) * dz(jk) 
     88            zsolcpsi = zsolcpsi + solcp(ji,jk,jsopal) * vols3d(ji,jk) 
     89            zsolcpcl = zsolcpcl + solcp(ji,jk,jsclay) * vols3d(ji,jk) 
    9090         END DO 
    9191         zsolcpsi = MAX( zsolcpsi, rtrn ) 
    … …  
    135135      DO jk = 2, jpksed 
    136136         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) 
     137            IF ( zundsat(ji,jk) > 0. ) THEN 
     138               zsolid1 = zvolc(ji,jk,jsopal) * solcp(ji,jk,jsopal) 
     139               zsatur = zundsat(ji,jk) / zsieq(ji) 
     140               zsatur2 = (1.0 + temp(ji) / 400.0 )**37 
     141               znusil = ( 0.225 * ( 1.0 + temp(ji) / 15.) + 0.775 * zsatur2 * zsatur**8.25 ) / zsieq(ji) 
     142               zbeta = 1.0 - reac_sil * znusil * dtsed2 * zundsat(ji,jk) + reac_sil * znusil * dtsed2 * zsolid1 
     143               zgamma = - reac_sil * znusil * dtsed2 * zundsat(ji,jk) 
     144               zundsat(ji,jk) = ( - zbeta + SQRT( zbeta**2 - 4.0 * zgamma ) ) / ( 2.0 * reac_sil * znusil * dtsed2 ) 
     145               solcp(ji,jk,jsopal ) = zsolid1 / ( 1. + reac_sil * znusil * dtsed2 * zundsat(ji,jk) ) / zvolc(ji,jk,jsopal) 
     146               pwcp(ji,jk,jwsil) = zsieq(ji) - zundsat(ji,jk)  
     147            ENDIF 
    161148         ENDDO 
    162149      ENDDO 
    … …  
    175162            zco3eq(ji)     = aksps(ji) * densSW(ji) * densSW(ji) / ( calcon2(ji) + rtrn ) 
    176163            zco3eq(ji)     = MAX( rtrn, zco3eq(ji) )  
    177             zundsat(ji,jk) = MAX(0., zco3eq(ji) - co3por(ji,jk) ) 
     164            zundsat(ji,jk) = ( zco3eq(ji) - co3por(ji,jk) ) / zco3eq(ji) 
    178165         ENDDO 
    179166      ENDDO 
    … …  
    181168      DO jk = 2, jpksed 
    182169         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) ) 
     170            IF ( zundsat(ji,jk) > 0. ) THEN 
     171               zsolid1 = zvolc(ji,jk,jscal) * solcp(ji,jk,jscal) 
     172               zbeta = 1.0 - reac_cal * dtsed2 * zundsat(ji,jk) + reac_cal * dtsed2 * zsolid1 
     173               zgamma = -reac_cal * dtsed2 * zundsat(ji,jk) 
     174               zundsat(ji,jk) = ( -zbeta + SQRT( zbeta**2 - 4.0 * zgamma ) ) / ( 2.0 * reac_sil * znusil * dtsed2 ) 
     175               saturco3(ji,jk) = zundsat(ji,jk) 
     176               zreasat = reac_cal * dtsed2 * zundsat(ji,jk) * zsolid1 / ( 1. + reac_cal * dtsed2 * zundsat(ji,jk) ) 
     177               solcp(ji,jk,jscal) = solcp(ji,jk,jscal) - zreasat / zvolc(ji,jk,jscal) 
     178               ! For DIC 
     179               pwcp(ji,jk,jwdic)  = pwcp(ji,jk,jwdic) + zreasat 
     180               ! For alkalinity 
     181               pwcp(ji,jk,jwalk)  = pwcp(ji,jk,jwalk) + 2.0 * zreasat 
     182            ENDIF 
    186183         END DO 
    187184      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 
    296185 
    297186      IF( ln_timing )  CALL timing_stop('sed_inorg') 

  • NEMO/branches/2019/dev_r11708_aumont_PISCES_QUOTA/src/TOP/PISCES/SED/sedstp.F90

    r10222 r14323  
    88   USE sedchem  ! chemical constant 
    99   USE sedco3   ! carbonate in sediment pore water 
    10    USE sedorg   ! Organic reactions and diffusion 
    11    USE sedinorg ! Inorganic dissolution 
     10   USE sedsol   ! Organic reactions and diffusion 
    1211   USE sedbtb   ! bioturbation 
    1312   USE sedadv   ! vertical advection 
    14    USE sedmbc   ! mass balance calculation 
    1513   USE sedsfc   ! sediment surface data 
    1614   USE sedrst   ! restart 
    … …  
    4543      !!---------------------------------------------------------------------- 
    4644      INTEGER, INTENT(in) ::   kt       ! number of iteration 
    47       INTEGER :: ji,jk,js,jn,jw 
    4845      !!---------------------------------------------------------------------- 
    49       IF( ln_timing )      CALL timing_start('sed_stp') 
     46      IF( ln_timing )           CALL timing_start('sed_stp') 
    5047        ! 
    5148                                CALL sed_rst_opn  ( kt )       ! Open tracer restart file  
    … …  
    6562         ENDIF 
    6663 
    67          CALL sed_btb( kt )         ! 1st pass of bioturbation at t+1/2 
    68          CALL sed_org( kt )         ! Organic related reactions and diffusion 
    69          CALL sed_inorg( kt )       ! Dissolution reaction 
    70          CALL sed_btb( kt )         ! 2nd pass of bioturbation at t+1 
    71          tokbot(:,:) = 0.0 
    72          DO jw = 1, jpwat 
    73             DO ji = 1, jpoce 
    74                tokbot(ji,jw) = pwcp(ji,1,jw) * 1.e-3 * dzkbot(ji) 
    75             END DO 
    76          ENDDO 
     64         CALL sed_btb( kt )        ! 1st pass of bioturbation at t+1/2 
     65         CALL sed_sol( kt )        ! Solute diffusion and reactions  
     66         CALL sed_btb( kt )        ! 2nd pass of bioturbation at t+1 
    7767         CALL sed_adv( kt )         ! advection 
    7868         CALL sed_co3( kt )         ! pH actualization for saving 
    79          ! This routine is commented out since it does not work at all 
    80          CALL sed_mbc( kt )         ! cumulation for mass balance calculation 
    81  
    8269         IF (ln_sed_2way) CALL sed_sfc( kt )         ! Give back new bottom wat chem to tracer model 
    8370      ENDIF 

  • NEMO/branches/2019/dev_r11708_aumont_PISCES_QUOTA/src/TOP/PISCES/SED/sedwri.F90

    r14306 r14323  
    9090      CALL unpack_arr( jpoce, flxsedi3d(1:jpi,1:jpj,1:jpksed,3)  , iarroce(1:jpoce), & 
    9191         &                   sedligand(1:jpoce,1:jpksed)  ) 
     92 
     93      CALL unpack_arr( jpoce, flxsedi3d(1:jpi,1:jpj,1:jpksed,4)  , iarroce(1:jpoce), & 
     94         &                   saturco3(1:jpoce,1:jpksed)  ) 
     95 
    9296       
    9397!      flxsedi3d = 0. 
    … …  
    97101         DO ji = 1, jpoce 
    98102            zflx(ji,jw) = ( pwcp(ji,1,jw) - pwcp_dta(ji,jw) ) & 
    99                &         * 1.e3 / ( 1.e2 * dzkbot(ji) ) / 1.E4 / r2dttrc 
     103               &         * 1.e3 * ( 1.e-2 * dzkbot(ji) ) / 1.E4 / r2dttrc 
    100104         ENDDO 
    101105      ENDDO 
Note: See TracChangeset for help on using the changeset viewer.