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

Changeset 15450


Ignore:
Timestamp:
2021-10-27T16:32:08+02:00 (7 months ago)
Author:
cetlod
Message:

Some updates to make the PISCES/SED module usable. Totally orthogonal with no effect on other parts of the code

Location:
NEMO/trunk/src/TOP/PISCES
Files:
5 added
1 deleted
20 edited

Legend:

Unmodified
Added
Removed
  • NEMO/trunk/src/TOP/PISCES/SED/par_sed.F90

    r15033 r15450  
    1818      jp_sal   =>   jp_sal     !: indice of salintity 
    1919 
    20    INTEGER, PUBLIC, PARAMETER :: jpdta = 17 
     20   INTEGER, PARAMETER :: jpdta = 18 
    2121 
    2222   ! Vertical sediment geometry 
    23    INTEGER, PUBLIC  :: jpksed  = 11  
    24    INTEGER, PUBLIC  :: jpksedm1  
     23   INTEGER, PUBLIC   ::      & 
     24      jpksed   = 11 ,        & 
     25      jpksedm1 = 10 
    2526 
    2627   ! sediment tracer species 
    27    INTEGER, PUBLIC, PARAMETER ::    & 
     28   INTEGER, PARAMETER ::    & 
    2829      jpsol =  8,           &  !: number of solid component 
    29       jpwat = 10,           &   !: number of pore water component 
     30      jpwat = 11,           &   !: number of pore water component 
     31      jpads = 2 ,           &   !! number adsorbed species 
    3032      jpwatp1 = jpwat +1,   & 
    3133      jpsol1 = jpsol - 1 
     
    3335    
    3436   ! pore water components        
    35    INTEGER, PUBLIC, PARAMETER :: & 
    36       jwsil  = 1,        & !: silic acid 
    37       jwoxy  = 2,        & !: oxygen 
    38       jwdic  = 3,        & !: dissolved inorganic carbon 
    39       jwno3  = 4,        & !: nitrate 
    40       jwpo4  = 5,        & !: phosphate 
    41       jwalk  = 6,        & !: alkalinity 
    42       jwnh4  = 7,        & !: Ammonium 
    43       jwh2s  = 8,        & !: Sulfate 
    44       jwso4  = 9,        & !: H2S 
    45       jwfe2  = 10          !: Fe2+ 
     37   INTEGER, PARAMETER :: & 
     38      jwoxy  = 1,        & !: oxygen 
     39      jwno3  = 2,        & !: nitrate 
     40      jwpo4  = 3,        & !: phosphate 
     41      jwnh4  = 4,        & !: Ammonium 
     42      jwh2s  = 5,        & !: Sulfate 
     43      jwso4  = 6,        & !: H2S 
     44      jwfe2  = 7,        & !: Fe2+ 
     45      jwalk  = 8,        & !: Alkalinity 
     46      jwlgw  = 9,        & !: Alkalinity 
     47      jwdic  = 10,       & !: DIC 
     48      jwsil  = 11          !: Silicate 
    4649 
    4750   ! solid components        
    48    INTEGER, PUBLIC, PARAMETER ::  & 
    49       jsopal  = 1,        & !: opal sediment 
    50       jsclay  = 2,        & !: clay 
     51   INTEGER, PARAMETER ::  & 
     52      jsfeo   = 1,        & !: iron hydroxides 
     53      jsfes   = 2,        & !: FeS 
    5154      jspoc   = 3,        & !: organic carbon 
    52       jscal   = 4,        & !: calcite 
    53       jspos   = 5,        & !: semi-ref POC 
    54       jspor   = 6,        & !: refractory POC 
    55       jsfeo   = 7,        & !: iron hydroxides 
    56       jsfes   = 8           !: FeS 
     55      jspos   = 4,        & !: semi-ref POC 
     56      jspor   = 5,        & !: refractory POC 
     57      jscal   = 6,        & !: Calcite 
     58      jsopal  = 7,        & !: Opal 
     59      jsclay  = 8           !: clay 
    5760 
    58    INTEGER, PUBLIC, PARAMETER ::  & 
     61   INTEGER, PARAMETER ::  & 
    5962      jptrased   = jpsol + jpwat , & 
    60       jpdia3dsed = 2             , & 
    61       jpdia2dsed = 12 
     63      jpvode     = jptrased - 11  , & 
     64      jpdia3dsed = 3             , & 
     65      jpdia2dsed = 22  
    6266 
    6367END MODULE par_sed 
  • NEMO/trunk/src/TOP/PISCES/SED/sed.F90

    r14086 r15450  
    77   !!        !  06-12  (C. Ethe)  Orignal 
    88   !!---------------------------------------------------------------------- 
     9   USE par_sed 
    910   USE oce_sed 
    1011   USE in_out_manager 
     
    2728   REAL(wp), PUBLIC               ::  reac_fe2            !: reactivity of Fe2+ in  [l.mol-1.s-1] 
    2829   REAL(wp), PUBLIC               ::  reac_feh2s          !: reactivity of Fe2+ in  [l.mol-1.s-1] 
    29    REAL(wp), PUBLIC               ::  reac_fes            !: reactivity of Fe with H2S in  [l.mol-1.s-1] 
    3030   REAL(wp), PUBLIC               ::  reac_feso           !: reactivity of FeS with O2 in  [l.mol-1.s-1] 
     31   REAL(wp), PUBLIC               ::  reac_fesp           !: precipitation of FeS  [mol.l-1.s-1] 
     32   REAL(wp), PUBLIC               ::  reac_fesd           !: Dissolution  of FeS  [s-1] 
    3133   REAL(wp), PUBLIC               ::  reac_cal            !: reactivity of cal in  [l.mol-1.s-1] 
    3234   REAL(wp), PUBLIC               ::  adsnh4              !: adsorption coefficient of NH4 
     35   REAL(wp), PUBLIC               ::  adsfe2              !: adsorption coefficient of Fe2 
    3336   REAL(wp), PUBLIC               ::  ratligc             !: C/L ratio in POC 
    3437   REAL(wp), PUBLIC               ::  so2ut  
     
    3740   REAL(wp), PUBLIC               ::  srDnit  
    3841   REAL(wp), PUBLIC               ::  dtsed               !: sedimentation time step 
    39    REAL(wp), PUBLIC               ::  dtsed2              !: sedimentation time step 
    4042   INTEGER , PUBLIC               ::  nitsed000 
    4143   INTEGER , PUBLIC               ::  nitsedend 
    42    INTEGER, PUBLIC                ::  nrseddt 
    43    REAL    , PUBLIC               ::  sedmask 
    44    REAL(wp), PUBLIC               ::  denssol                !: density of solid material 
    4544   LOGICAL , PUBLIC               ::  lrst_sed       !: logical to control the trc restart write 
    4645   LOGICAL , PUBLIC               ::  ln_rst_sed  = .TRUE.     !: initialisation from a restart file or not 
     
    5554   CHARACTER(len = 80) , PUBLIC   ::  cn_sedrst_out  !: suffix of pass. tracer restart name (output) 
    5655   CHARACTER(len = 256), PUBLIC   ::  cn_sedrst_outdir  !: restart output directory 
     56   INTEGER, PUBLIC                ::  nrosorder  !: order of the rosenbrock method 
     57   REAL(wp), PUBLIC               ::  rosatol   !: Tolerance for absolute error 
     58   REAL(wp), PUBLIC               ::  rosrtol   !: Tolerance for relative error 
    5759 
    5860   ! 
    59    REAL(wp), PUBLIC, DIMENSION(:,:,:), ALLOCATABLE ::  pwcp       !: pore water sediment data at given time-step 
    60    REAL(wp), PUBLIC, DIMENSION(:,:,:), ALLOCATABLE ::  pwcp0      !: pore water sediment data at initial time 
    61    REAL(wp), PUBLIC, DIMENSION(:,:,:), ALLOCATABLE ::  solcp      !: solid sediment data at given time-step 
    62    REAL(wp), PUBLIC, DIMENSION(:,:,:), ALLOCATABLE ::  solcp0     !: solid sediment data at initial time 
    63    REAL(wp), PUBLIC, DIMENSION(:,:,:), ALLOCATABLE ::  diff 
     61   REAL(wp), PUBLIC, DIMENSION(:,:,:),  ALLOCATABLE ::  pwcp, pwcpa 
     62   REAL(wp), PUBLIC, DIMENSION(:,:,:),  ALLOCATABLE ::  solcp, solcpa 
     63   REAL(wp), PUBLIC, DIMENSION(:,:,:),  ALLOCATABLE ::  diff 
    6464 
    6565   !! * Shared module variables 
    6666   REAL(wp), PUBLIC, DIMENSION(:,:  ), ALLOCATABLE ::  pwcp_dta   !: pore water data at given time-step 
    6767   REAL(wp), PUBLIC, DIMENSION(:,:  ), ALLOCATABLE ::  rainrm_dta !: rain data at at initial time 
    68    REAL(wp), PUBLIC, DIMENSION(:,:  ), ALLOCATABLE ::  rainrm     !: rain data at given time-step 
    6968   REAL(wp), PUBLIC, DIMENSION(:,:  ), ALLOCATABLE ::  rainrg     !: rain of each solid component in [g/(cm**2.s)] 
    7069   REAL(wp), PUBLIC, DIMENSION(:,:  ), ALLOCATABLE ::  fromsed    !: 
    7170   REAL(wp), PUBLIC, DIMENSION(:,:  ), ALLOCATABLE ::  tosed      !: 
    72    REAL(wp), PUBLIC, DIMENSION(:,:  ), ALLOCATABLE ::  rloss      !:  
    73    REAL(wp), PUBLIC, DIMENSION(:,:  ), ALLOCATABLE ::  tokbot         
     71   REAL(wp), PUBLIC, DIMENSION(:,:  ), ALLOCATABLE ::  rearatpom  !:  
     72   REAL(wp), PUBLIC, DIMENSION(:,:  ), ALLOCATABLE ::  apluss, aminuss  !:  
    7473   ! 
    7574   REAL(wp), PUBLIC, DIMENSION(:    ), ALLOCATABLE ::  temp       !: temperature 
    7675   REAL(wp), PUBLIC, DIMENSION(:    ), ALLOCATABLE ::  salt       !: salinity 
    77    REAL(wp), PUBLIC, DIMENSION(:    ), ALLOCATABLE ::  press      !: pressure 
    78    REAL(wp), PUBLIC, DIMENSION(:    ), ALLOCATABLE ::  raintg     !: total massic flux rained in each cell (sum of sol. comp.) 
    7976   REAL(wp), PUBLIC, DIMENSION(:    ), ALLOCATABLE ::  fecratio   !: Fe/C ratio in falling particles to the sediments 
    80    REAL(wp), PUBLIC, DIMENSION(:    ), ALLOCATABLE ::  dzdep      !: total thickness of solid material rained [cm] in each cell 
     77   REAL(wp), PUBLIC, DIMENSION(:    ), ALLOCATABLE ::  dzdep, slatit, slongit   !: total thickness of solid material rained [cm] in each cell 
    8178   REAL(wp), PUBLIC, DIMENSION(:    ), ALLOCATABLE ::  zkbot      !: total thickness of solid material rained [cm] in each cell 
    8279   REAL(wp), PUBLIC, DIMENSION(:    ), ALLOCATABLE ::  wacc       !: total thickness of solid material rained [cm] in each cell 
     
    8784   REAL(wp), PUBLIC, DIMENSION(:,:  ), ALLOCATABLE ::  volw3d     !:  ??? 
    8885   REAL(wp), PUBLIC, DIMENSION(:,:  ), ALLOCATABLE ::  vols3d     !:  ??? 
     86   REAL(wp), PUBLIC, DIMENSION(:,:,:), ALLOCATABLE ::  volc 
     87   REAL(wp), PUBLIC, DIMENSION(:    ), ALLOCATABLE ::  dens_sol   !: Density of each solid fraction 
    8988 
    9089 
     
    104103   REAL(wp), PUBLIC, DIMENSION(:    ), ALLOCATABLE ::  ak123ps 
    105104   REAL(wp), PUBLIC, DIMENSION(:    ), ALLOCATABLE ::  aksis  
     105   REAL(wp), PUBLIC, DIMENSION(:    ), ALLOCATABLE ::  aknh3 
    106106   REAL(wp), PUBLIC, DIMENSION(:    ), ALLOCATABLE ::  aksps  
     107   REAL(wp), PUBLIC, DIMENSION(:    ), ALLOCATABLE ::  akh2s 
    107108   REAL(wp), PUBLIC, DIMENSION(:    ), ALLOCATABLE ::  sieqs 
    108109   REAL(wp), PUBLIC, DIMENSION(:    ), ALLOCATABLE ::  aks3s 
     
    117118   INTEGER , PUBLIC, SAVE                          ::  jpoce, indoce !: Ocean points ( number/indices ) 
    118119   INTEGER , PUBLIC, DIMENSION(:    ), ALLOCATABLE ::  iarroce       !: Computation of 1D array of sediments points 
     120   INTEGER , PUBLIC, DIMENSION(:, : ), ALLOCATABLE ::  jarr       !: Computation of 1D array of sediments points 
     121   INTEGER , PUBLIC, DIMENSION(:    ), ALLOCATABLE ::  jsvode, isvode   !: Computation of 1D array of sediments points 
     122 
    119123   REAL(wp), PUBLIC, DIMENSION(:,:  ), ALLOCATABLE ::  epkbot        !: ocean bottom layer thickness 
    120124   REAL(wp), PUBLIC, DIMENSION(:,:  ), ALLOCATABLE ::  gdepbot       !: Depth of the sediment 
     
    127131   REAL(wp), PUBLIC, DIMENSION(:,:  ), ALLOCATABLE ::  db            !: bioturbation ceofficient 
    128132   REAL(wp), PUBLIC, DIMENSION(:,:  ), ALLOCATABLE ::  irrig        !: bioturbation ceofficient 
     133   REAL(wp), PUBLIC, DIMENSION(:,:  ), ALLOCATABLE ::  radssol, rads1sol 
     134   REAL(wp), PUBLIC, DIMENSION(:,:  ), ALLOCATABLE ::  saturco3 
    129135   REAL(wp), PUBLIC, DIMENSION(:    ), ALLOCATABLE ::  rdtsed        !:  sediment model time-step 
    130    REAL(wp), PUBLIC, DIMENSION(:,:  ), ALLOCATABLE :: sedligand 
     136   REAL(wp), PUBLIC, DIMENSION(:    ), ALLOCATABLE ::  rstepros      !:  Number of iteration of rosenbrock method 
    131137   REAL(wp)  ::   dens               !: density of solid material 
    132138   !! Inputs / Outputs 
     
    138144   CHARACTER( len = 20 ), DIMENSION(jpdia2dsed) ::  seddia2d, seddia2u 
    139145   ! 
    140    REAL(wp), PUBLIC, DIMENSION(:,:,:,:), ALLOCATABLE ::  trcsedi 
    141    REAL(wp), PUBLIC, DIMENSION(:,:,:,:), ALLOCATABLE ::  flxsedi3d 
    142    REAL(wp), PUBLIC, DIMENSION(:,:,:  ), ALLOCATABLE ::  flxsedi2d 
    143146 
    144147   INTEGER, PUBLIC ::  numsed = 27    ! units 
     
    158161         &      dz(jpksed)  , por(jpksed) , por1(jpksed)                  ,   & 
    159162         &      volw(jpksed), vols(jpksed), rdtsed(jpksed)  ,   & 
    160          &      trcsedi  (jpi,jpj,jpksed,jptrased)                        ,   & 
    161          &      flxsedi3d(jpi,jpj,jpksed,jpdia3dsed)                      ,   & 
    162          &      flxsedi2d(jpi,jpj,jpdia2dsed)                             ,   & 
    163163         &      mol_wgt(jpsol),                                           STAT=sed_alloc ) 
    164164 
  • NEMO/trunk/src/TOP/PISCES/SED/sedadv.F90

    r10425 r15450  
    1515 
    1616   PUBLIC sed_adv 
    17    PUBLIC sed_adv_alloc 
    18  
    19    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:) :: dvolsp, dvolsm 
    20    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:) :: c2por, ckpor 
    2117 
    2218   REAL(wp) :: cpor 
     
    4945      ! * local variables 
    5046      INTEGER :: ji, jk, js  
    51       INTEGER :: jn, ntimes, nztime, ikwneg 
     47      INTEGER :: jn 
    5248       
    53       REAL(wp), DIMENSION(jpksed,jpsol) :: zsolcpno 
    54       REAL(wp), DIMENSION(jpksed)       :: zfilled, zfull, zfromup, zempty 
    55       REAL(wp), DIMENSION(jpoce,jpksed) :: zgap, zwb 
    56       REAL(wp), DIMENSION(jpoce,jpsol) :: zrainrf 
    57       REAL(wp), DIMENSION(:  ), ALLOCATABLE :: zraipush 
    58  
    59       REAL(wp) :: zkwnup, zkwnlo, zfrac,  zfromce, zrest, sumtot, zsumtot1 
     49      REAL(wp), DIMENSION(jpoce,jpksed,jpsol+jpads) :: zsolcp 
     50      REAL(wp) :: solfu, zfilled, zwb, fulsed, uebers, seddef 
    6051 
    6152      !------------------------------------------------------------------------ 
    62  
    6353 
    6454      IF( ln_timing )  CALL timing_start('sed_adv') 
     
    7060            WRITE(numsed,*) ' ' 
    7161         ENDIF 
    72          por1clay = denssol * por1(jpksed) * dz(jpksed) 
    73          cpor     = por1(jpksed) / por1(2) 
    74          DO jk = 2, jpksed 
    75             c2por(jk) = por1(2)      / por1(jk) 
    76             ckpor(jk) = por1(jpksed) / por1(jk) 
    77          ENDDO 
    78          DO jk = jpksedm1, 2, -1 
    79             dvolsp(jk) = vols(jk+1) / vols(jk) 
    80          ENDDO 
    81          DO jk = 3, jpksed 
    82            dvolsm(jk) = vols(jk-1) / vols(jk) 
    83         ENDDO 
    8462      ENDIF 
     63 
     64      ! Allocation of the temporary arrays 
     65      ! ---------------------------------- 
     66      zsolcp(:,:,:) = 0._wp 
     67      DO js = 1, jpsol 
     68         zsolcp(:,:,js) = solcp(:,:,js) 
     69      END DO 
     70      DO jk = 2, jpksed 
     71         zsolcp(:,jk,jpsol+1) = pwcp(:,jk,jwnh4) * adsnh4  
     72         zsolcp(:,jk,jpsol+2) = pwcp(:,jk,jwfe2) * adsfe2 
     73      END DO 
    8574 
    8675      ! Initialization of data for mass balance calculation 
     
    8877      fromsed(:,:) = 0. 
    8978      tosed  (:,:) = 0.  
    90       rloss  (:,:) = 0. 
    91       ikwneg = 1 
    92       nztime = jpksed 
    9379 
    94       ALLOCATE( zraipush(nztime) ) 
    95  
     80      solfu = 0.0 
     81      DO jk = 2, jpksed 
     82         solfu = solfu + dz(jk) * por1(jk) 
     83      END DO 
    9684      ! Initiate gap  
    9785      !-------------- 
    98       zgap(:,:) = 0. 
    99       DO js = 1, jpsol 
    100          DO jk = 1, jpksed 
    101             DO ji = 1, jpoce 
    102                zgap(ji,jk) = zgap(ji,jk) + solcp(ji,jk,js) 
    103             END DO 
    104          ENDDO 
    105       ENDDO 
    106  
    107       zgap(1:jpoce,1:jpksed) = 1. - zgap(1:jpoce,1:jpksed)    
    10886 
    10987      ! Initiate burial rates 
    11088      !----------------------- 
    111       zwb(:,:) = 0. 
    112       DO jk = 2, jpksed 
    113          zfrac =  dtsed / ( denssol * por1(jk) )      
    114          DO ji = 1, jpoce 
    115             zwb(ji,jk) = zfrac * raintg(ji) 
    116          END DO 
    117       ENDDO 
     89      DO ji = 1, jpoce 
     90         DO jk = 2, jpksed-1 
     91            zfilled = 0._wp 
     92            DO js = 1, jpsol 
     93               zfilled = zfilled + zsolcp(ji,jk,js) / dens_sol(js) 
     94            END DO 
     95            zwb = MAX(0._wp, (zfilled - 1._wp) / (zfilled + rtrn) ) 
    11896 
    11997 
    120       DO ji = 1, jpoce 
    121          zwb(ji,2) = zwb(ji,2) - zgap(ji,2) * dz(2) 
    122       ENDDO 
     98            DO js = 1, jpsol + jpads 
     99               uebers = zwb * zsolcp(ji,jk,js) 
     100               zsolcp(ji,jk,js) = zsolcp(ji,jk,js) - uebers 
     101               zsolcp(ji,jk+1,js) = zsolcp(ji,jk+1,js) + uebers * dz(jk) * por1(jk) / ( dz(jk+1) * por1(jk+1) ) 
     102            END DO 
     103         END DO 
    123104 
    124       DO jk = 3, jpksed 
    125          zfrac = por1(jk-1) / por1(jk) 
    126          DO ji = 1, jpoce 
    127             zwb(ji,jk) = zwb(ji,jk-1) * zfrac - zgap(ji,jk) * dz(jk) 
     105         zfilled = 0._wp 
     106         DO js = 1, jpsol 
     107            zfilled = zfilled + zsolcp(ji,jpksed,js) / dens_sol(js) 
    128108         END DO 
    129       ENDDO 
    130  
    131       zrainrf(:,:) = 0. 
    132       DO ji = 1, jpoce 
    133          IF( raintg(ji) /= 0. )  & 
    134             &   zrainrf(ji,:) = rainrg(ji,:) / raintg(ji) 
    135       ENDDO 
    136  
    137  
    138       ! Computation of full and empty solid fraction in each layer 
    139       ! for all 'burial' case 
    140       !---------------------------------------------------------- 
    141  
     109         zwb = MAX(0._wp, (zfilled - 1._wp) / (zfilled + rtrn) ) 
     110         DO js = 1, jpsol + jpads 
     111            uebers = zwb * zsolcp(ji,jpksed,js) 
     112            zsolcp(ji,jpksed,js) = zsolcp(ji,jpksed,js) - uebers 
     113            tosed(ji,js) = uebers * dz(jpksed) * por1(jpksed) 
     114         END DO 
     115      END DO 
    142116 
    143117      DO ji = 1, jpoce 
     118         fulsed = 0._wp 
     119         DO jk = 2, jpksed  
     120            zfilled = 0._wp 
     121            DO js = 1, jpsol 
     122               zfilled = zfilled + zsolcp(ji,jk,js) / dens_sol(js) 
     123            END DO 
     124            fulsed = fulsed + zfilled * dz(jk) * por1(jk) 
     125         END DO  
    144126 
    145          ! computation of total weight fraction in sediment 
    146          !------------------------------------------------- 
    147          zfilled(:) = 0. 
    148          DO js = 1, jpsol 
    149             DO jk = 2, jpksed 
    150                zfilled(jk) = zfilled(jk) + solcp(ji,jk,js) 
    151             ENDDO 
    152          ENDDO 
    153           
    154          DO js = 1, jpsol 
    155             DO jk = 2, jpksed 
    156                zsolcpno(jk,js) = solcp(ji,jk,js) / zfilled(jk) 
    157             ENDDO 
    158          ENDDO 
     127         seddef = solfu - fulsed 
    159128 
    160          ! burial  3 cases:  
    161          ! zwb > 0 ==> rain > total rection loss  
    162          ! zwb = 0 ==> rain = 0 
    163          ! zwb < 0 ==> rain > 0 and rain < total reaction loss 
    164          !---------------------------------------------------------------- 
     129         zsolcp(ji,jpksed,jsclay) = zsolcp(ji,jpksed,jsclay) + seddef / ( por1(jpksed) * dz(jpksed) )    & 
     130         &         * dens_sol(jsclay) 
     131         fromsed(ji,jsclay) = seddef * dens_sol(jsclay)  
    165132 
    166          IF( zwb(ji,jpksed) > 0. ) THEN 
     133         DO jk = jpksed, 3, -1 
     134            zfilled = 0._wp 
     135            DO js = 1, jpsol 
     136               zfilled = zfilled + zsolcp(ji,jk,js) / dens_sol(js) 
     137            END DO 
     138            zwb = MAX(0._wp, (zfilled - 1._wp) / (zfilled + rtrn) ) 
     139            DO js = 1, jpsol + jpads 
     140               uebers = zwb * zsolcp(ji,jk,js) 
     141               zsolcp(ji,jk,js) = zsolcp(ji,jk,js) - uebers 
     142               zsolcp(ji,jk-1,js) = zsolcp(ji,jk-1,js) + uebers * dz(jk) * por1(jk) / ( dz(jk-1) * por1(jk-1) ) 
     143            END DO 
     144         END DO 
    167145 
    168             zfull (jpksed) = zfilled(jpksed) 
    169             zempty(jpksed) = 1. - zfull(jpksed) 
    170             DO jk = jpksedm1, 2, -1 
    171                zfull (jk) = zfilled(jk) 
    172                zfull (jk) = zfull(jk) - zempty(jk+1) * dvolsp(jk) 
    173                zempty(jk) = 1. - zfull(jk) 
    174             ENDDO 
     146      END DO 
    175147 
    176             ! Computation of solid sediment species 
    177             !-------------------------------------- 
    178             ! push entire sediment column downward to account rest of rain 
    179             DO js = 1, jpsol 
    180                DO jk = jpksed, 3, -1 
    181                   solcp(ji,jk,js) = zfull(jk) * zsolcpno(jk,js) + zempty(jk) * zsolcpno(jk-1,js) 
    182                ENDDO 
     148      DO js = 1, jpsol 
     149         solcp(:,:,js) = zsolcp(:,:,js) 
     150      END DO 
     151      DO jk = 2, jpksed 
     152         pwcp(:,jk,jwnh4) = (pwcp(:,jk,jwnh4) + zsolcp(:,jk,jpsol+1) * por1(jk) / por(jk) ) * radssol(jk,jwnh4) 
     153         IF (jpoce == 146 .and. slatit(145) > 0.) write(0,*) 'plante advection ',pwcp(145,jk,jwfe2)*1E6,zsolcp(145,jk,jpsol+2)*1E6,    & 
     154         &         (pwcp(145,jk,jwfe2) + zsolcp(145,jk,jpsol+2) * por1(jk) / por(jk) ) * radssol(jk,jwfe2) * 1E6 
     155         pwcp(:,jk,jwfe2) = (pwcp(:,jk,jwfe2) + zsolcp(:,jk,jpsol+2) * por1(jk) / por(jk) ) * radssol(jk,jwfe2) 
     156      END DO 
    183157 
    184                solcp(ji,2,js) = zfull(2) * zsolcpno(2,js) + zempty(2) * zrainrf(ji,js) 
    185  
    186                DO jk = 2, jpksed 
    187                   zsolcpno(jk,js) = solcp(ji,jk,js) 
    188                END DO 
    189             ENDDO 
    190  
    191             zrest = zwb(ji,jpksed) * cpor 
    192             ! what is remaining is less than dz(2) 
    193             IF( zrest <= dz(2) ) THEN            
    194  
    195                zfromup(2) = zrest / dz(2) 
    196                DO jk = 3, jpksed 
    197                   zfromup(jk) = zwb(ji,jpksed) * ckpor(jk) / dz(jk) 
    198                ENDDO 
    199                DO js = 1, jpsol 
    200                   zfromce = 1. - zfromup(2) 
    201                   solcp(ji,2,js) = zfromce * zsolcpno(2,js) + zfromup(2) * zrainrf(ji,js) 
    202                   DO jk = 3, jpksed 
    203                      zfromce = 1. - zfromup(jk) 
    204                      solcp(ji,jk,js) = zfromce * zsolcpno(jk,js) + zfromup(jk) * zsolcpno(jk-1,js) 
    205                   ENDDO 
    206                   fromsed(ji,js) = 0. 
    207                   ! quantities to push in deeper sediment 
    208                   tosed  (ji,js) = zsolcpno(jpksed,js) & 
    209                      &           * zwb(ji,jpksed) * denssol * por1(jpksed) 
    210                ENDDO 
    211  
    212             ELSE ! what is remaining is great than dz(2) 
    213  
    214                ntimes = INT( zrest / dz(2) ) + 1 
    215                IF( ntimes > nztime ) CALL ctl_stop( 'STOP', 'sed_adv : rest too large ' ) 
    216                zraipush(1) = dz(2) 
    217                zrest = zrest - zraipush(1) 
    218                DO jn = 2, ntimes 
    219                   IF( zrest >= dz(2) ) THEN 
    220                      zraipush(jn) = dz(2) 
    221                      zrest = zrest - zraipush(jn) 
    222                   ELSE 
    223                      zraipush(jn) = zrest 
    224                      zrest = 0. 
    225                   ENDIF 
    226                ENDDO 
    227  
    228                DO jn = 1, ntimes 
    229                   DO js = 1, jpsol 
    230                      DO jk = 2, jpksed 
    231                         zsolcpno(jk,js) = solcp(ji,jk,js) 
    232                      END DO 
    233                   ENDDO 
    234                    
    235                   zfromup(2) = zraipush(jn) / dz(2) 
    236                   DO jk = 3, jpksed 
    237                      zfromup(jk) = ( zraipush(jn) / dz(jk) ) * c2por(jk) 
    238                   ENDDO 
    239  
    240                   DO js = 1, jpsol 
    241                      zfromce = 1. - zfromup(2) 
    242                      solcp(ji,2,js) = zfromce * zsolcpno(2,js) + zfromup(2) * zrainrf(ji,js) 
    243                      DO jk = 3, jpksed 
    244                         zfromce = 1. - zfromup(jk) 
    245                         solcp(ji,jk,js) = zfromce * zsolcpno(jk,js) + zfromup(jk) * zsolcpno(jk-1,js) 
    246                      ENDDO 
    247                      fromsed(ji,js) = 0. 
    248                      tosed  (ji,js) = tosed(ji,js) + zsolcpno(jpksed,js) * zraipush(jn) & 
    249                         &             * denssol * por1(2)  
    250                   ENDDO 
    251                ENDDO 
    252   
    253             ENDIF 
    254  
    255          ELSE IF( raintg(ji) < eps ) THEN ! rain  = 0 
    256 !! Nadia    rloss(:,:) = rainrm(:,:)   bug ??????            
    257  
    258             rloss(ji,1:jpsol) = rainrm(ji,1:jpsol) 
    259  
    260             zfull (2) = zfilled(2) 
    261             zempty(2) = 1. - zfull(2) 
    262             DO jk = 3, jpksed 
    263                zfull (jk) = zfilled(jk) 
    264                zfull (jk) = zfull (jk) - zempty(jk-1) * dvolsm(jk) 
    265                zempty(jk) = 1. - zfull(jk) 
    266             ENDDO 
    267  
    268             ! fill boxes with weight fraction from underlying box 
    269             DO js = 1, jpsol 
    270                DO jk = 2, jpksedm1 
    271                   solcp(ji,jk,js) = zfull(jk) * zsolcpno(jk,js) + zempty(jk) * zsolcpno(jk+1,js) 
    272                END DO 
    273                solcp(ji,jpksed,js) = zsolcpno(jpksed,js) * zfull(jpksed) 
    274                tosed  (ji,js) = 0. 
    275                fromsed(ji,js) = 0. 
    276             ENDDO 
    277             ! for the last layer, one make go up clay  
    278             solcp(ji,jpksed,jsclay) = solcp(ji,jpksed,jsclay) + zempty(jpksed) * 1. 
    279             fromsed(ji,jsclay) = zempty(jpksed) * 1. * por1clay 
    280          ELSE  ! rain > 0 and rain < total reaction loss 
    281  
    282  
    283             DO jk = 2, jpksed 
    284                zfull (jk) = zfilled(jk) 
    285                zempty(jk) = 1. - zfull(jk) 
    286             ENDDO 
    287  
    288             ! Determination of indice of layer - ikwneg - where advection is reversed 
    289             !------------------------------------------------------------------------ 
    290  iflag:     DO jk = 2, jpksed 
    291                IF( zwb(ji,jk) < 0.  ) THEN 
    292                   ikwneg = jk 
    293                   EXIT iflag 
    294                ENDIF 
    295             ENDDO iflag 
    296  
    297             ! computation of zfull and zempty  
    298             ! 3 cases : a/ ikwneg=2, b/ikwneg=3...jpksedm1, c/ikwneg=jpksed     
    299             !-------------------------------------------------------------       
    300             IF( ikwneg == 2 ) THEN ! advection is reversed in the first sediment layer 
    301  
    302                zkwnup = rdtsed(ikwneg) * raintg(ji) / dz(ikwneg) 
    303                zkwnlo = ABS( zwb(ji,ikwneg) ) / dz(ikwneg) 
    304                zfull (ikwneg+1) = zfilled(ikwneg+1) - zkwnlo * dvolsm(ikwneg+1) 
    305                zempty(ikwneg+1) = 1. - zfull(ikwneg+1) 
    306                DO jk = ikwneg+2, jpksed 
    307                   zfull (jk) = zfilled(jk) - zempty(jk-1) * dvolsm(jk) 
    308                   zempty(jk) = 1. - zfull(jk) 
    309                ENDDO 
    310                DO js = 1, jpsol 
    311                   solcp(ji,2,js) = zfull(2) * zsolcpno(2,js)+ zkwnlo * zsolcpno(3,js) & 
    312                      &                                      + zkwnup * zrainrf(ji,js) 
    313                   DO jk = 3, jpksedm1 
    314                      solcp(ji,jk,js) = zfull(jk) * zsolcpno(jk,js) + zempty(jk) * zsolcpno(jk+1,js) 
    315                   ENDDO 
    316                   solcp(ji,jpksed,js) = zfull(jpksed) * zsolcpno(jpksed,js) 
    317                   tosed(ji,js)   = 0. 
    318                   fromsed(ji,js) = 0. 
    319                ENDDO 
    320                solcp(ji,jpksed,jsclay) =  solcp(ji,jpksed,jsclay) + zempty(jpksed) * 1. 
    321                !! C. Heinze  fromsed(ji,jsclay) = zempty(jpksed) * 1. * denssol * por1(jpksed) / mol_wgt(jsclay) 
    322                fromsed(ji,jsclay) = zempty(jpksed) * 1. * por1clay 
    323                 
    324             ELSE IF( ikwneg == jpksed ) THEN 
    325  
    326                zkwnup = ABS( zwb(ji,ikwneg-1) ) * dvolsm(ikwneg) / dz(ikwneg) 
    327                zkwnlo = ABS( zwb(ji,ikwneg) ) / dz(ikwneg) 
    328                zfull (ikwneg-1) = zfilled(ikwneg-1) - zkwnup * dvolsp(ikwneg-1) 
    329                zempty(ikwneg-1) = 1. - zfull(ikwneg-1) 
    330                DO jk = ikwneg-2, 2, -1 
    331                   zfull (jk) = zfilled(jk) - zempty(jk+1) * dvolsp(jk) 
    332                   zempty(jk) = 1. - zfull(jk)  
    333                ENDDO 
    334                DO  js = 1, jpsol 
    335                   solcp(ji,2,js) = zfull(2) * zsolcpno(2,js) + zempty(2) * zrainrf(ji,js) 
    336                ENDDO 
    337                DO  js = 1, jpsol 
    338                   DO jk = jpksedm1, 3, -1 
    339                      solcp(ji,jk,js) = zfull(jk) * zsolcpno(jk,js) + zempty(jk) * zsolcpno(jk-1,js) 
    340                   ENDDO 
    341                   solcp(ji,jpksed,js) = zfull(jpksed) * zsolcpno(jpksed,js) & 
    342                      &                       + zkwnup * zsolcpno(jpksedm1,js) 
    343                   tosed(ji,js)   = 0. 
    344                   fromsed(ji,js) = 0. 
    345                ENDDO 
    346                solcp(ji,jpksed,jsclay) = solcp(ji,jpksed,jsclay) + zkwnlo * 1. 
    347                ! Heinze  fromsed(ji,jsclay) = zkwnlo * 1. * denssol * por1(jpksed) / mol_wgt(jsclay) 
    348                fromsed(ji,jsclay) = zkwnlo * 1.* por1clay 
    349             ELSE   ! 2 < ikwneg(ji) <= jpksedm1 
    350  
    351                zkwnup = ABS( zwb(ji,ikwneg-1) ) * por1(ikwneg-1) / ( dz(ikwneg) * por1(ikwneg) ) 
    352                zkwnlo = ABS( zwb(ji,ikwneg) ) / dz(ikwneg) 
    353  
    354                IF( ikwneg > 3 ) THEN 
    355  
    356                   zfull (ikwneg-1) = zfilled(ikwneg-1) - zkwnup * dvolsp(ikwneg-1) 
    357                   zempty(ikwneg-1) = 1. - zfull(ikwneg-1)  
    358                   DO jk = ikwneg-2, 2, -1 
    359                      zfull (jk) = zfilled(jk) - zempty(jk+1) * dvolsp(jk) 
    360                      zempty(jk)    = 1. - zfull(jk)  
    361                   ENDDO 
    362                   DO js = 1, jpsol 
    363                      solcp(ji,2,js) = zfull(2) * zsolcpno(2,js) + zempty(2) * zrainrf(ji,js) 
    364                   ENDDO 
    365                   DO js = 1, jpsol 
    366                      DO jk = ikwneg-1, 3, -1 
    367                         solcp(ji,jk,js) = zfull(jk) * zsolcpno(jk,js) + zempty(jk) * zsolcpno(jk-1,js) 
    368                      ENDDO 
    369                   ENDDO 
    370                ELSE ! ikw = 3 
    371  
    372  
    373                   zfull (2) = zfilled(2) - zkwnup * dvolsm(3) 
    374                   zempty(2) = 1. - zfull(2) 
    375                   DO js = 1, jpsol 
    376                      solcp(ji,2,js) = zfull(2) * zsolcpno(2,js) + zempty(2) * zrainrf(ji,js) 
    377                   ENDDO 
    378                ENDIF 
    379  
    380                IF( ikwneg < jpksedm1) THEN 
    381  
    382                   zfull (ikwneg+1) = zfilled(ikwneg+1) - zkwnlo * dvolsm(ikwneg+1) 
    383                   zempty(ikwneg+1) = 1. - zfull(ikwneg+1)  
    384                   DO jk = ikwneg+2, jpksed 
    385                      zfull (jk) = zfilled(jk) - zempty(jk-1) * dvolsm(jk) 
    386                      zempty(jk) = 1. - zfull(jk) 
    387                   ENDDO 
    388                   DO js = 1, jpsol 
    389                      DO jk = ikwneg+1, jpksedm1 
    390                         solcp(ji,jk,js) = zfull(jk) * zsolcpno(jk,js) + zempty(jk) * zsolcpno(jk+1,js)  
    391                      ENDDO 
    392                      solcp(ji,jpksed,js) = zfull(jpksed) * zsolcpno(jpksed,js) 
    393                   ENDDO 
    394                   solcp(ji,jpksed,jsclay) = solcp(ji,jpksed,jsclay) + zempty(jpksed) * 1. 
    395                ELSE 
    396  
    397                   zfull (jpksed) = zfilled(jpksed) - zkwnlo * dvolsm(jpksed) 
    398                   zempty(jpksed) = 1. - zfull(jpksed)                 
    399                   DO js = 1, jpsol 
    400                      solcp(ji,jpksed,js) = zfull(jpksed) * zsolcpno(jpksed,js) 
    401                   ENDDO 
    402                   solcp(ji,jpksed,jsclay) = solcp(ji,jpksed,jsclay) + zempty(jpksed) * 1. 
    403                ENDIF ! jpksedm1 
    404                 
    405                ! ikwneg = jpksedm1  ; ikwneg+1 = jpksed ; ikwneg-1 = jpksed - 2 
    406                DO js = 1, jpsol 
    407                   solcp(ji,ikwneg,js) =  zfull(ikwneg) * zsolcpno(ikwneg  ,js) & 
    408                      &                +  zkwnup        * zsolcpno(ikwneg-1,js) & 
    409                      &                +  zkwnlo        * zsolcpno(ikwneg+1,js)    
    410                   tosed  (ji,js)   = 0. 
    411                   fromsed(ji,js)   = 0. 
    412                ENDDO 
    413                ! Heinze  fromsed(ji,jsclay) = zempty * 1. * denssol * por1(jpksed) / mol_wgt(jsclay) 
    414                fromsed(ji,jsclay) = zempty(jpksed) * 1. * por1clay 
    415  
    416             ENDIF ! ikwneg(ji) = 2 
    417          ENDIF  ! zwb > 0 
    418       ENDDO  ! ji = 1, jpoce 
    419  
    420       rainrm(:,:) = 0. 
    421158      rainrg(:,:) = 0. 
    422       raintg(:)   = 0. 
    423  
    424       DEALLOCATE( zraipush ) 
    425159 
    426160      IF( ln_timing )  CALL timing_stop('sed_adv') 
     
    429163 
    430164 
    431    INTEGER FUNCTION sed_adv_alloc() 
    432       !!---------------------------------------------------------------------- 
    433       !!                     ***  ROUTINE p4z_prod_alloc  *** 
    434       !!---------------------------------------------------------------------- 
    435       ALLOCATE( dvolsp(jpksed), dvolsm(jpksed), c2por(jpksed),         & 
    436       &         ckpor(jpksed) ,           STAT = sed_adv_alloc ) 
    437       ! 
    438       IF( sed_adv_alloc /= 0 ) CALL ctl_stop( 'STOP', 'sed_adv_alloc : failed to allocate arrays.' ) 
    439       ! 
    440    END FUNCTION sed_adv_alloc 
    441  
    442165END MODULE sedadv 
  • NEMO/trunk/src/TOP/PISCES/SED/sedarr.F90

    r14086 r15450  
    1010   !!---------------------------------------------------------------------- 
    1111   !! * Modules used 
    12    USE par_oce 
    1312   USE par_sed 
    14    USE in_out_manager, ONLY : ln_timing 
    15    USE timing 
     13   USE dom_oce 
     14   USE sed 
    1615 
    1716   IMPLICIT NONE 
     
    5251         jjd        = ( tab_ind(jn) - 1 ) / jpi + 1 
    5352         tab1d(jn)  = tab2d(jid, jjd) 
     53!         IF (mig(jid) == 150 .and. mjg(jjd) == 136) write(0,*) 'plante indices ',jn,ndim1d,slatit(jn),slongit(jn)  
    5454      END DO  
    5555 
  • NEMO/trunk/src/TOP/PISCES/SED/sedbtb.F90

    r10222 r15450  
    66   !! * Modules used 
    77   USE sed     ! sediment global variable 
     8   USE sed_oce 
    89   USE sedmat  ! linear system of equations 
     10   USE sedinorg 
     11   USE sedorg 
     12   USE sedini 
    913   USE lib_mpp         ! distribued memory computing library 
    1014 
     
    3135      !!        !  06-04 (C. Ethe)  Re-organization 
    3236      !!---------------------------------------------------------------------- 
    33       !!* Arguments 
    34       INTEGER, INTENT(in) ::  kt              ! time step 
    35  
     37      !! * Arguments 
     38      INTEGER, INTENT(in)  :: kt   ! time step 
    3639      ! * local variables 
    3740      INTEGER :: ji, jk, js 
    38       REAL(wp), DIMENSION(jpoce,jpksedm1,jpsol) ::  zsol  !   solution 
     41      REAL(wp), DIMENSION(jpoce,jpksed,jpsol) ::  zrearat  !   solution 
     42      REAL(wp) :: zsolid1, zsolid2, zsolid3, zsumtot, zlimo2 
    3943      !------------------------------------------------------------------------ 
    4044 
     
    4246 
    4347      IF( kt == nitsed000 ) THEN 
    44          IF (lwp) WRITE(numsed,*) ' sed_btb : Bioturbation  ' 
     48         IF (lwp) WRITE(numsed,*) ' sed_btb : bioturbation of solid and adsorbed species  ' 
    4549         IF (lwp) WRITE(numsed,*) ' ' 
    4650      ENDIF 
    4751 
     52 
    4853      ! Initializations 
    4954      !---------------- 
    50       zsol(:,:,:) = 0. 
     55      zrearat = 0. 
     56 
     57      ! Remineralization rates of the different POC pools 
     58      zrearat(:,:,jspoc) = -reac_pocl 
     59      zrearat(:,:,jspos) = -reac_pocs 
     60      zrearat(:,:,jspor) = -reac_pocr 
    5161 
    5262      ! right hand side of coefficient matrix 
    5363      !-------------------------------------- 
    54       DO js = 1, jpsol 
    55          DO jk = 1, jpksedm1 
    56             DO ji = 1, jpoce 
    57                zsol(ji,jk,js) = solcp(ji,jk+1,js) 
    58             ENDDO 
    59          ENDDO 
    60       ENDDO 
     64      CALL sed_mat_btbi( jpksed, jpsol, solcp, zrearat(:,:,:), dtsed ) 
    6165 
    62       CALL sed_mat( jpsol, jpoce, jpksedm1, zsol, dtsed / 2.0 ) 
     66      DO ji = 1, jpoce 
    6367 
     68         zsumtot = 0. 
     69         DO jk = 2, jpksed 
     70            zsolid1 = volc(ji,jk,jspoc) * solcp(ji,jk,jspoc) 
     71            zsolid2 = volc(ji,jk,jspos) * solcp(ji,jk,jspos) 
     72            zsolid3 = volc(ji,jk,jspor) * solcp(ji,jk,jspor) 
     73            rearatpom(ji,jk)  = ( reac_pocl * zsolid1 + reac_pocs * zsolid2 + reac_pocr * zsolid3 ) 
     74            zsumtot = zsumtot + rearatpom(ji,jk) * volw3d(ji,jk) * 1.e-3 * 86400. * 365. * 1E3 
     75         END DO 
    6476 
    65       ! store solution of the tridiagonal system 
    66       !------------------------ 
    67       DO js = 1, jpsol 
    68          DO jk = 1, jpksedm1 
    69             DO ji = 1, jpoce 
    70                solcp(ji,jk+1,js) = zsol(ji,jk,js) 
    71             ENDDO 
    72          ENDDO 
    73       ENDDO 
     77         !    4/ Computation of the bioturbation coefficient 
     78         !       This parameterization is taken from Archer et al. (2002) 
     79         ! -------------------------------------------------------------- 
     80         zlimo2   = max(0.01, pwcp(ji,1,jwoxy) / (pwcp(ji,1,jwoxy) + 20.E-6) ) 
     81         db(ji,:) = dbiot * zsumtot**0.85 * zlimo2 / (365.0 * 86400.0) 
     82 
     83         ! ------------------------------------------------------ 
     84         !    Vertical variations of the bioturbation coefficient 
     85         ! ------------------------------------------------------ 
     86         IF (ln_btbz) THEN 
     87            DO jk = 1, jpksed 
     88                  db(ji,jk) = db(ji,jk) * exp( -(profsedw(jk) / dbtbzsc)**2 ) 
     89            END DO 
     90         ELSE 
     91            DO jk = 1, jpksed 
     92               IF (profsedw(jk) > dbtbzsc) THEN 
     93                  db(ji,jk) = 0.0 
     94               ENDIF 
     95            END DO 
     96         ENDIF 
     97 
     98         ! Computation of the bioirrigation factor (from Archer, MUDS model) 
     99         irrig(ji,:) = 0.0 
     100         IF (ln_irrig) THEN 
     101            DO jk = 1, jpksed 
     102               irrig(ji,jk) = ( 7.63752 - 7.4465 * exp( -0.89603 * zsumtot ) ) * zlimo2 
     103               irrig(ji,jk) = irrig(ji,jk) * exp( -(profsedw(jk) / xirrzsc) ) 
     104            END DO 
     105         ENDIF 
     106      END DO 
     107  
     108      ! CALL inorganic and organic slow redow/chemistry processes 
     109      ! --------------------------------------------------------- 
     110      CALL sed_inorg( kt ) 
    74111 
    75112      IF( ln_timing )  CALL timing_stop('sed_btb') 
  • NEMO/trunk/src/TOP/PISCES/SED/sedchem.F90

    r15237 r15450  
    139139         CALL sed_chem_cst 
    140140      ELSE 
    141          DO_2D( 1, 1, 1, 1 ) 
     141         DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) 
    142142            ikt = mbkt(ji,jj)  
    143143            IF ( tmask(ji,jj,ikt) == 1 ) THEN 
  • NEMO/trunk/src/TOP/PISCES/SED/seddiff.F90

    r10225 r15450  
    77   !! * Modules used 
    88   USE sed     ! sediment global variable 
    9    USE sed_oce 
    109   USE sedmat  ! linear system of equations 
    1110   USE sedini 
     
    4241      !! Arguments 
    4342      INTEGER, INTENT(in) ::   kt, knt       ! number of iteration 
    44       ! --- local variables 
    45       INTEGER :: ji, jk, js   ! dummy looop indices 
     43!       
    4644 
    47       REAL(wp), DIMENSION(jpoce,jpksed) :: zrearat1, zrearat2   ! reaction rate in pore water 
    48       !! 
    49       !!---------------------------------------------------------------------- 
    50  
    51       IF( ln_timing )  CALL timing_start('sed_diff') 
    52 ! 
    53       IF( kt == nitsed000 .AND. knt == 1 ) THEN 
    54          IF (lwp) THEN 
    55             WRITE(numsed,*) ' sed_diff : pore-water diffusion ' 
    56             WRITE(numsed,*) ' ' 
    57          ENDIF 
    58       ENDIF 
    59  
    60      ! Initializations 
    61      !---------------------- 
    62       zrearat1(:,:)   = 0. 
    63       zrearat2(:,:) = 0. 
    64  
    65       !--------------------------- 
    66       ! Solves PO4 diffusion  
    67       !---------------------------- 
    68  
    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 ) 
    113  
    114       CALL sed_mat( jwdic, jpoce, jpksed, zrearat1, zrearat2, sedligand(:,:), dtsed2 / 2.0 ) 
    115  
    116       IF( ln_timing )  CALL timing_stop('sed_diff') 
    117 !       
    11845   END SUBROUTINE sed_diff 
    11946 
  • NEMO/trunk/src/TOP/PISCES/SED/seddsr.F90

    r10362 r15450  
    1919   !! * Module variables 
    2020 
    21    REAL(wp) :: zadsnh4 
    2221   REAL(wp), DIMENSION(jpsol), PUBLIC      :: dens_mol_wgt  ! molecular density  
    23    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: zvolc    ! temp. variables 
    24  
    2522 
    2623   !! $Id$ 
    2724CONTAINS 
    2825    
    29    SUBROUTINE sed_dsr( kt, knt )  
     26   SUBROUTINE sed_dsr( accmask ) 
    3027      !!---------------------------------------------------------------------- 
    3128      !!                   ***  ROUTINE sed_dsr  *** 
     
    4946      !!---------------------------------------------------------------------- 
    5047      !! Arguments 
    51       INTEGER, INTENT(in) ::   kt, knt       ! number of iteration 
    5248      ! --- local variables 
     49      INTEGER, DIMENSION(jpoce), INTENT(in) :: accmask 
    5350      INTEGER :: ji, jk, js, jw, jn   ! dummy looop indices 
    5451 
    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 
    59       REAL(wp)  ::  zsolid1, zsolid2, zsolid3, zvolw, zreasat 
    60       REAL(wp)  ::  zsatur, zsatur2, znusil, zkpoca, zkpocb, zkpocc 
    61       REAL(wp)  ::  zratio, zgamma, zbeta, zlimtmp, zundsat2 
     52      REAL(wp), DIMENSION(jpoce,jpksed) ::zlimo2, zlimno3, zlimso4, zlimfeo    ! undersaturation ; indice jpwatp1 is for calcite    
     53      REAL(wp)  ::  zreasat 
    6254      !! 
    6355      !!---------------------------------------------------------------------- 
     
    6557      IF( ln_timing )  CALL timing_start('sed_dsr') 
    6658! 
    67       IF( kt == nitsed000 .AND. knt == 1 ) THEN 
    68          IF (lwp) THEN 
    69             WRITE(numsed,*) ' sed_dsr : Dissolution reaction ' 
    70             WRITE(numsed,*) ' ' 
    71          ENDIF 
    72       ENDIF 
    73  
    7459     ! Initializations 
    7560     !---------------------- 
    7661       
    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 
     62      zlimo2 (:,:)    = 0.    ;   zlimno3(:,:)  = 0. 
     63      zlimso4(:,:)    = 0. 
    8264   
    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 
    96  
    97       ! Conversion of volume units 
    98       !---------------------------- 
    99       DO js = 1, jpsol 
    100          DO jk = 1, jpksed 
    101             DO ji = 1, jpoce 
    102                zvolc(ji,jk,js) = ( vols3d(ji,jk) * dens_mol_wgt(js) ) /  & 
    103                   &              ( volw3d(ji,jk) * 1.e-3 ) 
    104             ENDDO 
    105          ENDDO 
    106       ENDDO 
    107  
    10865      !---------------------------------------------------------- 
    10966      ! 5.  Beginning of solid reaction 
    11067      !--------------------------------------------------------- 
    11168       
    112       ! Definition of reaction rates [rearat]=sans dim  
    113       ! For jk=1 no reaction (pure water without solid) for each solid compo 
    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 
    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          
    186             ! For DIC 
    187             pwcp(ji,jk,jwdic)  = pwcp(ji,jk,jwdic) + zreasat 
    188             zsumtot(ji) = zsumtot(ji) + zreasat / dtsed2 * volw3d(ji,jk) * 1.e-3 * 86400. * 365. * 1E3 
    189             ! For Phosphate (in mol/l) 
    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 
    193             ! For alkalinity 
    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 
     69      ! Computes the SMS of the species which do not affect the remin 
     70      ! processes (Alkalinity, PO4, NH4, and the release of Fe from  
     71      ! biogenic Fe 
     72      ! In the following, all loops start from jpk = 2 as reactions  
     73      ! only occur in the sediment 
     74      ! -------------------------------------------------------------- 
     75      DO jk = 2, jpksed 
     76         DO ji = 1, jpoce 
     77            IF (accmask(ji) == 0 ) THEN 
     78               zreasat = rearatpom(ji,jk) 
     79               ! For alkalinity 
     80               pwcpa(ji,jk,jwalk)  = pwcpa(ji,jk,jwalk) + zreasat * ( srno3 - 2.* spo4r ) 
     81               ! For Phosphate (in mol/l) 
     82               pwcpa(ji,jk,jwpo4)  = pwcpa(ji,jk,jwpo4) + zreasat * spo4r 
     83               ! For Ammonium 
     84               pwcpa(ji,jk,jwnh4)  = pwcpa(ji,jk,jwnh4) + zreasat * srno3 * radssol(jk,jwnh4) 
     85               ! For iron (in mol/l) 
     86               pwcpa(ji,jk,jwfe2)  = pwcpa(ji,jk,jwfe2) + fecratio(ji) * zreasat * radssol(jk,jwfe2) 
     87            ENDIF 
     88         END DO 
     89      ENDDO 
     90 
     91      ! Computes SMS of oxygen 
     92      DO jk = 2, jpksed 
     93         DO ji = 1, jpoce 
     94            IF (accmask(ji) == 0 ) THEN 
     95               zlimo2(ji,jk) = pwcp(ji,jk,jwoxy) / ( pwcp(ji,jk,jwoxy) + xksedo2 ) 
     96               zlimo2(ji,jk) = MAX( 0., zlimo2(ji,jk) ) 
     97               zreasat = rearatpom(ji,jk) * zlimo2(ji,jk) 
     98               ! Acid Silicic  
     99               pwcpa(ji,jk,jwoxy)  = pwcpa(ji,jk,jwoxy) - so2ut * zreasat 
     100            ENDIF 
     101         END DO 
     102      ENDDO 
     103 
     104      !-------------------------------------------------------------------- 
     105      ! Denitrification 
     106      !-------------------------------------------------------------------- 
     107      DO jk = 2, jpksed 
     108         DO ji = 1, jpoce 
     109            IF (accmask(ji) == 0 ) THEN 
     110               zlimno3(ji,jk) = ( 1.0 - zlimo2(ji,jk) ) * pwcp(ji,jk,jwno3) / ( pwcp(ji,jk,jwno3) + xksedno3 )  
     111               zlimno3(ji,jk) = MAX(0., zlimno3(ji,jk) ) 
     112               zreasat = rearatpom(ji,jk) * zlimno3(ji,jk) * srDnit 
     113               ! For nitrates 
     114               pwcpa(ji,jk,jwno3) = pwcpa(ji,jk,jwno3) - zreasat 
     115               ! For alkalinity 
     116               pwcpa(ji,jk,jwalk) = pwcpa(ji,jk,jwalk) + zreasat 
     117            ENDIF 
     118         END DO 
     119      ENDDO 
     120 
     121 
     122      !-------------------------------------------------------------------- 
     123      ! Begining POC iron reduction 
     124      !-------------------------------------------------------------------- 
     125 
     126      DO jk = 2, jpksed 
     127         DO ji = 1, jpoce 
     128            IF (accmask(ji) == 0 ) THEN 
     129               zlimfeo(ji,jk) = ( 1.0 - zlimno3(ji,jk) - zlimo2(ji,jk) ) * solcp(ji,jk,jsfeo) / ( solcp(ji,jk,jsfeo) + xksedfeo ) 
     130               zlimfeo(ji,jk) = MAX(0., zlimfeo(ji,jk) ) 
     131               zreasat = 4.0 * rearatpom(ji,jk) * zlimfeo(ji,jk) 
     132               ! For FEOH 
     133               solcpa(ji,jk,jsfeo) = solcpa(ji,jk,jsfeo) - zreasat / volc(ji,jk,jsfeo) 
     134               ! For PO4 
     135               pwcpa(ji,jk,jwpo4)  = pwcpa(ji,jk,jwpo4) + zreasat * redfep 
     136               ! For alkalinity 
     137               pwcpa(ji,jk,jwalk)  = pwcpa(ji,jk,jwalk) + 2.0 * zreasat 
     138               ! Iron 
     139               pwcpa(ji,jk,jwfe2)  = pwcpa(ji,jk,jwfe2) + zreasat * radssol(jk,jwfe2) 
     140            ENDIF 
     141         END DO 
    200142      ENDDO 
    201143 
    202144      !-------------------------------------------------------------------- 
    203145      ! Begining POC denitrification and NO3- diffusion 
    204       ! (indice n°5 for couple POC/NO3- ie solcp(:,:,jspoc)/pwcp(:,:,jwno3)) 
    205       !-------------------------------------------------------------------- 
    206  
    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 
    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 
    277             ! 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 
    289             ! 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 
    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 
    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) 
    375             ! 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 
    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 
    386             ! 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 
    390             pwcp(ji,jk,jwfe2)  = pwcp(ji,jk,jwfe2) + zreasat * 4.0 
    391          ENDDO 
    392       ENDDO 
    393  
    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) 
    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 
    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 
    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  
    506       DO ji = 1, jpoce 
    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 ) 
     146      !-------------------------------------------------------------------- 
     147 
     148      DO jk = 2, jpksed 
     149         DO ji = 1, jpoce 
     150            IF (accmask(ji) == 0 ) THEN 
     151               zlimso4(ji,jk) = ( 1.0 - zlimno3(ji,jk) - zlimo2(ji,jk) - zlimfeo(ji,jk) ) * pwcp(ji,jk,jwso4) / ( pwcp(ji,jk,jwso4) + xksedso4 ) 
     152               zlimso4(ji,jk) = MAX(0., zlimso4(ji,jk) ) 
     153               zreasat = 0.5 * rearatpom(ji,jk) * zlimso4(ji,jk) 
     154               ! For sulfur 
     155               pwcpa(ji,jk,jwso4)  =  pwcpa(ji,jk,jwso4) - zreasat  
     156               pwcpa(ji,jk,jwh2s)  = pwcpa(ji,jk,jwh2s) + zreasat 
     157               ! For alkalinity 
     158               pwcpa(ji,jk,jwalk)  = pwcpa(ji,jk,jwalk) + 2.0 * zreasat 
     159            ENDIF 
     160         END DO 
     161      ENDDO 
     162 
     163      ! Secondary redox reactions 
     164      ! ------------------------- 
     165 
     166      call sed_dsr_redoxb( accmask ) 
    538167 
    539168      IF( ln_timing )  CALL timing_stop('sed_dsr') 
     
    541170   END SUBROUTINE sed_dsr 
    542171 
    543    SUBROUTINE sed_dsr_redoxb 
     172   SUBROUTINE sed_dsr_redoxb( accmask ) 
    544173      !!---------------------------------------------------------------------- 
    545174      !!                   ***  ROUTINE sed_dsr_redox  *** 
     
    547176      !!  ** Purpose :  computes secondary redox reactions 
    548177      !! 
    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       !!  
    561178      !!   History : 
    562179      !!        !  18-08 (O. Aumont)  Original code 
     
    564181      !! Arguments 
    565182      ! --- 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 
     183      INTEGER, DIMENSION(jpoce), INTENT(in) :: accmask 
     184      INTEGER   ::  ji, jni, jnj, jk   ! dummy looop indices 
     185 
     186      REAL(wp)  ::  zalpha, zexcess, zh2seq, zsedfer, zreasat 
    570187      !! 
    571188      !!---------------------------------------------------------------------- 
     
    573190      IF( ln_timing )  CALL timing_start('sed_dsr_redoxb') 
    574191 
    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             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 )  & 
    597                &            + zalpha * exp( reac_fe2 * zalpha * dtsed2 / 2. ) ) 
    598             ENDIF 
    599             zsedtra(1) = zalpha + 0.25 * zsedtra(4)  
    600             zsedtra(5) = zbeta  - zsedtra(4) 
    601             pwcp(ji,jk,jwalk) = zgamma + 2.0 * zsedtra(4) 
    602             pwcp(ji,jk,jwpo4) = zdelta + redfep * zsedtra(4) 
    603             ! H2S + O2 
    604             zalpha = zsedtra(1) - 2.0 * zsedtra(2) 
    605             zbeta  = pwcp(ji,jk,jwso4) + zsedtra(2) 
    606             zgamma = pwcp(ji,jk,jwalk) - 2.0 * zsedtra(2) 
    607             IF ( zalpha == 0. ) THEN 
    608                zsedtra(2) = zsedtra(2) / ( 1.0 + zsedtra(2) * reac_h2s * dtsed2 / 2.0 ) 
    609             ELSE 
    610                zsedtra(2) = ( zsedtra(2) * zalpha ) / ( 2.0 * zsedtra(2) * ( exp( reac_h2s * zalpha * dtsed2 / 2. ) - 1.0 )  & 
    611                &            + zalpha * exp( reac_h2s * zalpha * dtsed2 / 2. ) ) 
    612             ENDIF 
    613             zsedtra(1) = zalpha + 2.0 * zsedtra(2) 
    614             pwcp(ji,jk,jwalk) = zgamma + 2.0 * zsedtra(2) 
    615             pwcp(ji,jk,jwso4) = zbeta - zsedtra(2) 
    616             ! 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) 
    627             ! FeS - O2 
    628             zalpha = zsedtra(1) - 2.0 * zsedtra(6) 
    629             zbeta  = zsedtra(4) + zsedtra(6) 
    630             zgamma = pwcp(ji,jk,jwso4) + zsedtra(6) 
    631             IF ( zalpha == 0. ) THEN 
    632                zsedtra(6) = zsedtra(6) / ( 1.0 + zsedtra(6) * reac_feso * dtsed2 / 2.0 ) 
    633             ELSE 
    634                zsedtra(6) = ( zsedtra(6) * zalpha ) / ( 2.0 * zsedtra(6) * ( exp( reac_feso * zalpha * dtsed2 / 2. ) - 1.0 )  & 
    635                &            + zalpha * exp( reac_feso * zalpha * dtsed2 /2. ) ) 
    636             ENDIF 
    637             zsedtra(1) = zalpha + 2.0 * zsedtra(6) 
    638             zsedtra(4) = zbeta  - zsedtra(6) 
    639             pwcp(ji,jk,jwso4) = zgamma - zsedtra(6) 
    640 !            ! Fe - H2S 
    641             zalpha = zsedtra(2) - zsedtra(4) 
    642             zbeta  = zsedtra(4) + zsedtra(6) 
    643             zgamma = pwcp(ji,jk,jwalk) - 2.0 * zsedtra(4) 
    644             IF ( zalpha == 0. ) THEN 
    645                zsedtra(4) = zsedtra(4) / ( 1.0 + zsedtra(4) * reac_fes * dtsed2 / 2.0 ) 
    646             ELSE 
    647                zsedtra(4) = ( zsedtra(4) * zalpha ) / ( zsedtra(4) * ( exp( reac_fes * zalpha * dtsed2 / 2. ) - 1.0 )  & 
    648                &            + zalpha * exp( reac_fes * zalpha * dtsed2 /2. ) ) 
    649             ENDIF 
    650             zsedtra(2) = zalpha + zsedtra(4) 
    651             zsedtra(6) = zbeta  - zsedtra(4) 
    652             pwcp(ji,jk,jwalk) = zgamma + 2.0 * zsedtra(4) 
    653             ! FEOH + H2S 
    654             zalpha = zsedtra(5) - 2.0 * zsedtra(2) 
    655             zbeta  = zsedtra(5) + zsedtra(4) 
    656             zgamma = pwcp(ji,jk,jwalk) - 2.0 * zsedtra(4) 
    657             zdelta = pwcp(ji,jk,jwso4) + zsedtra(2) 
    658             zepsi  = pwcp(ji,jk,jwpo4) + redfep * zsedtra(5) 
    659             IF ( zalpha == 0. ) THEN 
    660                zsedtra(2) = zsedtra(2) / ( 1.0 + zsedtra(2) * reac_feh2s * dtsed2 ) 
    661             ELSE 
    662                zsedtra(2) = ( zsedtra(2) * zalpha ) / ( 2.0 * zsedtra(2) * ( exp( reac_feh2s * zalpha * dtsed2 ) - 1.0 )  & 
    663                &            + zalpha * exp( reac_feh2s * zalpha * dtsed2 ) ) 
    664             ENDIF 
    665             zsedtra(5) = zalpha + 2.0 * zsedtra(2) 
    666             zsedtra(4) = zbeta  - zsedtra(5) 
    667             pwcp(ji,jk,jwso4) = zdelta - zsedtra(2) 
    668             pwcp(ji,jk,jwalk) = zgamma + 2.0 * zsedtra(4) 
    669             pwcp(ji,jk,jwpo4) = zepsi - redfep * zsedtra(5) 
    670             ! Fe - H2S 
    671             zalpha = zsedtra(2) - zsedtra(4) 
    672             zbeta  = zsedtra(4) + zsedtra(6) 
    673             zgamma = pwcp(ji,jk,jwalk) - 2.0 * zsedtra(4) 
    674             IF ( zalpha == 0. ) THEN 
    675                zsedtra(4) = zsedtra(4) / ( 1.0 + zsedtra(4) * reac_fes * dtsed2 / 2.0 ) 
    676             ELSE 
    677                zsedtra(4) = ( zsedtra(4) * zalpha ) / ( zsedtra(4) * ( exp( reac_fes * zalpha * dtsed2 / 2. ) - 1.0 )  & 
    678                &            + zalpha * exp( reac_fes * zalpha * dtsed2 /2. ) ) 
    679             ENDIF 
    680             zsedtra(2) = zalpha + zsedtra(4) 
    681             zsedtra(6) = zbeta  - zsedtra(4) 
    682             pwcp(ji,jk,jwalk) = zgamma + 2.0 * zsedtra(4) 
    683             ! FeS - O2 
    684             zalpha = zsedtra(1) - 2.0 * zsedtra(6) 
    685             zbeta  = zsedtra(4) + zsedtra(6) 
    686             zgamma = pwcp(ji,jk,jwso4) + zsedtra(6) 
    687             IF (zalpha == 0.) THEN 
    688                zsedtra(6) = zsedtra(6) / ( 1.0 + zsedtra(6) * reac_feso * dtsed2 / 2. ) 
    689             ELSE 
    690                zsedtra(6) = ( zsedtra(6) * zalpha ) / ( 2.0 * zsedtra(6) * ( exp( reac_feso * zalpha * dtsed2 / 2. ) - 1.0 )  & 
    691                &            + zalpha * exp( reac_feso * zalpha * dtsed2 /2. ) ) 
    692             ENDIF 
    693             zsedtra(1) = zalpha + 2.0 * zsedtra(6) 
    694             zsedtra(4) = zbeta  - zsedtra(6) 
    695             pwcp(ji,jk,jwso4) = zgamma - zsedtra(6) 
    696             ! 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) 
    707             ! H2S + O2 
    708             zalpha = zsedtra(1) - 2.0 * zsedtra(2) 
    709             zbeta  = pwcp(ji,jk,jwso4) + zsedtra(2) 
    710             zgamma = pwcp(ji,jk,jwalk) - 2.0 * zsedtra(2) 
    711             IF ( zalpha == 0. ) THEN 
    712                zsedtra(2) = zsedtra(2) / ( 1.0 + zsedtra(2) * reac_h2s * dtsed2 / 2.0 ) 
    713             ELSE 
    714                zsedtra(2) = ( zsedtra(2) * zalpha ) / ( 2.0 * zsedtra(2) * ( exp( reac_h2s * zalpha * dtsed2 / 2. ) - 1.0 )  & 
    715                &            + zalpha * exp( reac_h2s * zalpha * dtsed2 / 2. ) ) 
    716             ENDIF 
    717             zsedtra(1) = zalpha + 2.0 * zsedtra(2) 
    718             pwcp(ji,jk,jwso4) = zbeta - zsedtra(2) 
    719             pwcp(ji,jk,jwalk) = zgamma + 2.0 * zsedtra(2) 
    720             ! Fe + O2 
    721             zalpha = zsedtra(1) - 0.25 * zsedtra(4) 
    722             zbeta  = zsedtra(4) + zsedtra(5) 
    723             zgamma = pwcp(ji,jk,jwalk) - 2.0 * zsedtra(4) 
    724             zdelta = pwcp(ji,jk,jwpo4) - redfep * zsedtra(4) 
    725             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 )  & 
    729                &            + zalpha * exp( reac_fe2 * zalpha * dtsed2 / 2. ) ) 
    730             ENDIF 
    731             zsedtra(1) = zalpha + 0.25 * zsedtra(4) 
    732             zsedtra(5) = zbeta  - zsedtra(4) 
    733             pwcp(ji,jk,jwpo4) = zdelta + redfep * zsedtra(4) 
    734             pwcp(ji,jk,jwalk) = zgamma + 2.0 * zsedtra(4) 
    735             pwcp(ji,jk,jwoxy)  = zsedtra(1) 
    736             pwcp(ji,jk,jwh2s)  = zsedtra(2) 
    737             pwcp(ji,jk,jwnh4)  = zsedtra(3) 
    738             pwcp(ji,jk,jwfe2)  = zsedtra(4) + sedligand(ji,jk) + zsedfer 
    739             pwcp(ji,jk,jwno3)  = pwcp(ji,jk,jwno3)  + ( zsedtrn(3) - pwcp(ji,jk,jwnh4) ) 
    740             solcp(ji,jk,jsfeo) = zsedtra(5) / zvolc(ji,jk,jsfeo) 
    741             solcp(ji,jk,jsfes) = zsedtra(6) / zvolc(ji,jk,jsfes) 
    742          END DO 
    743       END DO 
     192      DO jk = 2, jpksed 
     193         DO ji = 1, jpoce 
     194            IF (accmask(ji) == 0 ) THEN 
     195               zalpha = ( pwcp(ji,jk,jwfe2) - pwcp(ji,jk,jwlgw) ) * 1E9 
     196               zsedfer = ( zalpha + SQRT( zalpha**2 + 1.E-10 ) ) /2.0 * 1E-9 
     197 
     198               ! First pass of the scheme. At the end, it is 1st order  
     199               ! ----------------------------------------------------- 
     200               ! Fe (both adsorbed and solute) + O2 
     201               zreasat = reac_fe2 * zsedfer / radssol(jk,jwfe2) * pwcp(ji,jk,jwoxy) 
     202               pwcpa(ji,jk,jwfe2) = pwcpa(ji,jk,jwfe2) - zreasat * radssol(jk,jwfe2) 
     203               pwcpa(ji,jk,jwoxy) = pwcpa(ji,jk,jwoxy) - 0.25 * zreasat  
     204               pwcpa(ji,jk,jwpo4) = pwcpa(ji,jk,jwpo4) - redfep * zreasat  
     205               pwcpa(ji,jk,jwalk) = pwcpa(ji,jk,jwalk) - 2.0 * zreasat 
     206               solcpa(ji,jk,jsfeo) = solcpa(ji,jk,jsfeo) + zreasat / volc(ji,jk,jsfeo) 
     207 
     208               ! H2S + O2 
     209               zreasat = reac_h2s * pwcp(ji,jk,jwoxy) * pwcp(ji,jk,jwh2s) 
     210               pwcpa(ji,jk,jwh2s) = pwcpa(ji,jk,jwh2s) - zreasat 
     211               pwcpa(ji,jk,jwoxy) = pwcpa(ji,jk,jwoxy) - 2.0 * zreasat 
     212               pwcpa(ji,jk,jwso4) = pwcpa(ji,jk,jwso4) + zreasat 
     213               pwcpa(ji,jk,jwalk) = pwcpa(ji,jk,jwalk) - 2.0 * zreasat 
     214 
     215               ! NH4 + O2 
     216               zreasat = reac_nh4 * pwcp(ji,jk,jwnh4) / radssol(jk,jwnh4) * pwcp(ji,jk,jwoxy) 
     217               pwcpa(ji,jk,jwnh4) = pwcpa(ji,jk,jwnh4) - zreasat * radssol(jk,jwnh4)  
     218               pwcpa(ji,jk,jwoxy) = pwcpa(ji,jk,jwoxy) - 2.0 * zreasat  
     219               pwcpa(ji,jk,jwno3) = pwcpa(ji,jk,jwno3) + zreasat 
     220               pwcpa(ji,jk,jwalk) = pwcpa(ji,jk,jwalk) - 2.0 * zreasat 
     221 
     222               ! FeS - O2 
     223               zreasat = reac_feso * solcp(ji,jk,jsfes) * volc(ji,jk,jsfes) * pwcp(ji,jk,jwoxy)  
     224               solcpa(ji,jk,jsfes) = solcpa(ji,jk,jsfes) - zreasat / volc(ji,jk,jsfes) 
     225               pwcpa(ji,jk,jwoxy) = pwcpa(ji,jk,jwoxy) - 2.0 * zreasat 
     226               pwcpa(ji,jk,jwfe2) = pwcpa(ji,jk,jwfe2) + zreasat * radssol(jk,jwfe2) 
     227               pwcpa(ji,jk,jwso4) = pwcpa(ji,jk,jwso4) + zreasat 
     228 
     229               ! FEOH + H2S 
     230               zreasat = reac_feh2s * solcp(ji,jk,jsfeo) * volc(ji,jk,jsfeo) * pwcp(ji,jk,jwh2s)  
     231               solcpa(ji,jk,jsfeo) = solcpa(ji,jk,jsfeo) - 8.0 * zreasat / volc(ji,jk,jsfeo) 
     232               pwcpa(ji,jk,jwh2s) = pwcpa(ji,jk,jwh2s) - zreasat 
     233               pwcpa(ji,jk,jwfe2) = pwcpa(ji,jk,jwfe2) + 8.0 * zreasat * radssol(jk,jwfe2)  
     234               pwcpa(ji,jk,jwalk) = pwcpa(ji,jk,jwalk) + 14.0 * zreasat  
     235               pwcpa(ji,jk,jwso4) = pwcpa(ji,jk,jwso4) + zreasat 
     236               pwcpa(ji,jk,jwpo4) = pwcpa(ji,jk,jwpo4) + 8.0 * redfep * zreasat 
     237 
     238               ! Fe + H2S 
     239               zh2seq     = MAX(rtrn, 2.1E-3 * hipor(ji,jk)) 
     240               zexcess = pwcp(ji,jk,jwh2s) * zsedfer / zh2seq - 1.0 
     241               IF ( zexcess >= 0.0 ) THEN 
     242                  zreasat = reac_fesp * zexcess 
     243                  pwcpa (ji,jk,jwfe2) = pwcpa(ji,jk,jwfe2) - zreasat * radssol(jk,jwfe2) 
     244                  solcpa(ji,jk,jsfes) = solcpa(ji,jk,jsfes) + zreasat / volc(ji,jk,jsfes) 
     245                  pwcpa(ji,jk,jwh2s)  = pwcpa(ji,jk,jwh2s) - zreasat 
     246               ELSE 
     247                  zreasat = reac_fesd * zexcess * solcp(ji,jk,jsfes) * volc(ji,jk,jsfes) 
     248                  pwcpa (ji,jk,jwfe2) = pwcpa(ji,jk,jwfe2) - zreasat * radssol(jk,jwfe2) 
     249                  solcpa(ji,jk,jsfes) = solcpa(ji,jk,jsfes) + zreasat / volc(ji,jk,jsfes) 
     250                  pwcpa(ji,jk,jwh2s)  = pwcpa(ji,jk,jwh2s) - zreasat 
     251               ENDIF 
     252               ! For alkalinity 
     253               pwcpa(ji,jk,jwalk)  = pwcpa(ji,jk,jwalk) - zreasat * 2.0 
     254            ENDIF 
     255         END DO 
     256     END DO 
    744257 
    745258      IF( ln_timing )  CALL timing_stop('sed_dsr_redoxb') 
  • NEMO/trunk/src/TOP/PISCES/SED/seddta.F90

    r15237 r15450  
    88   USE sed 
    99   USE sedarr 
    10    USE par_pisces 
     10   USE sedini 
    1111   USE phycst, ONLY : rday 
    1212   USE iom 
     
    2020 
    2121   !! *  Module variables 
    22    REAL(wp) ::  rsecday  ! number of second per a day 
    2322   REAL(wp) ::  conv2    ! [kg/m2/month]-->[g/cm2/s] ( 1 month has 30 days ) 
    2423 
     
    2625#  include "do_loop_substitute.h90" 
    2726#  include "domzgr_substitute.h90" 
    28    !! $Id$ 
     27 
    2928CONTAINS 
    3029 
     
    4746 
    4847      !! Arguments 
    49       INTEGER, INTENT(in) ::  kt         ! time-step 
    50       INTEGER, INTENT(in) ::  Kbb, Kmm  ! time level indices 
     48      INTEGER, INTENT( in ) ::   kt    ! time-step 
     49      INTEGER, INTENT( in ) ::   Kbb, Kmm ! time level indices 
    5150 
    5251      !! * Local declarations 
     
    5453 
    5554      REAL(wp), DIMENSION(jpoce) :: zdtap, zdtag 
    56       REAL(wp), DIMENSION(jpi,jpj) :: zwsbio4, zwsbio3 
     55      REAL(wp), DIMENSION(jpi,jpj) :: zwsbio4, zwsbio3, zddust 
    5756      REAL(wp) :: zf0, zf1, zf2, zkapp, zratio, zdep 
     57      REAL(wp) :: zzf0, zf0s, zf0b, zzf1, zf1s, zf1b, zzf2, zf2s, zf2b 
    5858 
    5959      !---------------------------------------------------------------------- 
     
    7878         IF (lwp) WRITE(numsed,*) ' sed_dta : Sediment fields' 
    7979         dtsed = rDt_trc 
    80          rsecday = 60.* 60. * 24. 
    81 !         conv2   = 1.0e+3 / ( 1.0e+4 * rsecday * 30. ) 
    8280         conv2 = 1.0e+3 /  1.0e+4  
    83          rdtsed(2:jpksed) = dtsed / ( denssol * por1(2:jpksed) ) 
    8481      ENDIF 
    8582 
     
    8784      zdtap(:)    = 0.  
    8885      zdtag(:)    = 0.   
     86      zddust(:,:) = 0.0 
    8987 
    9088      ! reading variables 
     
    9795      !    ----------------------------------------------------------- 
    9896      IF (ln_sediment_offline) THEN 
    99          DO_2D( 1, 1, 1, 1 ) 
     97         DO_2D( 0, 0, 0, 0 ) 
    10098            ikt = mbkt(ji,jj) 
    10199            zwsbio4(ji,jj) = wsbio2 / rday 
     
    103101         END_2D 
    104102      ELSE 
    105          DO_2D( 1, 1, 1, 1 ) 
     103         DO_2D( 0, 0, 0, 0 ) 
    106104            ikt = mbkt(ji,jj) 
    107105            zdep = e3t(ji,jj,ikt,Kmm) / rDt_trc 
     
    112110 
    113111      trc_data(:,:,:) = 0. 
    114       DO_2D( 1, 1, 1, 1 ) 
     112      DO_2D( 0, 0, 0, 0 ) 
    115113         ikt = mbkt(ji,jj) 
    116          IF ( tmask(ji,jj,ikt) == 1 ) THEN 
    117             trc_data(ji,jj,1)   = tr(ji,jj,ikt,jpsil,Kbb) 
    118             trc_data(ji,jj,2)   = tr(ji,jj,ikt,jpoxy,Kbb) 
    119             trc_data(ji,jj,3)   = tr(ji,jj,ikt,jpdic,Kbb) 
    120             trc_data(ji,jj,4)   = tr(ji,jj,ikt,jpno3,Kbb) / 7.625 
    121             trc_data(ji,jj,5)   = tr(ji,jj,ikt,jppo4,Kbb) / 122. 
    122             trc_data(ji,jj,6)   = tr(ji,jj,ikt,jptal,Kbb) 
    123             trc_data(ji,jj,7)   = tr(ji,jj,ikt,jpnh4,Kbb) / 7.625 
    124             trc_data(ji,jj,8)   = 0.0 
    125             trc_data(ji,jj,9)   = 28.0E-3 
    126             trc_data(ji,jj,10)  = tr(ji,jj,ikt,jpfer,Kbb) 
    127             trc_data(ji,jj,11 ) = MIN(tr(ji,jj,ikt,jpgsi,Kbb), 1E-4) * zwsbio4(ji,jj) * 1E3 
    128             trc_data(ji,jj,12 ) = MIN(tr(ji,jj,ikt,jppoc,Kbb), 1E-4) * zwsbio3(ji,jj) * 1E3 
    129             trc_data(ji,jj,13 ) = MIN(tr(ji,jj,ikt,jpgoc,Kbb), 1E-4) * zwsbio4(ji,jj) * 1E3 
    130             trc_data(ji,jj,14)  = MIN(tr(ji,jj,ikt,jpcal,Kbb), 1E-4) * zwsbio4(ji,jj) * 1E3 
    131             trc_data(ji,jj,15)  = ts(ji,jj,ikt,jp_tem,Kmm) 
    132             trc_data(ji,jj,16)  = ts(ji,jj,ikt,jp_sal,Kmm) 
    133             trc_data(ji,jj,17 ) = ( tr(ji,jj,ikt,jpsfe,Kbb) * zwsbio3(ji,jj) + tr(ji,jj,ikt,jpbfe,Kbb)  & 
    134             &                     * zwsbio4(ji,jj)  ) * 1E3 / ( trc_data(ji,jj,12 ) + trc_data(ji,jj,13 ) + rtrn ) 
    135             trc_data(ji,jj,17 ) = MIN(1E-3, trc_data(ji,jj,17 ) ) 
     114         IF ( tmask(ji,jj,ikt) == 1.0 ) THEN 
     115            trc_data(ji,jj,jwsil) = tr(ji,jj,ikt,jpsil,Kbb) 
     116            trc_data(ji,jj,jwoxy) = tr(ji,jj,ikt,jpoxy,Kbb) 
     117            trc_data(ji,jj,jwdic) = tr(ji,jj,ikt,jpdic,Kbb) 
     118            trc_data(ji,jj,jwno3) = tr(ji,jj,ikt,jpno3,Kbb) * redNo3 / redC 
     119            trc_data(ji,jj,jwpo4) = tr(ji,jj,ikt,jppo4,Kbb) / redC 
     120            trc_data(ji,jj,jwalk) = tr(ji,jj,ikt,jptal,Kbb)  
     121            trc_data(ji,jj,jwnh4) = tr(ji,jj,ikt,jpnh4,Kbb) * redNo3 / redC  
     122            trc_data(ji,jj,jwh2s) = 0.0 
     123            trc_data(ji,jj,jwso4) = 0.14 * ts(ji,jj,ikt,jp_sal,Kmm) / 1.80655 / 96.062 
     124            trc_data(ji,jj,jwfe2) = tr(ji,jj,ikt,jpfer,Kbb) 
     125            trc_data(ji,jj,jwlgw) = 1E-9 
     126            trc_data(ji,jj,12 )   = MIN(tr(ji,jj,ikt,jpgsi,Kbb), 1E-4) * zwsbio4(ji,jj) * 1E3 
     127            trc_data(ji,jj,13 )   = MIN(tr(ji,jj,ikt,jppoc,Kbb), 1E-4) * zwsbio3(ji,jj) * 1E3 
     128            trc_data(ji,jj,14 )   = MIN(tr(ji,jj,ikt,jpgoc,Kbb), 1E-4) * zwsbio4(ji,jj) * 1E3 
     129            trc_data(ji,jj,15)    = MIN(tr(ji,jj,ikt,jpcal,Kbb), 1E-4) * zwsbio4(ji,jj) * 1E3 
     130            trc_data(ji,jj,16)    = ts(ji,jj,ikt,jp_tem,Kmm) 
     131            trc_data(ji,jj,17)    = ts(ji,jj,ikt,jp_sal,Kmm) 
     132            trc_data(ji,jj,18 )   = ( tr(ji,jj,ikt,jpsfe,Kbb) * zwsbio3(ji,jj) + tr(ji,jj,ikt,jpbfe,Kbb)  & 
     133            &                       * zwsbio4(ji,jj)  ) * 1E3 / ( trc_data(ji,jj,13 ) + trc_data(ji,jj,14 ) + rtrn ) 
     134            trc_data(ji,jj,18 )   = MIN(1E-3, trc_data(ji,jj,18 ) ) 
    136135         ENDIF 
    137136      END_2D 
     
    142141         CALL pack_arr ( jpoce,  pwcp_dta(1:jpoce,jw), trc_data(1:jpi,1:jpj,jw), iarroce(1:jpoce) ) 
    143142      END DO 
     143 
    144144      !  Solid components :  
    145145      !----------------------- 
    146146      !  Sinking fluxes for OPAL in mol.m-2.s-1 ; conversion in mol.cm-2.s-1 
    147       CALL pack_arr ( jpoce, rainrm_dta(1:jpoce,jsopal), trc_data(1:jpi,1:jpj,11), iarroce(1:jpoce) )  
     147      CALL pack_arr ( jpoce, rainrm_dta(1:jpoce,jsopal), trc_data(1:jpi,1:jpj,12), iarroce(1:jpoce) )  
    148148      rainrm_dta(1:jpoce,jsopal) = rainrm_dta(1:jpoce,jsopal) * 1e-4 
     149 
    149150      !  Sinking fluxes for POC in mol.m-2.s-1 ; conversion in mol.cm-2.s-1 
    150       CALL pack_arr ( jpoce, zdtap(1:jpoce), trc_data(1:jpi,1:jpj,12) , iarroce(1:jpoce) )       
    151       CALL pack_arr ( jpoce, zdtag(1:jpoce), trc_data(1:jpi,1:jpj,13) , iarroce(1:jpoce) ) 
     151      CALL pack_arr ( jpoce, zdtap(1:jpoce), trc_data(1:jpi,1:jpj,13) , iarroce(1:jpoce) )       
     152      CALL pack_arr ( jpoce, zdtag(1:jpoce), trc_data(1:jpi,1:jpj,14) , iarroce(1:jpoce) ) 
    152153      DO ji = 1, jpoce 
    153 !        zkapp  = MIN( (1.0 - 0.02 ) * reac_poc, 3731.0 * max(100.0, zkbot(ji) )**(-1.011) / ( 365.0 * 24.0 * 3600.0 ) ) 
    154 !        zkapp   = MIN( 0.98 * reac_poc, 100.0 * max(100.0, zkbot(ji) )**(-0.6) / ( 365.0 * 24.0 * 3600.0 ) ) 
    155 !        zratio = ( ( 1.0 - 0.02 ) * reac_poc + 0.02 * reac_poc * 0. - zkapp) / ( ( 0.02 - 1.0 ) * reac_poc / 100. - 0.02 * reac_poc * 0. + zkapp ) 
    156 !        zf1    = ( 0.02 * (reac_poc - reac_poc * 0.) + zkapp - reac_poc ) / ( reac_poc / 100. - reac_poc ) 
    157 !        zf1    = MIN(0.98, MAX(0., zf1 ) ) 
    158          zf1    = 0.48 
    159          zf0    = 1.0 - 0.02 - zf1 
    160          zf2    = 0.02 
    161          rainrm_dta(ji,jspoc) =   ( zdtap(ji) +  zdtag(ji) ) * 1e-4 * zf0 
    162          rainrm_dta(ji,jspos) =   ( zdtap(ji) +  zdtag(ji) ) * 1e-4 * zf1 
    163          rainrm_dta(ji,jspor) =   ( zdtap(ji) +  zdtag(ji) ) * 1e-4 * zf2 
     154         zzf2 = 2E-2 
     155         zzf1 = 0.3 
     156         zzf0 = 1.0 - zzf1 - zzf2 
     157         zf0s = zzf0 
     158         zf1s = zzf1 
     159         zf2s = 1.0 - zf1s - zf0s 
     160         zf0b = zzf0 
     161         zf1b = zzf1 
     162         zf2b = 1.0 - zf1b - zf0b 
     163         rainrm_dta(ji,jspoc) =   ( zdtap(ji) * zf0s +  zdtag(ji) * zf0b ) * 1e-4 
     164         rainrm_dta(ji,jspos) =   ( zdtap(ji) * zf1s +  zdtag(ji) * zf1b ) * 1e-4 
     165         rainrm_dta(ji,jspor) =   ( zdtap(ji) * zf2s +  zdtag(ji) * zf2b ) * 1e-4 
    164166      END DO 
     167 
    165168      !  Sinking fluxes for Calcite in mol.m-2.s-1 ; conversion in mol.cm-2.s-1 
    166       CALL pack_arr ( jpoce,  rainrm_dta(1:jpoce,jscal), trc_data(1:jpi,1:jpj,14), iarroce(1:jpoce) ) 
     169      CALL pack_arr ( jpoce,  rainrm_dta(1:jpoce,jscal), trc_data(1:jpi,1:jpj,15), iarroce(1:jpoce) ) 
    167170      rainrm_dta(1:jpoce,jscal) = rainrm_dta(1:jpoce,jscal) * 1e-4 
    168       ! vector temperature [�C] and salinity  
    169       CALL pack_arr ( jpoce,  temp(1:jpoce), trc_data(1:jpi,1:jpj,15), iarroce(1:jpoce) ) 
    170       CALL pack_arr ( jpoce,  salt(1:jpoce), trc_data(1:jpi,1:jpj,16), iarroce(1:jpoce) ) 
     171 
     172      ! vector temperature [°C] and salinity  
     173      CALL pack_arr ( jpoce,  temp(1:jpoce), trc_data(1:jpi,1:jpj,16), iarroce(1:jpoce) ) 
     174      CALL pack_arr ( jpoce,  salt(1:jpoce), trc_data(1:jpi,1:jpj,17), iarroce(1:jpoce) ) 
    171175       
    172176      ! Clay rain rate in [mol/(cm**2.s)]  
     
    175179      CALL pack_arr ( jpoce,  rainrm_dta(1:jpoce,jsclay), dust(1:jpi,1:jpj), iarroce(1:jpoce) ) 
    176180      rainrm_dta(1:jpoce,jsclay) = rainrm_dta(1:jpoce,jsclay) * conv2 / mol_wgt(jsclay)   & 
    177       &                            + wacc(1:jpoce) * por1(2) * denssol / mol_wgt(jsclay) / ( rsecday * 365.0 ) 
    178       rainrm_dta(1:jpoce,jsclay) = rainrm_dta(1:jpoce,jsclay) * 0.965 
    179       rainrm_dta(1:jpoce,jsfeo)  = rainrm_dta(1:jpoce,jsclay) * mol_wgt(jsclay) / mol_wgt(jsfeo) * 0.035 / 0.965 
     181      &                            + wacc(1:jpoce) * por1(2) * dens_sol(jsclay) / mol_wgt(jsclay) / ( rday * 365.0 ) 
     182      rainrm_dta(1:jpoce,jsfeo)  = rainrm_dta(1:jpoce,jsclay) * mol_wgt(jsclay) / mol_wgt(jsfeo) * 0.035 * 0.5 
     183      rainrm_dta(1:jpoce,jsclay) = rainrm_dta(1:jpoce,jsclay) * ( 1.0 - 0.035 * 0.5 )  
     184      CALL unpack_arr ( jpoce, zddust(1:jpi,1:jpj), iarroce(1:jpoce), wacc(1:jpoce) ) 
     185      zddust(:,:) = dust(:,:) + zddust(:,:) / ( rday * 365.0 ) * por1(2) * dens_sol(jsclay) / conv2 
     186 
    180187!    rainrm_dta(1:jpoce,jsclay) = 1.0E-4 * conv2 / mol_wgt(jsclay) 
    181188 
     
    184191 
    185192      ! Fe/C ratio in sinking particles that fall to the sediments 
    186       CALL pack_arr ( jpoce,  fecratio(1:jpoce), trc_data(1:jpi,1:jpj,17), iarroce(1:jpoce) ) 
    187  
    188       sedligand(:,1) = 1.E-9 
     193      CALL pack_arr ( jpoce,  fecratio(1:jpoce), trc_data(1:jpi,1:jpj,18), iarroce(1:jpoce) ) 
    189194 
    190195      ! sediment pore water at 1st layer (k=1) 
    191       DO jw = 1, jpwat 
    192          pwcp(1:jpoce,1,jw) = pwcp_dta(1:jpoce,jw) 
    193       ENDDO 
    194  
    195       !  rain 
    196       DO js = 1, jpsol 
    197          rainrm(1:jpoce,js) = rainrm_dta(1:jpoce,js) 
    198       ENDDO 
     196      pwcp(1:jpoce,1,1:jpwat) = pwcp_dta(1:jpoce,1:jpwat) 
    199197 
    200198      ! Calculation of raintg of each sol. comp.: rainrm in [g/(cm**2.s)] 
    201199      DO js = 1, jpsol 
    202          rainrg(1:jpoce,js) = rainrm(1:jpoce,js) * mol_wgt(js) 
     200         rainrg(1:jpoce,js) = rainrm_dta(1:jpoce,js) * mol_wgt(js) 
    203201      ENDDO 
    204202 
    205       ! Calculation of raintg = total massic flux rained in each cell (sum of sol. comp.) 
    206       raintg(:) = 0. 
     203      ! computation of dzdep = total thickness of solid material rained [cm] in each cell 
     204      dzdep(:) = 0. 
    207205      DO js = 1, jpsol 
    208          raintg(1:jpoce) = raintg(1:jpoce) + rainrg(1:jpoce,js) 
    209       ENDDO 
    210  
    211       ! computation of dzdep = total thickness of solid material rained [cm] in each cell 
    212       dzdep(1:jpoce) = raintg(1:jpoce) * rdtsed(2)  
     206         dzdep(1:jpoce) = dzdep(1:jpoce) + rainrg(1:jpoce,js) * dtsed / ( dens_sol(js) * por1(2) ) 
     207      END DO 
    213208 
    214209      IF( lk_iomput ) THEN 
    215           IF( iom_use("sflxclay" ) ) CALL iom_put( "sflxclay", dust(:,:) * conv2 * 1E4 ) 
    216           IF( iom_use("sflxcal" ) )  CALL iom_put( "sflxcal", trc_data(:,:,13) ) 
    217           IF( iom_use("sflxbsi" ) )  CALL iom_put( "sflxbsi", trc_data(:,:,10) ) 
    218           IF( iom_use("sflxpoc" ) )  CALL iom_put( "sflxpoc", trc_data(:,:,11) + trc_data(:,:,12) ) 
     210          IF( iom_use("sflxclay" ) ) CALL iom_put( "sflxclay", zddust(:,:) * 1E3 / 1.E4 ) 
     211          IF( iom_use("sflxcal" ) )  CALL iom_put( "sflxcal", trc_data(:,:,15) / 1.E4 ) 
     212          IF( iom_use("sflxbsi" ) )  CALL iom_put( "sflxbsi", trc_data(:,:,12) / 1.E4 ) 
     213          IF( iom_use("sflxpoc" ) )  CALL iom_put( "sflxpoc", ( trc_data(:,:,13) + trc_data(:,:,14) ) / 1.E4 ) 
    219214      ENDIF 
    220215 
  • NEMO/trunk/src/TOP/PISCES/SED/sedini.F90

    r14086 r15450  
    66 
    77   !!---------------------------------------------------------------------- 
    8    !!   sed_init    : initialization, namelist read, and parameters control 
     8   !!   sed_ini    : initialization, namelist read, and parameters control 
    99   !!---------------------------------------------------------------------- 
    1010   !! * Modules used 
    11    USE par_trc        ! need jptra, number of passive tracers 
    12    USE par_sed 
    1311   USE sed     ! sediment global variable 
    1412   USE sed_oce 
     
    2422   PRIVATE 
    2523 
    26    !! * Substitutions 
    27 #  include "do_loop_substitute.h90" 
    2824   !! Module variables 
     25   REAL(wp), PUBLIC :: sedmask  
     26 
    2927   REAL(wp)    ::  & 
    3028      sedzmin = 0.3    ,  &  !: Minimum vertical spacing 
     
    3937 
    4038   REAL(wp)    ::  & 
    41       rcopal  =   40.  ,  &  !: reactivity for si    [l.mol-1.an-1] 
    42       dcoef   =  8.e-6       !: diffusion coefficient (*por)   [cm**2/s] 
     39      rcopal  =   40.        !: reactivity for si    [l.mol-1.an-1] 
    4340 
    4441   REAL(wp), PUBLIC    ::  & 
    45       redO2    =  172.  ,  &  !: Redfield coef for Oxygen 
     42      redO2    =  138.  ,  &  !: Redfield coef for Oxygen 
    4643      redNo3   =   16.  ,  &  !: Redfield coef for Nitrate 
    4744      redPo4   =    1.  ,  &  !: Redfield coef for Phosphate 
    48       redC     =  122.  ,  &  !: Redfield coef for Carbon 
     45      redC     =  117.  ,  &  !: Redfield coef for Carbon 
    4946      redfep   =  0.175 ,  &  !: Ratio for iron bound phosphorus 
    5047      rcorgl   =   50.  ,  &  !: reactivity for POC/O2 [l.mol-1.an-1]     
     
    5552      rcfe2    =   5.E8 ,  &  !: reactivity for O2/Fe2+ [l.mol-1.an-1] 
    5653      rcfeh2s  =   1.E4 ,  &  !: Reactivity for FEOH/H2S [l.mol-1.an-1] 
    57       rcfes    =   1.E5 ,  &  !: Reactivity for FE2+/H2S [l.mol-1.an-1] 
    5854      rcfeso   =   3.E5 ,  &  !: Reactivity for FES/O2 [l.mol-1.an-1] 
     55      rcfesp   =   5E-6 ,  &  !: Precipitation of FeS [mol/l-1.an-1] 
     56      rcfesd   =   1E-3 ,  &  !: Dissolution of FeS [an-1] 
    5957      xksedo2  =   5E-6 ,  &  !: half-sturation constant for oxic remin. 
    6058      xksedno3 =   5E-6 ,  &  !: half-saturation constant for denitrification 
     
    7371 
    7472   REAL(wp), DIMENSION(jpwat), PUBLIC  :: diff1 
    75    DATA diff1/4.59E-6, 1.104E-5, 4.81E-6 , 9.78E-6, 3.58E-6, 4.01E-6, 9.8E-6, 9.73E-6, 5.0E-6, 3.31E-6 / 
     73   DATA diff1/ 1.104E-5, 9.78E-6, 3.58E-6, 9.8E-6, 9.73E-6, 5.0E-6, 3.31E-6, 4.81E-6, 4.81E-6, 4.81E-6, 4.59E-6 / 
     74 
    7675 
    7776   REAL(wp), DIMENSION(jpwat), PUBLIC  :: diff2 
    78    DATA diff2/1.74E-7, 4.47E-7, 2.51E-7, 3.89E-7, 1.77E-7, 2.5E-7, 3.89E-7, 3.06E-7, 2.5E-7, 1.5E-7 / 
    79  
    80  
     77   DATA diff2/ 4.47E-7, 3.89E-7, 1.77E-7, 3.89E-7, 3.06E-7, 2.5E-7, 1.5E-7, 2.51E-7, 2.51E-7, 2.51E-7, 1.74E-7 / 
    8178 
    8279   !! *  Routine accessibility 
    83    PUBLIC sed_init          ! routine called by opa.F90 
    84  
     80   PUBLIC sed_ini          ! routine called by opa.F90 
     81 
     82   !! * Substitutions 
     83#  include "do_loop_substitute.h90" 
    8584   !! $Id$ 
    8685CONTAINS 
    8786 
    8887 
    89    SUBROUTINE sed_init 
     88   SUBROUTINE sed_ini 
    9089      !!---------------------------------------------------------------------- 
    91       !!                   ***  ROUTINE sed_init  *** 
     90      !!                   ***  ROUTINE sed_ini  *** 
    9291      !! 
    9392      !! ** Purpose :  Initialization of sediment module 
     
    104103      !!        !  06-07  (C. Ethe)  Re-organization 
    105104      !!---------------------------------------------------------------------- 
    106       INTEGER :: ji, jj, ikt, ierr 
     105      INTEGER  :: ji, jj, js, jn, jk, ikt, ierr 
     106      REAL(wp) :: ztmp1, ztmp2  
    107107      !!---------------------------------------------------------------------- 
    108  
    109108 
    110109      ! Reading namelist.sed variables 
     
    122121      ENDIF 
    123122 
    124       IF(lwp) WRITE(numsed,*) ' sed_init : Initialization of sediment module  ' 
     123      IF(lwp) WRITE(numsed,*) ' sed_ini : Initialization of sediment module  ' 
    125124      IF(lwp) WRITE(numsed,*) ' ' 
    126125 
    127126      ! Read sediment Namelist 
    128127      !------------------------- 
    129       CALL sed_init_nam 
     128      CALL sed_ini_nam 
    130129 
    131130      ! Allocate SEDIMENT arrays 
    132131      ierr =        sed_alloc() 
    133132      ierr = ierr + sed_oce_alloc() 
    134       ierr = ierr + sed_adv_alloc()  
    135133      IF( ierr /= 0 )   CALL ctl_stop( 'STOP', 'sed_ini: unable to allocate sediment model arrays' ) 
    136134 
    137135      ! Determination of sediments number of points and allocate global variables 
    138136      epkbot(:,:) = 0. 
    139       DO_2D( 1, 1, 1, 1 ) 
     137      gdepbot(:,:) = 0. 
     138      DO_2D( 0, 0, 0, 0 ) 
    140139         ikt = mbkt(ji,jj)  
    141          IF( tmask(ji,jj,ikt) == 1 ) epkbot(ji,jj) = e3t_1d(ikt) 
    142          gdepbot(ji,jj) = gdepw_0(ji,jj,ikt) 
     140         IF( tmask(ji,jj,ikt) == 1 ) epkbot(ji,jj) = e3t_0(ji,jj,ikt) 
     141         gdepbot(ji,jj) = gdepw_0(ji,jj,ikt+1) 
    143142      END_2D 
    144143 
     
    154153 
    155154      ! Allocate memory size of global variables 
    156       ALLOCATE( pwcp (jpoce,jpksed,jpwat) )  ;  ALLOCATE( pwcp0 (jpoce,jpksed,jpwat) ) ;  ALLOCATE( pwcp_dta  (jpoce,jpwat) ) 
    157       ALLOCATE( solcp(jpoce,jpksed,jpsol) )  ;  ALLOCATE( solcp0(jpoce,jpksed,jpsol) ) ;  ALLOCATE( rainrm_dta(jpoce,jpsol) ) 
    158       ALLOCATE( rainrm(jpoce,jpsol) )        ;  ALLOCATE( rainrg(jpoce,jpsol) )        ;  ALLOCATE( raintg(jpoce) )  
     155      ALLOCATE( pwcp (jpoce,jpksed,jpwat) )  ;  ALLOCATE( pwcp_dta  (jpoce,jpwat) ) 
     156      ALLOCATE( pwcpa(jpoce,jpksed,jpwat) )  ;  ALLOCATE( solcpa(jpoce,jpksed,jpsol) ) 
     157      ALLOCATE( solcp(jpoce,jpksed,jpsol) )  ;  ALLOCATE( rainrm_dta(jpoce,jpsol) ) 
     158      ALLOCATE( rainrg(jpoce,jpsol) ) 
    159159      ALLOCATE( dzdep(jpoce) )               ;  ALLOCATE( iarroce(jpoce) )             ;  ALLOCATE( dzkbot(jpoce) ) 
     160      ALLOCATE( slatit(jpoce) )              ;  ALLOCATE( slongit(jpoce) ) 
    160161      ALLOCATE( zkbot(jpoce) )               ;  ALLOCATE( db(jpoce,jpksed) ) 
    161162      ALLOCATE( temp(jpoce) )                ;  ALLOCATE( salt(jpoce) )   
    162163      ALLOCATE( diff(jpoce,jpksed,jpwat ) )  ;  ALLOCATE( irrig(jpoce, jpksed) ) 
    163164      ALLOCATE( wacc(jpoce) )                ;  ALLOCATE( fecratio(jpoce) ) 
    164       ALLOCATE( press(jpoce) )               ;  ALLOCATE( densSW(jpoce) )  
     165      ALLOCATE( densSW(jpoce) )   ;   ALLOCATE( saturco3(jpoce,jpksed) )  
    165166      ALLOCATE( hipor(jpoce,jpksed) )        ;  ALLOCATE( co3por(jpoce,jpksed) ) 
    166167      ALLOCATE( dz3d(jpoce,jpksed) )         ;  ALLOCATE( volw3d(jpoce,jpksed) )       ;  ALLOCATE( vols3d(jpoce,jpksed) ) 
    167       ALLOCATE( sedligand(jpoce, jpksed) ) 
     168      ALLOCATE( rearatpom(jpoce, jpksed) )   ;  ALLOCATE( volc(jpoce,jpksed,jpsol) ) 
     169      ALLOCATE( radssol(jpksed, jpwat) )     ;  ALLOCATE( rads1sol(jpksed, jpwat) ) 
     170      ALLOCATE( apluss(jpoce, jpksed) )      ;  ALLOCATE( aminuss(jpoce,jpksed) ) 
    168171 
    169172      ! Initialization of global variables 
    170       pwcp  (:,:,:) = 0.   ;  pwcp0 (:,:,:) = 0.  ; pwcp_dta  (:,:) = 0.   
    171       solcp (:,:,:) = 0.   ;  solcp0(:,:,:) = 0.  ; rainrm_dta(:,:) = 0. 
    172       rainrm(:,:  ) = 0.   ;  rainrg(:,:  ) = 0.  ; raintg    (:  ) = 0.  
     173      pwcp  (:,:,:) = 0.   ;  pwcp_dta  (:,:) = 0.   
     174      pwcpa (:,:,:) = 0.   ;  solcpa(:,:,:) = 0. 
     175      solcp (:,:,:) = 0.   ;  rainrm_dta(:,:) = 0. 
     176      rainrg(:,:  ) = 0. 
    173177      dzdep (:    ) = 0.   ;  iarroce(:   ) = 0   ; dzkbot    (:  ) = 0. 
    174178      temp  (:    ) = 0.   ;  salt   (:   ) = 0.  ; zkbot     (:  ) = 0. 
    175       press (:    ) = 0.   ;  densSW (:   ) = 0.  ; db        (:,:) = 0.  
     179      densSW (:   ) = 0.   ;  db     (:,:) = 0.  
    176180      hipor (:,:  ) = 0.   ;  co3por (:,: ) = 0.  ; irrig     (:,:) = 0.  
    177181      dz3d  (:,:  ) = 0.   ;  volw3d (:,: ) = 0.  ; vols3d    (:,:) = 0.  
    178       fecratio(:)   = 1E-5  
    179       sedligand(:,:) = 0.6E-9 
     182      fecratio(:)   = 1E-5 ;  rearatpom(:,:) = 0.  
     183      radssol(:,:)  = 1.0  ;  rads1sol(:,:) = 0. 
     184      apluss(:,:)   = 0.0  ;  aminuss(:,:)  = 0.0 
    180185 
    181186      ! Chemical variables       
     
    185190      ALLOCATE( borats(jpoce) )  ;  ALLOCATE( calcon2(jpoce) )  ;  ALLOCATE( sieqs (jpoce) )  
    186191      ALLOCATE( aks3s(jpoce) )   ;  ALLOCATE( akf3s(jpoce) )    ;  ALLOCATE( sulfats(jpoce) ) 
    187       ALLOCATE( fluorids(jpoce) ) 
     192      ALLOCATE( fluorids(jpoce) ) ; ALLOCATE( akh2s(jpoce) )    ;  ALLOCATE( aknh3(jpoce) ) 
    188193 
    189194      akbs  (:) = 0. ;   ak1s   (:) = 0. ;  ak2s  (:) = 0. ;   akws   (:) = 0. 
    190195      ak1ps (:) = 0. ;   ak2ps  (:) = 0. ;  ak3ps (:) = 0. ;   aksis  (:) = 0. 
    191196      aksps (:) = 0. ;   ak12s  (:) = 0. ;  ak12ps(:) = 0. ;   ak123ps(:) = 0. 
    192       borats(:) = 0. ;   calcon2(:) = 0. ;  sieqs (:) = 0. 
     197      borats(:) = 0. ;   calcon2(:) = 0. ;  sieqs (:) = 0. ;   akh2s  (:) = 0. 
    193198      aks3s(:)  = 0. ;   akf3s(:)   = 0. ;  sulfats(:) = 0. ;  fluorids(:) = 0. 
     199      aknh3(:)  = 0. 
    194200 
    195201      ! Mass balance calculation   
    196       ALLOCATE( fromsed(jpoce, jpsol) ) ; ALLOCATE( tosed(jpoce, jpsol) ) ;  ALLOCATE( rloss(jpoce, jpsol) ) 
    197       ALLOCATE( tokbot (jpoce, jpwat) )  
    198  
    199       fromsed(:,:) = 0.    ;   tosed(:,:) = 0. ;  rloss(:,:) = 0.  ;   tokbot(:,:) = 0.  
     202      ALLOCATE( fromsed(jpoce, jpsol+jpads) ) ; ALLOCATE( tosed(jpoce, jpsol+jpads) ) 
     203 
     204      fromsed(:,:) = 0.    ;   tosed(:,:) = 0. 
    200205 
    201206      ! Initialization of sediment geometry 
    202207      !------------------------------------ 
    203       CALL sed_init_geom 
     208      CALL sed_ini_geom 
    204209 
    205210      ! Offline specific mode 
    206211      ! --------------------- 
    207212      ln_sediment_offline = .FALSE. 
     213 
     214      !--------------------------------------------- 
     215      ! Molecular weight [g/mol] for solid species 
     216      !--------------------------------------------- 
     217 
     218      ! opal=sio2*0.4(h20)=28+2*16+0.4*(2+16) 
     219      !--------------------------------------- 
     220      mol_wgt(jsopal) = 28. + 2. * 16. + 0.4 * ( 2. + 16. ) 
     221 
     222      !  clay 
     223      !  some kind of Illit (according to Pape) 
     224      !  K0.58(Al 1.38 Fe(III)0.37Fe(II)0.04Mg0.34)[(OH)2|(Si3.41Al0.59)O10] 
     225      !-------------------------------------------------------------------- 
     226      mol_wgt(jsclay) = 0.58 * 39. + 1.38 * 27. + ( 0.37 + 0.04 ) * 56.+ & 
     227         &              0.34 * 24. + 2. * ( 16. + 1. ) + 3.41 * 38. +    & 
     228         &              0.59 * 27. + 10. * 16. 
     229 
     230      mol_wgt(jsfeo)  = 55.0 + 3.0 * ( 16.0 + 1.0) 
     231 
     232 
     233      mol_wgt(jsfes)  = 55.0 + 32.0 
     234 
     235      ! for chemistry Poc : C(122)H(244)O(86)N(16)P(1) 
     236      ! But den sity of Poc is an Hydrated material (= POC + 30H2O) 
     237      ! So C(122)H(355)O(120)N(16)P(1) 
     238      !------------------------------------------------------------ 
     239      mol_wgt(jspoc) = ( redC * 12. + 355. + 120. * 16. + redNo3 * 14. + 31. ) / redC 
     240      mol_wgt(jspos) = mol_wgt(jspoc) 
     241      mol_wgt(jspor) = mol_wgt(jspoc) 
     242 
     243      ! CaCO3 
     244      !--------- 
     245      mol_wgt(jscal) = 40. + 12. + 3. * 16. 
     246 
     247      ! Density of solid material in sediment [g/cm**3] 
     248      !------------------------------------------------ 
     249      ALLOCATE( dens_sol(jpsol) ) 
     250      dens_sol(jsclay) = 2.6 
     251      dens_sol(jscal)  = 2.7 
     252      dens_sol(jsopal) = 2.1  
     253      dens_sol(jsfeo)  = 3.4 
     254      dens_sol(jsfes)  = 4.8 
     255      dens_sol(jspoc)  = 1.0 
     256      dens_sol(jspos)  = 1.0 
     257      dens_sol(jspor)  = 1.0 
     258 
     259      ! Accumulation rate from Burwicz et al. (2011). This is used to 
     260      ! compute the flux of clays and minerals 
     261      ! -------------------------------------------------------------- 
     262      DO ji = 1, jpoce 
     263          ztmp1 = 0.0117 * 3.0 / ( 1.0 + ( zkbot(ji) / 200.)**3 ) 
     264          ztmp2 = 0.006 / ( 1.0 + ( zkbot(ji) / 5000.)**10 ) 
     265          wacc(ji) = ztmp2+ztmp1 
     266      END DO 
     267 
     268      ! Vertical profile of of the adsorption factor for adsorbed species 
     269      ! ----------------------------------------------------------------- 
     270      radssol(:,jwfe2) = 1.0 / ( 1.0 + adsfe2 * por1(:) / por(:) ) 
     271      radssol(:,jwnh4) = 1.0 / ( 1.0 + adsnh4 * por1(:) / por(:) ) 
     272      rads1sol(:,jwnh4) = adsnh4 * por1(:) / ( por(:) + adsnh4 * por1(:) ) 
     273      rads1sol(:,jwfe2) = adsfe2 * por1(:) / ( por(:) + adsfe2 * por1(:) ) 
     274 
     275      ! Initialization of the array for non linear solving 
     276      ! -------------------------------------------------- 
     277 
     278      ALLOCATE( jarr(jpksed*jpvode,2) ) 
     279      ALLOCATE( jsvode(jpvode), isvode(jptrased) ) 
     280      jsvode(1) = jwoxy ; jsvode(2) = jwno3 ; jsvode(3) = jwnh4 
     281      jsvode(4) = jwh2s ; jsvode(5) = jwso4 ; jsvode(6) = jwfe2 
     282      jsvode(7) = jpwat+jsfeo ; jsvode(8) = jpwat+jsfes 
     283      isvode(jwoxy) = 1 ; isvode(jwno3) = 2 ; isvode(jwnh4) = 3 
     284      isvode(jwh2s) = 4 ; isvode(jwso4) = 5 ; isvode(jwfe2) = 6 
     285      isvode(jpwat+jsfeo) = 7 ; isvode(jpwat+jsfes) = 8 
     286      DO js = 1, jpvode 
     287         DO jk = 1, jpksed 
     288            jn = (jk-1) * jpvode + js 
     289            jarr(jn,1) = jk 
     290            jarr(jn,2) = jsvode(js) 
     291         END DO 
     292      END DO 
     293 
     294      ALLOCATE( rstepros(jpoce) ) 
    208295 
    209296#if defined key_sed_off 
     
    218305      ! -------------------------------- 
    219306      IF (ln_sediment_offline) THEN 
    220          CALL trc_dta_ini(jptra) 
    221307         CALL trc_dmp_sed_ini 
    222308      ENDIF 
    223309 
    224    END SUBROUTINE sed_init 
    225  
    226    SUBROUTINE sed_init_geom 
     310   END SUBROUTINE sed_ini 
     311 
     312   SUBROUTINE sed_ini_geom 
    227313      !!---------------------------------------------------------------------- 
    228       !!                   ***  ROUTINE sed_init_geom  *** 
     314      !!                   ***  ROUTINE sed_ini_geom  *** 
    229315      !! 
    230316      !! ** Purpose :  Initialization of sediment geometry 
     
    244330      !----------------------------------------------------------  
    245331 
    246       IF(lwp) WRITE(numsed,*) ' sed_init_geom : Initialization of sediment geometry ' 
     332      IF(lwp) WRITE(numsed,*) ' sed_ini_geom : Initialization of sediment geometry ' 
    247333      IF(lwp) WRITE(numsed,*) ' ' 
    248334 
    249335      ! Computation of 1D array of sediments points 
    250336      indoce = 0 
    251       DO_2D( 1, 1, 1, 1 ) 
     337      DO_2D( 0, 0, 0, 0 ) 
    252338         IF (  epkbot(ji,jj) > 0. ) THEN 
    253339            indoce          = indoce + 1 
     
    272358      !------------------------------------------------     
    273359      CALL pack_arr ( jpoce, dzkbot(1:jpoce), epkbot(1:jpi,1:jpj), iarroce(1:jpoce) ) 
     360#if defined key_sed_off 
     361      dzkbot(1:jpoce) = 1.E8 
     362#else 
    274363      dzkbot(1:jpoce) = dzkbot(1:jpoce) * 1.e+2  
     364#endif 
    275365      CALL pack_arr ( jpoce, zkbot(1:jpoce), gdepbot(1:jpi,1:jpj), iarroce(1:jpoce) ) 
     366      CALL pack_arr ( jpoce, slatit(1:jpoce), gphit(1:jpi,1:jpj), iarroce(1:jpoce) ) 
     367      CALL pack_arr ( jpoce, slongit(1:jpoce), glamt(1:jpi,1:jpj), iarroce(1:jpoce) ) 
    276368 
    277369      ! Geometry and  constants  
     
    285377      zsur = - za0 - za1 * sedacr * LOG( COSH( (1-sedkth) / sedacr )  ) 
    286378 
     379      dz(1)       = 0.1 
    287380      profsedw(1) = 0.0 
    288       profsed(1) = -dz(1) / 2. 
     381      profsed(1)  = -dz(1) / 2. 
    289382      DO jk = 2, jpksed 
    290383         zw = REAL( jk , wp ) 
     
    292385         profsed(jk)  = ( zsur + za0 * zt + za1 * sedacr * LOG ( COSH( (zt-sedkth) / sedacr ) )  )  
    293386         profsedw(jk) = ( zsur + za0 * zw + za1 * sedacr * LOG ( COSH( (zw-sedkth) / sedacr ) )  ) 
    294       END DO 
    295  
    296       dz(1) = 0.1 
    297       DO jk = 2, jpksed 
    298387         dz(jk) = profsedw(jk) - profsedw(jk-1) 
    299388      END DO 
    300389 
    301       DO jk = 1, jpksed 
    302          DO ji = 1, jpoce 
    303             dz3d(ji,jk) = dz(jk) 
    304          END DO 
    305       ENDDO 
     390      DO ji = 1, jpoce 
     391         dz3d(ji,:) = dz(:) 
     392      END DO 
    306393 
    307394      !  Porosity profile [0] 
     
    309396      por(1) = 1.0 
    310397      DO jk = 2, jpksed 
    311          por(jk) = porinf + ( porsurf-porinf) * exp(-rhox * (profsed(jk) ) ) 
     398         por(jk) = porinf + ( porsurf-porinf) * exp(-rhox * profsed(jk) ) 
    312399      END DO 
    313400  
     
    333420      dz3d(1:jpoce,1) = dz(1) 
    334421 
    335       !--------------------------------------------- 
    336       ! Molecular weight [g/mol] for solid species 
    337       !--------------------------------------------- 
    338  
    339       ! opal=sio2*0.4(h20)=28+2*16+0.4*(2+16) 
    340       !--------------------------------------- 
    341       mol_wgt(jsopal) = 28. + 2. * 16. + 0.4 * ( 2. + 16. )   
    342  
    343       !  clay 
    344       !  some kind of Illit (according to Pape) 
    345       !  K0.58(Al 1.38 Fe(III)0.37Fe(II)0.04Mg0.34)[(OH)2|(Si3.41Al0.59)O10] 
    346       !-------------------------------------------------------------------- 
    347       mol_wgt(jsclay) = 0.58 * 39. + 1.38 * 27. + ( 0.37 + 0.04 ) * 56.+ & 
    348          &              0.34 * 24. + 2. * ( 16. + 1. ) + 3.41 * 38. +    & 
    349          &              0.59 * 27. + 10. * 16. 
    350  
    351       mol_wgt(jsfeo)  = 55.0 + 3.0 * ( 16.0 + 1.0) 
    352  
    353       mol_wgt(jsfes)  = 55.0 + 32.0 
    354  
    355       ! for chemistry Poc : C(122)H(244)O(86)N(16)P(1) 
    356       ! But den sity of Poc is an Hydrated material (= POC + 30H2O) 
    357       ! So C(122)H(355)O(120)N(16)P(1) 
    358       !------------------------------------------------------------ 
    359       mol_wgt(jspoc) = ( 122. * 12. + 355. + 120. * 16.+ & 
    360          &                16. * 14. + 31. ) / 122. 
    361       mol_wgt(jspos) = mol_wgt(jspoc) 
    362       mol_wgt(jspor) = mol_wgt(jspoc) 
    363  
    364       ! CaCO3 
    365       !--------- 
    366       mol_wgt(jscal) = 40. + 12. + 3. * 16. 
    367  
    368       ! Density of solid material in sediment [g/cm**3] 
    369       !------------------------------------------------ 
    370       denssol = 2.6 
    371  
    372       ! Initialization of diffusion coefficient as function of porosity [cm**2/s] 
    373       !-------------------------------------------------------------------- 
    374 !      DO jn = 1, jpsol 
    375 !         DO jk = 1, jpksed 
    376 !            DO ji = 1, jpoce 
    377 !               diff(ji,jk,jn) = dcoef / ( 1.0 - 2.0 * log(por(jk)) ) 
    378 !            END DO 
    379 !         END DO 
    380 !      END DO 
    381  
    382       ! Accumulation rate from Burwicz et al. (2011). This is used to 
    383       ! compute the flux of clays and minerals 
    384       ! -------------------------------------------------------------- 
    385       DO ji = 1, jpoce 
    386           ztmp1 = 0.117 / ( 1.0 + ( zkbot(ji) / 200.)**3 ) 
    387           ztmp2 = 0.006 / ( 1.0 + ( zkbot(ji) / 4000.)**10 ) 
    388           wacc(ji) = ztmp1 + ztmp2 
    389       END DO 
    390  
    391  
    392       ! Initialization of time step as function of porosity [cm**2/s] 
    393       !------------------------------------------------------------------ 
    394    END SUBROUTINE sed_init_geom 
    395  
    396    SUBROUTINE sed_init_nam 
     422   END SUBROUTINE sed_ini_geom 
     423 
     424   SUBROUTINE sed_ini_nam 
    397425      !!---------------------------------------------------------------------- 
    398       !!                   ***  ROUTINE sed_init_nam  *** 
     426      !!                   ***  ROUTINE sed_ini_nam  *** 
    399427      !! 
    400428      !! ** Purpose :  Initialization of sediment geometry 
     
    405433      !!---------------------------------------------------------------------- 
    406434 
    407       CHARACTER(:), ALLOCATABLE ::   numnamsed_ref           !! Character buffer for reference namelist sediment 
    408       CHARACTER(:), ALLOCATABLE ::   numnamsed_cfg           !! Character buffer for configuration namelist sediment 
     435      INTEGER ::   numnamsed_ref = -1           !! Logical units for namelist sediment 
     436      INTEGER ::   numnamsed_cfg = -1           !! Logical units for namelist sediment 
    409437      INTEGER :: ios                 ! Local integer output status for namelist read 
    410438      CHARACTER(LEN=20)   ::   clname 
     
    421449      TYPE(PSED) , DIMENSION(jpdia2dsed) :: seddiag2d 
    422450 
    423       NAMELIST/nam_run/nrseddt,ln_sed_2way 
     451      NAMELIST/nam_run/ln_sed_2way,nrosorder,rosatol,rosrtol 
    424452      NAMELIST/nam_geom/jpksed, sedzmin, sedhmax, sedkth, sedacr, porsurf, porinf, rhox 
    425453      NAMELIST/nam_trased/sedsol, sedwat 
    426454      NAMELIST/nam_diased/seddiag3d, seddiag2d 
    427       NAMELIST/nam_inorg/rcopal, dcoef, rccal, ratligc, rcligc 
     455      NAMELIST/nam_inorg/rcopal, rccal, ratligc, rcligc 
    428456      NAMELIST/nam_poc/redO2, redNo3, redPo4, redC, redfep, rcorgl, rcorgs,  & 
    429          &             rcorgr, rcnh4, rch2s, rcfe2, rcfeh2s, rcfes, rcfeso, & 
    430          &             xksedo2, xksedno3, xksedfeo, xksedso4 
    431       NAMELIST/nam_btb/dbiot, ln_btbz, dbtbzsc, adsnh4, ln_irrig, xirrzsc 
     457         &             rcorgr, rcnh4, rch2s, rcfe2, rcfeh2s, rcfeso, rcfesp, & 
     458         &             rcfesd, xksedo2, xksedno3, xksedfeo, xksedso4 
     459      NAMELIST/nam_btb/dbiot, ln_btbz, dbtbzsc, adsnh4, adsfe2, ln_irrig, xirrzsc 
    432460      NAMELIST/nam_rst/ln_rst_sed, cn_sedrst_indir, cn_sedrst_outdir, cn_sedrst_in, cn_sedrst_out 
    433461 
     
    435463      !------------------------------------------------------- 
    436464 
    437       IF(lwp) WRITE(numsed,*) ' sed_init_nam : Read namelists ' 
     465      IF(lwp) WRITE(numsed,*) ' sed_ini_nam : Read namelists ' 
    438466      IF(lwp) WRITE(numsed,*) ' ' 
    439467 
     
    449477      !--------------------------------- 
    450478      clname = 'namelist_sediment' 
    451       IF(lwp) WRITE(numsed,*) ' sed_init_nam : read SEDIMENT namelist' 
     479      IF(lwp) WRITE(numsed,*) ' sed_ini_nam : read SEDIMENT namelist' 
    452480      IF(lwp) WRITE(numsed,*) ' ~~~~~~~~~~~~~~' 
    453       CALL load_nml( numnamsed_ref, TRIM( clname )//'_ref', numout, lwm ) 
    454       CALL load_nml( numnamsed_cfg, TRIM( clname )//'_cfg', numout, lwm ) 
     481      CALL ctl_opn( numnamsed_ref, TRIM( clname )//'_ref', 'OLD'    , 'FORMATTED', 'SEQUENTIAL', -1, numout, .FALSE. ) 
     482      CALL ctl_opn( numnamsed_cfg, TRIM( clname )//'_cfg', 'OLD'    , 'FORMATTED', 'SEQUENTIAL', -1, numout, .FALSE. ) 
    455483 
    456484      nitsed000 = nittrc000 
    457485      nitsedend = nitend 
    458486      ! Namelist nam_run 
     487      REWIND( numnamsed_ref )              ! Namelist nam_run in reference namelist : Pisces variables 
    459488      READ  ( numnamsed_ref, nam_run, IOSTAT = ios, ERR = 901) 
    460489901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_run in reference namelist' ) 
    461490 
     491      REWIND( numnamsed_cfg )              ! Namelist nam_run in reference namelist : Pisces variables 
    462492      READ  ( numnamsed_cfg, nam_run, IOSTAT = ios, ERR = 902) 
    463493902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_run in configuration namelist' ) 
     
    465495      IF (lwp) THEN 
    466496         WRITE(numsed,*) ' namelist nam_run' 
    467          WRITE(numsed,*) ' Nb of iterations for fast species    nrseddt = ', nrseddt 
    468497         WRITE(numsed,*) ' 2-way coupling between PISCES and Sed ln_sed_2way = ', ln_sed_2way 
     498         WRITE(numsed,*) ' Order of the Rosenbrock method (2,3,4) = ', nrosorder 
     499         WRITE(numsed,*) ' Tolerance for absolute error = ', rosatol 
     500         WRITE(numsed,*) ' Tolerance for relative order = ', rosrtol 
    469501      ENDIF 
    470502 
    471503      IF ( ln_p5z .AND. ln_sed_2way ) CALL ctl_stop( '2 ways coupling with sediment cannot be activated with PISCES-QUOTA' ) 
    472504 
     505      REWIND( numnamsed_ref )              ! Namelist nam_geom in reference namelist : Pisces variables 
    473506      READ  ( numnamsed_ref, nam_geom, IOSTAT = ios, ERR = 903) 
    474507903   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_geom in reference namelist' ) 
    475508 
     509      REWIND( numnamsed_cfg )              ! Namelist nam_geom in reference namelist : Pisces variables 
    476510      READ  ( numnamsed_cfg, nam_geom, IOSTAT = ios, ERR = 904) 
    477511904   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_geom in configuration namelist' ) 
     
    492526      dtsed = rDt_trc 
    493527 
     528      REWIND( numnamsed_ref )              ! Namelist nam_trased in reference namelist : Pisces variables 
    494529      READ  ( numnamsed_ref, nam_trased, IOSTAT = ios, ERR = 905) 
    495530905   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_trased in reference namelist' ) 
    496531 
     532      REWIND( numnamsed_cfg )              ! Namelist nam_trased in reference namelist : Pisces variables 
    497533      READ  ( numnamsed_cfg, nam_trased, IOSTAT = ios, ERR = 906) 
    498534906   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_trased in configuration namelist' ) 
     
    523559      ENDIF 
    524560 
     561      REWIND( numnamsed_ref )              ! Namelist nam_diased in reference namelist : Pisces variables 
    525562      READ  ( numnamsed_ref, nam_diased, IOSTAT = ios, ERR = 907) 
    526563907   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_diased in reference namelist' ) 
    527564 
     565      REWIND( numnamsed_cfg )              ! Namelist nam_diased in reference namelist : Pisces variables 
    528566      READ  ( numnamsed_cfg, nam_diased, IOSTAT = ios, ERR = 908) 
    529567908   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_diased in configuration namelist' ) 
     
    563601      ! Inorganic chemistry parameters 
    564602      !---------------------------------- 
     603      REWIND( numnamsed_ref )              ! Namelist nam_inorg in reference namelist : Pisces variables 
    565604      READ  ( numnamsed_ref, nam_inorg, IOSTAT = ios, ERR = 909) 
    566605909   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_inorg in reference namelist' ) 
    567606 
     607      REWIND( numnamsed_cfg )              ! Namelist nam_inorg in reference namelist : Pisces variables 
    568608      READ  ( numnamsed_cfg, nam_inorg, IOSTAT = ios, ERR = 910) 
    569609910   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_inorg in configuration namelist' ) 
     
    572612         WRITE(numsed,*) ' namelist nam_inorg' 
    573613         WRITE(numsed,*) ' reactivity for Si      rcopal  = ', rcopal 
    574          WRITE(numsed,*) ' diff. coef for por.    dcoef   = ', dcoef 
    575614         WRITE(numsed,*) ' reactivity for calcite rccal   = ', rccal 
    576615         WRITE(numsed,*) ' L/C ratio in POC       ratligc = ', ratligc 
     
    587626      ! Additional parameter linked to POC/O2/No3/Po4 
    588627      !---------------------------------------------- 
     628      REWIND( numnamsed_ref )              ! Namelist nam_poc in reference namelist : Pisces variables 
    589629      READ  ( numnamsed_ref, nam_poc, IOSTAT = ios, ERR = 911) 
    590630911   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_poc in reference namelist' ) 
    591631 
     632      REWIND( numnamsed_cfg )              ! Namelist nam_poc in reference namelist : Pisces variables 
    592633      READ  ( numnamsed_cfg, nam_poc, IOSTAT = ios, ERR = 912) 
    593634912   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_poc in configuration namelist' ) 
     
    607648         WRITE(numsed,*) ' reactivity for Fe2+              rcfe2    = ', rcfe2 
    608649         WRITE(numsed,*) ' reactivity for FeOH/H2S          rcfeh2s  = ', rcfeh2s 
    609          WRITE(numsed,*) ' reactivity for Fe2+/H2S          rcfes    = ', rcfes 
    610650         WRITE(numsed,*) ' reactivity for FeS/O2            rcfeso   = ', rcfeso 
     651         WRITE(numsed,*) ' Precipitation of FeS             rcfesp   = ', rcfesp 
     652         WRITE(numsed,*) ' Dissolution of FeS               rcfesd   = ', rcfesd 
    611653         WRITE(numsed,*) ' Half-sat. cste for oxic remin    xksedo2  = ', xksedo2 
    612654         WRITE(numsed,*) ' Half-sat. cste for denit.        xksedno3 = ', xksedno3 
     
    629671      reac_fe2   = rcfe2  / ryear 
    630672      reac_feh2s = rcfeh2s/ ryear 
    631       reac_fes   = rcfes  / ryear 
    632673      reac_feso  = rcfeso / ryear 
     674      reac_fesp  = rcfesp / ryear 
     675      reac_fesd  = rcfesd / ryear 
    633676 
    634677      ! reactivity rc in  [l.mol-1.s-1]       
     
    637680      ! Bioturbation parameter 
    638681      !------------------------ 
     682      REWIND( numnamsed_ref )              ! Namelist nam_btb in reference namelist : Pisces variables 
    639683      READ  ( numnamsed_ref, nam_btb, IOSTAT = ios, ERR = 913) 
    640684913   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_btb in reference namelist' ) 
    641685 
     686      REWIND( numnamsed_cfg )              ! Namelist nam_btb in reference namelist : Pisces variables 
    642687      READ  ( numnamsed_cfg, nam_btb, IOSTAT = ios, ERR = 914) 
    643688914   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_btb in configuration namelist' ) 
     
    649694         WRITE(numsed,*) ' coefficient for btb attenuation   dbtbzsc  = ', dbtbzsc 
    650695         WRITE(numsed,*) ' Adsorption coefficient of NH4     adsnh4   = ', adsnh4 
     696         WRITE(numsed,*) ' Adsorption coefficient of Fe2     adsfe2   = ', adsfe2 
    651697         WRITE(numsed,*) ' Bioirrigation in sediment         ln_irrig = ', ln_irrig 
    652698         WRITE(numsed,*) ' coefficient for irrig attenuation xirrzsc  = ', xirrzsc 
     
    656702      ! Initial value (t=0) for sediment pore water and solid components 
    657703      !---------------------------------------------------------------- 
     704      REWIND( numnamsed_ref )              ! Namelist nam_rst in reference namelist : Pisces variables 
    658705      READ  ( numnamsed_ref, nam_rst, IOSTAT = ios, ERR = 915) 
    659706915   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_rst in reference namelist' ) 
    660707 
     708      REWIND( numnamsed_cfg )              ! Namelist nam_rst in reference namelist : Pisces variables 
    661709      READ  ( numnamsed_cfg, nam_rst, IOSTAT = ios, ERR = 916) 
    662710916   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_rst in configuration namelist' ) 
     
    667715         WRITE(numsed,*) ' ' 
    668716      ENDIF 
    669       nn_dtsed = 1 
    670  
    671  
    672    END SUBROUTINE sed_init_nam 
     717 
     718      CLOSE( numnamsed_cfg ) 
     719      CLOSE( numnamsed_ref ) 
     720 
     721   END SUBROUTINE sed_ini_nam 
    673722 
    674723END MODULE sedini 
  • NEMO/trunk/src/TOP/PISCES/SED/sedinitrc.F90

    r12377 r15450  
    1010   !! * Modules used 
    1111   USE sed     ! sediment global variable 
    12    USE sed_oce 
    1312   USE sedini 
    1413   USE seddta 
     
    1716   USE sedchem 
    1817   USE sedarr 
     18   USE sed_oce 
    1919   USE lib_mpp         ! distribued memory computing library 
    2020 
     
    5050      !!        !  06-07  (C. Ethe)  Re-organization 
    5151      !!---------------------------------------------------------------------- 
    52       INTEGER, INTENT(in)  ::  Kbb, Kmm      ! time level indices 
     52      INTEGER, INTENT(in) ::   Kbb, Kmm  ! time level indices 
     53 
    5354      INTEGER :: ji, jj, ikt 
    5455      !!---------------------------------------------------------------------- 
     
    8687      !!        !  06-07  (C. Ethe)  original 
    8788      !!---------------------------------------------------------------------- 
    88       INTEGER, INTENT(in)  ::  Kbb, Kmm      ! time level indices 
    8989  
     90      INTEGER, INTENT(in) ::   Kbb, Kmm  ! time level indices 
     91 
    9092      ! local variables 
    91       INTEGER :: & 
    92          ji, jk, zhipor 
     93      INTEGER :: ji, jk, zhipor 
    9394 
    9495      !-------------------------------------------------------------------- 
     
    134135      ! Initialization of chemical constants 
    135136      CALL sed_chem ( nitsed000 ) 
    136  
    137       ! Stores initial sediment data for mass balance calculation 
    138       pwcp0 (1:jpoce,1:jpksed,1:jpwat ) = pwcp (1:jpoce,1:jpksed,1:jpwat )  
    139       solcp0(1:jpoce,1:jpksed,1:jpsol ) = solcp(1:jpoce,1:jpksed,1:jpsol)  
    140137 
    141138      ! Conversion of [h+] in mol/Kg to get it in mol/l ( multiplication by density) 
  • NEMO/trunk/src/TOP/PISCES/SED/sedinorg.F90

    r12839 r15450  
    88   USE sed     ! sediment global variable 
    99   USE sed_oce 
    10    USE sedmat  ! linear system of equations 
    11    USE sedco3  ! carbonate ion and proton concentration  
    1210   USE sedini 
    13    USE seddsr 
     11   USE sedmat 
    1412   USE lib_mpp         ! distribued memory computing library 
    1513   USE lib_fortran 
     
    2321CONTAINS 
    2422    
    25    SUBROUTINE sed_inorg( kt )  
     23   SUBROUTINE sed_inorg( kt ) 
    2624      !!---------------------------------------------------------------------- 
    2725      !!                   ***  ROUTINE sed_inorg  *** 
     
    4644      !!---------------------------------------------------------------------- 
    4745      !! Arguments 
    48       INTEGER, INTENT(in) ::   kt       ! number of iteration 
     46      INTEGER, INTENT(in)  :: kt   ! time step 
    4947      ! --- local variables 
    50       INTEGER :: ji, jk, js, jw         ! dummy looop indices 
    51       REAL(wp), DIMENSION(jpoce,jpksed) :: zrearat1, zrearat2    ! reaction rate in pore water 
    52       REAL(wp), DIMENSION(jpoce,jpksed) :: zundsat    ! undersaturation ; indice jpwatp1 is for calcite    
    53       REAL(wp), DIMENSION(jpoce) :: zco3eq 
    54       REAL(wp), DIMENSION(jpoce,jpksed,jpsol) :: zvolc    ! temp. variables 
    55       REAL(wp), DIMENSION(jpoce) :: zsieq 
    56       REAL(wp)  ::  zsolid1, zvolw, zreasat 
    57       REAL(wp)  ::  zsatur, zsatur2, znusil, zsolcpcl, zsolcpsi 
     48      INTEGER   ::  ji,jk          ! dummy looop indices 
     49      REAL(wp), DIMENSION(jpoce) ::  zsieq, reac_silf 
     50      REAL(wp)  ::  zsolid1, zreasat, zco3sat 
     51      REAL(wp)  ::  zsatur, zsatur2, znusil, zsolcpcl, zsolcpsi, zexcess 
     52      REAL(wp), DIMENSION(jpoce,jpksed) :: zundsat, zrearat, psms 
    5853      !! 
    5954      !!---------------------------------------------------------------------- 
    6055 
    6156      IF( ln_timing )  CALL timing_start('sed_inorg') 
     57 
     58      IF( kt == nitsed000 ) THEN 
     59         IF (lwp) WRITE(numsed,*) ' sed_inorg : Dissolution of CaCO3 and BSi  ' 
     60         IF (lwp) WRITE(numsed,*) ' ' 
     61      ENDIF 
    6262! 
    63       IF( kt == nitsed000 ) THEN 
    64          IF (lwp) THEN 
    65             WRITE(numsed,*) ' sed_inorg : Dissolution reaction ' 
    66             WRITE(numsed,*) ' ' 
    67          ENDIF 
    68 !         !  
    69       ENDIF 
     63      DO ji = 1, jpoce 
     64         ! ----------------------------------------------- 
     65         ! Computation of Si solubility 
     66         ! Param of Ridgwell et al. 2002 
     67         ! ----------------------------------------------- 
    7068 
    71      ! Initializations 
    72      !---------------------- 
    73        
    74       zrearat1(:,:) = 0.    ;   zundsat(:,:)  = 0.  
    75       zrearat2(:,:) = 0.    ;   zrearat2(:,:) = 0. 
    76       zco3eq(:)     = rtrn 
    77       zvolc(:,:,:)  = 0. 
    78  
    79       ! ----------------------------------------------- 
    80       ! Computation of Si solubility 
    81       ! Param of Ridgwell et al. 2002 
    82       ! ----------------------------------------------- 
    83  
    84       DO ji = 1, jpoce 
    8569         zsolcpcl = 0.0 
    8670         zsolcpsi = 0.0 
    8771         DO jk = 1, jpksed 
    88             zsolcpsi = zsolcpsi + solcp(ji,jk,jsopal) * dz(jk) 
    89             zsolcpcl = zsolcpcl + solcp(ji,jk,jsclay) * dz(jk) 
     72            zsolcpsi = zsolcpsi + solcp(ji,jk,jsopal) * vols3d(ji,jk) 
     73            zsolcpcl = zsolcpcl + solcp(ji,jk,jsclay) * vols3d(ji,jk) 
    9074         END DO 
    9175         zsolcpsi = MAX( zsolcpsi, rtrn ) 
    92          zsieq(ji) = sieqs(ji) * MAX(0.25, 1.0 - (0.045 * zsolcpcl / zsolcpsi )**0.58 ) 
    93          zsieq(ji) = MAX( rtrn, sieqs(ji) ) 
     76         zsieq(ji) = sieqs(ji) * MAX(0.2, 1.0 - (0.045 * zsolcpcl / zsolcpsi )**0.58 ) 
     77         reac_silf(ji) = reac_sil * ( 0.05 + 0.055 * ( 1.64 * ( zsolcpcl / zsolcpsi + 0.01 ) )**(-0.75) ) / 1.25  
    9478      END DO 
    9579 
    96       DO js = 1, jpsol 
    97          DO jk = 1, jpksed 
    98             DO ji = 1, jpoce     
    99                zvolc(ji,jk,js) = ( vols3d(ji,jk) * dens_mol_wgt(js) ) /  & 
    100                   &              ( volw3d(ji,jk) * 1.e-3  )      
    101             ENDDO 
    102          ENDDO 
    103       ENDDO 
    104  
    105       !---------------------------------------------------------- 
    106       ! 5.  Beginning of  Pore Water diffusion and solid reaction 
    107       !--------------------------------------------------------- 
    10880       
    109       !----------------------------------------------------------------------------- 
    110       ! For jk=2,jpksed, and for couple  
    111       !  1 : jwsil/jsopal  ( SI/Opal ) 
    112       !  2 : jsclay/jsclay ( clay/clay )  
    113       !  3 : jwoxy/jspoc   ( O2/POC ) 
    114       !  reaction rate is a function of solid=concentration in solid reactif in [mol/l]  
    115       !  and undersaturation in [mol/l]. 
    116       !  Solid weight fractions should be in ie [mol/l]) 
    117       !  second member and solution are in zundsat variable 
    118       !------------------------------------------------------------------------- 
    119  
    120       DO jk = 1, jpksed 
    121          DO ji = 1, jpoce 
    122             ! For Silicic Acid and clay 
    123             zundsat(ji,jk) = zsieq(ji) - pwcp(ji,jk,jwsil) 
    124          ENDDO 
    125       ENDDO 
    126        
    127       ! Definition of reaction rates [rearat]=sans dim  
    128       ! For jk=1 no reaction (pure water without solid) for each solid compo 
    12981      DO ji = 1, jpoce 
    130          zrearat1(ji,:) = 0. 
    131          zrearat2(ji,:) = 0. 
    132       ENDDO 
    133  
    134       ! left hand side of coefficient matrix 
    135       DO jk = 2, jpksed 
    136          DO ji = 1, jpoce 
    137             zsolid1 = zvolc(ji,jk,jsopal) * solcp(ji,jk,jsopal) 
    138             zsatur = MAX(0., zundsat(ji,jk) / zsieq(ji) ) 
     82         DO jk = 2, jpksed 
     83            zsolid1 = volc(ji,jk,jsopal) * solcp(ji,jk,jsopal) 
     84            zsatur = MAX(0., ( zsieq(ji) - pwcp(ji,jk,jwsil) ) / zsieq(ji) ) 
    13985            zsatur2 = (1.0 + temp(ji) / 400.0 )**37 
    140             znusil = ( 0.225 * ( 1.0 + temp(ji) / 15.) + 0.775 * zsatur2 * zsatur**2.25 ) / zsieq(ji) 
    141             zrearat1(ji,jk)  = ( reac_sil * znusil * dtsed * zsolid1 ) / & 
    142                &                ( 1. + reac_sil * znusil * dtsed * zundsat(ji,jk) ) 
    143          ENDDO 
    144       ENDDO 
    145  
    146       CALL sed_mat( jwsil, jpoce, jpksed, zrearat1, zrearat2, zundsat, dtsed ) 
    147  
    148       ! New solid concentration values (jk=2 to jksed) for each couple  
    149       DO jk = 2, jpksed 
    150          DO ji = 1, jpoce 
    151             zreasat = zrearat1(ji,jk) * zundsat(ji,jk) / ( zvolc(ji,jk,jsopal) ) 
    152             solcp(ji,jk,jsopal) = solcp(ji,jk,jsopal) - zreasat 
    153          ENDDO 
    154       ENDDO 
    155  
    156       ! New pore water concentrations     
    157       DO jk = 1, jpksed 
    158          DO ji = 1, jpoce 
    159             ! Acid Silicic  
    160             pwcp(ji,jk,jwsil)  = zsieq(ji) - zundsat(ji,jk) 
    161          ENDDO 
    162       ENDDO 
     86            znusil = ( 0.225 * ( 1.0 + temp(ji) / 15.) * zsatur + 0.775 * zsatur2 * zsatur**9.25 ) 
     87            solcp(ji,jk,jsopal) = solcp(ji,jk,jsopal) - reac_silf(ji) * znusil * dtsed * solcp(ji,jk,jsopal) 
     88            pwcp(ji,jk,jwsil) = pwcp(ji,jk,jwsil) + reac_silf(ji) * znusil * dtsed * zsolid1 
     89         END DO 
     90      END DO 
    16391 
    16492      !--------------------------------------------------------------- 
     
    16795 
    16896      ! computes co3por from the updated pwcp concentrations (note [co3por] = mol/l) 
    169  
    170       CALL sed_co3( kt ) 
    171  
    17297      ! *densSW(l)**2 converts aksps [mol2/kg sol2] into [mol2/l2] to get [undsat] in [mol/l] 
    173       DO jk = 1, jpksed 
    174          DO ji = 1, jpoce 
    175             zco3eq(ji)     = aksps(ji) * densSW(ji) * densSW(ji) / ( calcon2(ji) + rtrn ) 
    176             zco3eq(ji)     = MAX( rtrn, zco3eq(ji) )  
    177             zundsat(ji,jk) = MAX(0., zco3eq(ji) - co3por(ji,jk) ) 
    178          ENDDO 
    179       ENDDO 
    180  
    181       DO jk = 2, jpksed 
    182          DO ji = 1, jpoce 
    183             zsolid1 = zvolc(ji,jk,jscal) * solcp(ji,jk,jscal) 
    184             zrearat1(ji,jk) = ( reac_cal * dtsed * zsolid1 / zco3eq(ji) ) / & 
    185                   &               ( 1. + reac_cal * dtsed * zundsat(ji,jk) / zco3eq(ji) ) 
     98      DO ji = 1, jpoce 
     99         zco3sat = aksps(ji) * densSW(ji) * densSW(ji) / ( calcon2(ji) + rtrn ) 
     100         saturco3(ji,:) = 1.0 - co3por(ji,:) / ( rtrn + zco3sat ) 
     101         DO jk = 1, jpksed 
     102            zsolid1 = volc(ji,jk,jscal) * solcp(ji,jk,jscal) 
     103            zundsat(ji,jk) = MAX( 0., zco3sat - co3por(ji,jk) ) 
     104            zrearat(ji,jk) = ( reac_cal * zsolid1 / ( zco3sat + rtrn ) ) / & 
     105            &                ( 1. + reac_cal * dtsed * zundsat(ji,jk) / ( zco3sat + rtrn ) ) 
    186106         END DO 
    187107      END DO 
    188108 
     109      psms(:,:) = 0.0 
    189110      ! solves tridiagonal system 
    190       CALL sed_mat( jwdic, jpoce, jpksed, zrearat1, zrearat2, zundsat, dtsed ) 
     111      CALL sed_mat_dsri( jpksed, jwdic, -zrearat(:,:), psms(:,:), dtsed, zundsat ) 
    191112 
    192113      ! New solid concentration values (jk=2 to jksed) for cacO3 
    193114      DO jk = 2, jpksed 
    194115         DO ji = 1, jpoce 
    195             zreasat = zrearat1(ji,jk) * zundsat(ji,jk) / zvolc(ji,jk,jscal) 
     116            zreasat = zrearat(ji,jk) * dtsed * zundsat(ji,jk) / ( volc(ji,jk,jscal) + rtrn ) 
    196117            solcp(ji,jk,jscal) = solcp(ji,jk,jscal) - zreasat 
     118            zreasat = zrearat(ji,jk) * dtsed * zundsat(ji,jk) 
     119            ! For DIC 
     120            pwcp(ji,jk,jwdic)  = pwcp(ji,jk,jwdic) + zreasat 
     121            ! For alkalinity 
     122            pwcp(ji,jk,jwalk)  = pwcp(ji,jk,jwalk) + 2.* zreasat 
    197123         ENDDO 
    198124      ENDDO 
    199125 
    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 
    296126 
    297127      IF( ln_timing )  CALL timing_stop('sed_inorg') 
  • NEMO/trunk/src/TOP/PISCES/SED/sedmat.F90

    r10222 r15450  
    1414   PRIVATE 
    1515 
    16    PUBLIC sed_mat  
    17  
    18    INTERFACE sed_mat 
    19       MODULE PROCEDURE sed_mat_dsr, sed_mat_btb 
    20    END INTERFACE 
    21  
    22    INTEGER, PARAMETER :: nmax = 30  
    23  
     16   PUBLIC sed_mat_dsr  
     17   PUBLIC sed_mat_dsrjac 
     18   PUBLIC sed_mat_dsri 
     19   PUBLIC sed_mat_btb 
     20   PUBLIC sed_mat_btbjac 
     21   PUBLIC sed_mat_btbi 
     22   PUBLIC sed_mat_coef 
    2423 
    2524   !! $Id$ 
    2625 CONTAINS 
    2726 
    28     SUBROUTINE sed_mat_dsr( nvar, ndim, nlev, preac, psms, psol, dtsed_in ) 
     27    SUBROUTINE sed_mat_coef( nksed ) 
    2928       !!--------------------------------------------------------------------- 
    30        !!                  ***  ROUTINE sed_mat_dsr  *** 
     29       !!                  ***  ROUTINE sed_mat_coef  *** 
    3130       !! 
    3231       !! ** Purpose :  solves tridiagonal system of linear equations  
     
    4948       !!---------------------------------------------------------------------- 
    5049       !! * Arguments 
    51        INTEGER , INTENT(in) ::  nvar  ! number of variable 
    52        INTEGER , INTENT(in) ::  ndim  ! number of points 
    53        INTEGER , INTENT(in) ::  nlev  ! number of sediment levels 
    54  
    55        REAL(wp), DIMENSION(ndim,nlev), INTENT(in   ) :: preac  ! reaction rates 
    56        REAL(wp), DIMENSION(ndim,nlev), INTENT(in   ) :: psms  ! reaction rates 
    57        REAL(wp), DIMENSION(ndim,nlev), INTENT(inout) :: psol   ! solution ( undersaturation values ) 
    58        REAL(wp), INTENT(in) ::  dtsed_in 
    59   
     50       INTEGER, INTENT(in) :: nksed 
     51 
    6052       !---Local declarations 
    61        INTEGER  ::  ji, jk, jn 
    62        REAL(wp), DIMENSION(ndim,nlev) :: za, zb, zc, zr 
    63        REAL(wp), DIMENSION(ndim)      :: zbet 
    64        REAL(wp), DIMENSION(ndim,nmax) :: zgamm 
    65  
    66        REAL(wp) ::  aplus,aminus   
    67        REAL(wp) ::  rplus,rminus    
    68        REAL(wp) ::  dxplus,dxminus 
    69  
     53       INTEGER  ::  ji, jk 
     54       REAL(wp) ::  aplus, aminus, dxplus, dxminus 
    7055       !---------------------------------------------------------------------- 
    7156 
    72        IF( ln_timing )  CALL timing_start('sed_mat_dsr') 
     57       IF( ln_timing )  CALL timing_start('sed_mat_coef') 
    7358 
    7459       ! Computation left hand side of linear system of  
    7560       ! equations for dissolution reaction 
    7661       !--------------------------------------------- 
    77  
    78  
     62       ! first sediment level           
     63       DO ji = 1, jpoce 
     64          aplus  = ( por(1) + por(2) ) * 0.5 
     65          dxplus = ( dz3d(ji,1) + dz3d(ji,2) ) / 2. 
     66          apluss(ji,1) = ( 1.0 / ( volw3d(ji,1) ) ) * aplus / dxplus 
     67 
     68          DO jk = 2, nksed - 1 
     69             aminus  = ( por(jk-1) + por(jk) ) * 0.5 
     70             dxminus = ( dz3d(ji,jk-1) + dz3d(ji,jk) ) / 2. 
     71 
     72             aplus   = ( por(jk+1) + por(jk) ) * 0.5 
     73             dxplus  = ( dz3d(ji,jk) + dz3d(ji,jk+1) ) / 2 
     74             ! 
     75             aminuss(ji,jk) = ( 1.0 / volw3d(ji,jk) ) * aminus / dxminus 
     76             apluss (ji,jk) = ( 1.0 / volw3d(ji,jk) ) * aplus / dxplus 
     77          END DO 
     78 
     79          aminus  = ( por(nksed-1) + por(nksed) ) * 0.5 
     80          dxminus = ( dz3d(ji,nksed-1) + dz3d(ji,nksed) ) / 2. 
     81          aminuss(ji,nksed)  = ( 1.0 / volw3d(ji,nksed) ) * aminus / dxminus 
     82 
     83       END DO 
     84 
     85       IF( ln_timing )  CALL timing_stop('sed_mat_coef') 
     86 
     87    END SUBROUTINE sed_mat_coef 
     88 
     89    SUBROUTINE sed_mat_dsr( nksed, nvar, accmask ) 
     90       !!--------------------------------------------------------------------- 
     91       !!                  ***  ROUTINE sed_mat_dsr  *** 
     92       !! 
     93       !! ** Purpose :  solves tridiagonal system of linear equations  
     94       !! 
     95       !! ** Method  :  
     96       !!        1 - computes left hand side of linear system of equations 
     97       !!            for dissolution reaction 
     98       !!         For mass balance in kbot+sediment : 
     99       !!              dz3d  (:,1) = dz(1) = 0.5 cm 
     100       !!              volw3d(:,1) = dzkbot ( see sedini.F90 )  
     101       !!              dz(2)       = 0.3 cm  
     102       !!              dz3d(:,2)   = 0.3 + dzdep   ( see seddsr.F90 )      
     103       !!              volw3d(:,2) and vols3d(l,2) are thickened ( see seddsr.F90 )  
     104       !! 
     105       !!         2 - forward/backward substitution.  
     106       !! 
     107       !!   History : 
     108       !!        !  04-10 (N. Emprin, M. Gehlen ) original 
     109       !!        !  06-04 (C. Ethe)  Module Re-organization 
     110       !!---------------------------------------------------------------------- 
     111       !! * Arguments 
     112       INTEGER , INTENT(in) ::  nvar, nksed  ! number of variable 
     113       INTEGER, DIMENSION(jpoce) :: accmask 
     114       INTEGER :: ji 
     115 
     116       !---Local declarations 
     117       INTEGER  ::  jk, jn 
     118       REAL(wp), DIMENSION(nksed) :: za, zb, zc 
     119 
     120       REAL(wp) ::  rplus,rminus    
     121       !---------------------------------------------------------------------- 
     122 
     123       IF( ln_timing )  CALL timing_start('sed_mat_dsr') 
     124 
     125       ! Computation left hand side of linear system of  
     126       ! equations for dissolution reaction 
     127       !--------------------------------------------- 
    79128       jn = nvar 
    80129       ! first sediment level           
    81        DO ji = 1, ndim 
    82           aplus  = ( ( volw3d(ji,1) / ( dz3d(ji,1) ) ) + & 
    83                         ( volw3d(ji,2) / ( dz3d(ji,2) ) ) ) / 2. 
    84           dxplus = ( dz3d(ji,1) + dz3d(ji,2) ) / 2. 
    85           rplus  = ( dtsed_in / ( volw3d(ji,1) ) ) * diff(ji,1,jn) * aplus / dxplus  
     130       DO ji = 1, jpoce 
     131          IF (accmask(ji) == 0) THEN 
     132             rplus  = apluss(ji,1) * diff(ji,1,jn) * radssol(1,jn) 
     133 
     134             za(1) = 0. 
     135             zb(1) = rplus 
     136             zc(1) = -rplus 
     137  
     138             DO jk = 2, nksed - 1 
     139                rminus  = aminuss(ji,jk) * diff(ji,jk-1,jn) * radssol(jk,jn) 
     140                rplus   = apluss (ji,jk) * diff(ji,jk,jn) * radssol(jk,jn) 
     141                !      
     142                za(jk) = -rminus 
     143                zb(jk) = rminus + rplus  
     144                zc(jk) = -rplus 
     145             END DO 
     146 
     147             rminus  = aminuss(ji,nksed) * diff(ji,nksed-1,jn) * radssol(nksed,jn) 
     148             ! 
     149             za(nksed) = -rminus 
     150             zb(nksed) = rminus 
     151             zc(nksed) = 0. 
     152 
     153             ! solves tridiagonal system of linear equations  
     154             ! ----------------------------------------------- 
     155 
     156             pwcpa(ji,1,jn) = pwcpa(ji,1,jn) - ( zc(1) * pwcp(ji,2,jn) + zb(1) * pwcp(ji,1,jn) ) 
     157             DO jk = 2, nksed - 1 
     158                pwcpa(ji,jk,jn) =  pwcpa(ji,jk,jn) - ( zc(jk) * pwcp(ji,jk+1,jn) + za(jk) * pwcp(ji,jk-1,jn)    & 
     159                &                  + zb(jk) * pwcp(ji,jk,jn) ) 
     160             ENDDO 
     161             pwcpa(ji,nksed,jn) = pwcpa(ji,nksed,jn) - ( za(nksed) * pwcp(ji,nksed-1,jn)    & 
     162             &                     + zb(nksed) * pwcp(ji,nksed,jn) ) 
     163 
     164          ENDIF 
     165       END DO 
     166 
     167       IF( ln_timing )  CALL timing_stop('sed_mat_dsr') 
     168 
     169    END SUBROUTINE sed_mat_dsr 
     170 
     171    SUBROUTINE sed_mat_dsrjac( nksed, nvar, NEQ, NROWPD, jacvode, accmask ) 
     172       !!--------------------------------------------------------------------- 
     173       !!                  ***  ROUTINE sed_mat_dsrjac  *** 
     174       !! 
     175       !! ** Purpose :  solves tridiagonal system of linear equations  
     176       !! 
     177       !! ** Method  :  
     178       !!        1 - computes left hand side of linear system of equations 
     179       !!            for dissolution reaction 
     180       !!         For mass balance in kbot+sediment : 
     181       !!              dz3d  (:,1) = dz(1) = 0.5 cm 
     182       !!              volw3d(:,1) = dzkbot ( see sedini.F90 )  
     183       !!              dz(2)       = 0.3 cm  
     184       !!              dz3d(:,2)   = 0.3 + dzdep   ( see seddsr.F90 )      
     185       !!              volw3d(:,2) and vols3d(l,2) are thickened ( see 
     186       !seddsr.F90 )  
     187       !! 
     188       !!         2 - forward/backward substitution.  
     189       !! 
     190       !!   History : 
     191       !!        !  04-10 (N. Emprin, M. Gehlen ) original 
     192       !!        !  06-04 (C. Ethe)  Module Re-organization 
     193       !!---------------------------------------------------------------------- 
     194       !! * Arguments 
     195       INTEGER , INTENT(in) ::  nvar, nksed, NEQ, NROWPD  ! number of variable 
     196       REAL, DIMENSION(jpoce,NROWPD,NEQ), INTENT(inout) :: jacvode 
     197       INTEGER, DIMENSION(jpoce), INTENT(in) :: accmask 
     198 
     199       !---Local declarations 
     200       INTEGER  ::  ji,jk, jn, jnn, jni, jnj ,jnij 
     201       REAL(wp), DIMENSION(nksed) :: za, zb, zc 
     202 
     203       REAL(wp) ::  rplus,rminus 
     204       !---------------------------------------------------------------------- 
     205 
     206       IF( ln_timing )  CALL timing_start('sed_mat_dsrjac') 
     207 
     208       ! Computation left hand side of linear system of  
     209       ! equations for dissolution reaction 
     210       !--------------------------------------------- 
     211       jn = nvar 
     212       ! first sediment level           
     213       DO ji = 1, jpoce 
     214          IF (accmask(ji) == 0 ) THEN 
     215             rplus  = apluss(ji,1) * diff(ji,1,jn) * radssol(1,jn) 
     216 
     217             za(1) = 0. 
     218             zb(1) = rplus 
     219             zc(1) = -rplus 
     220 
     221             DO jk = 2, nksed - 1 
     222                rminus  = aminuss(ji,jk) * diff(ji,jk-1,jn) * radssol(jk,jn) 
     223                rplus   = apluss (ji,jk) * diff(ji,jk,jn) * radssol(jk,jn) 
     224                !      
     225                za(jk) = -rminus 
     226                zb(jk) = rminus + rplus 
     227                zc(jk) = -rplus 
     228             END DO 
     229 
     230             rminus  = aminuss(ji,nksed) * diff(ji,nksed-1,jn) * radssol(nksed,jn) 
     231             ! 
     232             za(nksed) = -rminus 
     233             zb(nksed) = rminus 
     234             zc(nksed) = 0. 
     235 
     236             ! solves tridiagonal system of linear equations  
     237 
     238             jnn = isvode(jn) 
     239             jnij = jpvode + 1 
     240             jacvode(ji, jnij, jnn) = jacvode(ji,jnij,jnn) - zb(1) 
     241             jnj = jpvode + jnn 
     242             jnij = jnn - jnj + jpvode + 1 
     243             jacvode(ji, jnij, jnj) = jacvode(ji, jnij, jnj) -zc(1) 
     244             DO jk = 2, nksed - 1 
     245                jni = (jk-1) * jpvode + jnn 
     246                jnj = (jk-2) * jpvode + jnn 
     247                jnij = jni - jnj + jpvode + 1 
     248                jacvode(ji, jnij, jnj) = jacvode(ji, jnij, jnj) - za(jk) 
     249                jnj = (jk-1) * jpvode + jnn 
     250                jnij = jni - jnj + jpvode + 1 
     251                jacvode(ji, jnij, jnj) = jacvode(ji, jnij, jnj) - zb(jk) 
     252                jnj = (jk) * jpvode + jnn 
     253                jnij = jni - jnj + jpvode + 1 
     254                jacvode(ji, jnij, jnj)   = jacvode(ji, jnij, jnj) - zc(jk) 
     255             END DO 
     256             jni = (nksed-1) * jpvode + jnn 
     257             jnj = (nksed-2) * jpvode + jnn 
     258             jnij = jni - jnj + jpvode + 1 
     259             jacvode(ji, jnij, jnj) = jacvode(ji, jnij, jnj) - za(nksed) 
     260             jnij = jpvode + 1 
     261             jacvode(ji, jnij, jni) = jacvode(ji, jnij, jni) - zb(nksed) 
     262          ENDIF 
     263       END DO 
     264 
     265       IF( ln_timing )  CALL timing_stop('sed_mat_dsrjac') 
     266 
     267    END SUBROUTINE sed_mat_dsrjac 
     268 
     269    SUBROUTINE sed_mat_btbi( nksed, nvar, psol, preac, dtsed_in ) 
     270       !!--------------------------------------------------------------------- 
     271       !!                  ***  ROUTINE sed_mat_btb  *** 
     272       !! 
     273       !! ** Purpose :  solves tridiagonal system of linear equations  
     274       !! 
     275       !! ** Method  :  
     276       !!        1 - computes left hand side of linear system of equations 
     277       !!            for dissolution reaction 
     278       !! 
     279       !! 
     280       !!         2 - forward/backward substitution.  
     281       !! 
     282       !!   History : 
     283       !!        !  04-10 (N. Emprin, M. Gehlen ) original 
     284       !!        !  06-04 (C. Ethe)  Module Re-organization 
     285       !!---------------------------------------------------------------------- 
     286       !! * Arguments 
     287       INTEGER , INTENT(in) :: nksed, nvar      ! number of sediment levels 
     288 
     289      REAL(wp), DIMENSION(jpoce,nksed,nvar), INTENT(inout) :: & 
     290          psol, preac 
     291 
     292      REAL(wp), INTENT(in) :: dtsed_in 
     293 
     294       !---Local declarations 
     295       INTEGER  ::  & 
     296          ji, jk, jn 
     297 
     298       REAL(wp) ::  & 
     299          aplus,aminus   ,  & 
     300          rplus,rminus   ,  & 
     301          dxplus,dxminus 
     302 
     303       REAL(wp), DIMENSION(nksed)    :: za, zb, zc 
     304       REAL(wp), DIMENSION(nksed)    :: zr, zgamm 
     305       REAL(wp) ::  zbet 
     306 
     307       !---------------------------------------------------------------------- 
     308 
     309      ! Computation left hand side of linear system of  
     310      ! equations for dissolution reaction 
     311      !--------------------------------------------- 
     312      IF( ln_timing )  CALL timing_start('sed_mat_btbi') 
     313 
     314      ! first sediment level           
     315      DO ji = 1, jpoce 
     316         aplus  = ( por1(2) + por1(3) ) / 2.0 
     317         dxplus = ( dz(2) + dz(3) ) / 2. 
     318         rplus  = ( dtsed_in / vols(2) ) * db(ji,2) * aplus / dxplus 
     319         za(2) = 0. 
     320         zb(2) = 1. + rplus 
     321         zc(2) = -rplus 
     322 
     323         DO jk = 3, nksed - 1 
     324            aminus  = ( por1(jk-1) + por1(jk) ) * 0.5 
     325            aminus  = ( ( vols(jk-1) / dz(jk-1) ) + ( vols(jk) / dz(jk) ) ) / 2. 
     326            dxminus = ( dz(jk-1) + dz(jk) ) / 2. 
     327            rminus  = ( dtsed_in / vols(jk) ) * db(ji,jk-1) * aminus / dxminus 
     328            ! 
     329            aplus   = ( por1(jk) + por1(jk+1) ) * 0.5 
     330            dxplus  = ( dz(jk) + dz(jk+1) ) / 2. 
     331            rplus   = ( dtsed_in / vols(jk) ) * db(ji,jk) * aplus / dxplus 
     332            !      
     333            za(jk) = -rminus 
     334            zb(jk) = 1. + rminus + rplus 
     335            zc(jk) = -rplus 
     336 
     337         ENDDO 
     338 
     339         aminus = ( por1(nksed-1) + por1(nksed) ) * 0.5 
     340         dxminus = ( dz(nksed-1) + dz(nksed) ) / 2. 
     341         rminus  = ( dtsed_in / vols(nksed) ) * db(ji,nksed-1) * aminus / dxminus 
     342         ! 
     343         za(nksed) = -rminus 
     344         zb(nksed) = 1. + rminus 
     345         zc(nksed) = 0. 
     346 
     347         ! solves tridiagonal system of linear equations  
     348         ! -----------------------------------------------     
     349         DO jn = 1, nvar 
     350            zr(:) = psol(ji,:,jn) 
     351            zbet     = zb(2) - preac(ji,2,jn) * dtsed_in 
     352            psol(ji,2,jn) = zr(2) / zbet 
     353            !  
     354            DO jk = 3, nksed 
     355               zgamm(jk) =  zc(jk-1) / zbet 
     356               zbet      =  zb(jk) - preac(ji,jk,jn) * dtsed_in - za(jk) * zgamm(jk) 
     357               psol(ji,jk,jn) = ( zr(jk) - za(jk) * psol(ji,jk-1,jn) ) / zbet 
     358            ENDDO 
     359            !  
     360            DO jk = nksed - 1, 2, -1 
     361               psol(ji,jk,jn) = psol(ji,jk,jn) - zgamm(jk+1) * psol(ji,jk+1,jn) 
     362            ENDDO 
     363         END DO 
     364      END DO 
     365      ! 
     366      IF( ln_timing )  CALL timing_stop('sed_mat_btbi') 
     367 
     368    END SUBROUTINE sed_mat_btbi 
     369 
     370 
     371    SUBROUTINE sed_mat_btb( nksed, nvar, accmask ) 
     372       !!--------------------------------------------------------------------- 
     373       !!                  ***  ROUTINE sed_mat_btb  *** 
     374       !! 
     375       !! ** Purpose :  solves tridiagonal system of linear equations  
     376       !! 
     377       !! ** Method  :  
     378       !!        1 - computes left hand side of linear system of equations 
     379       !!            for dissolution reaction 
     380       !! 
     381       !!         2 - forward/backward substitution.  
     382       !! 
     383       !!   History : 
     384       !!        !  04-10 (N. Emprin, M. Gehlen ) original 
     385       !!        !  06-04 (C. Ethe)  Module Re-organization 
     386       !!---------------------------------------------------------------------- 
     387       !! * Arguments 
     388       INTEGER , INTENT(in) :: & 
     389          nvar, nksed     ! number of sediment levels 
     390       INTEGER, DIMENSION(jpoce) :: accmask 
     391 
     392       !---Local declarations 
     393       INTEGER  ::  ji, jk, jn 
     394 
     395       REAL(wp) ::  & 
     396          aplus,aminus   ,  & 
     397          rplus,rminus   ,  & 
     398          dxplus,dxminus 
     399 
     400       REAL(wp), DIMENSION(nksed)      :: za, zb, zc 
     401 
     402       !---------------------------------------------------------------------- 
     403 
     404      ! Computation left hand side of linear system of  
     405      ! equations for dissolution reaction 
     406      !--------------------------------------------- 
     407      IF( ln_timing )  CALL timing_start('sed_mat_btb') 
     408 
     409      ! first sediment level           
     410      jn = nvar 
     411      DO ji = 1, jpoce 
     412         IF (accmask(ji) == 0) THEN 
     413         aplus  = ( por1(2) + por1(3) ) / 2.0 
     414         dxplus = ( dz(2) + dz(3) ) / 2. 
     415         rplus  = ( 1.0 / vols(2) ) * db(ji,2) * aplus / dxplus * rads1sol(2,jn) 
     416 
     417         za(2) = 0. 
     418         zb(2) = rplus  
     419         zc(2) = -rplus 
     420 
     421         DO jk = 3, nksed - 1 
     422            aminus  = ( por1(jk-1) + por1(jk) ) * 0.5 
     423            aminus  = ( ( vols(jk-1) / dz(jk-1) ) + ( vols(jk) / dz(jk) ) ) / 2. 
     424            dxminus = ( dz(jk-1) + dz(jk) ) / 2. 
     425            rminus  = ( 1.0 / vols(jk) ) * db(ji,jk-1) * aminus / dxminus * rads1sol(jk,jn) 
     426            ! 
     427            aplus   = ( por1(jk) + por1(jk+1) ) * 0.5 
     428            dxplus  = ( dz(jk) + dz(jk+1) ) / 2. 
     429            rplus   = ( 1.0 / vols(jk) ) * db(ji,jk) * aplus / dxplus * rads1sol(jk,jn) 
     430            !      
     431            za(jk) = -rminus 
     432            zb(jk) = rminus + rplus 
     433            zc(jk) = -rplus 
     434 
     435         ENDDO 
     436 
     437         aminus = ( por1(nksed-1) + por1(nksed) ) * 0.5 
     438         dxminus = ( dz(nksed-1) + dz(nksed) ) / 2. 
     439         rminus  = ( 1.0 / vols(nksed) ) * db(ji,nksed-1) * aminus / dxminus * rads1sol(nksed,jn) 
     440         ! 
     441         za(nksed) = -rminus 
     442         zb(nksed) = rminus 
     443         zc(nksed) = 0. 
     444 
     445         ! solves tridiagonal system of linear equations  
     446         ! -----------------------------------------------     
     447         pwcpa(ji,2,jn) = pwcpa(ji,2,jn) - ( zc(2) * pwcp(ji,3,jn) + zb(2) * pwcp(ji,2,jn) ) 
     448         DO jk = 3, nksed-1 
     449            pwcpa(ji,jk,jn) =  pwcpa(ji,jk,jn) - ( zc(jk) * pwcp(ji,jk+1,jn) + za(jk) * pwcp(ji,jk-1,jn)    & 
     450            &                  + zb(jk) * pwcp(ji,jk,jn) ) 
     451         ENDDO 
     452         pwcpa(ji,nksed,jn) = pwcpa(ji,nksed,jn) - ( za(nksed) * pwcp(ji,nksed-1,jn)    & 
     453         &                     + zb(nksed) * pwcp(ji,nksed,jn) ) 
     454         !  
     455         ENDIF 
     456      END DO 
     457      ! 
     458      IF( ln_timing )  CALL timing_stop('sed_mat_btb') 
     459        
     460    END SUBROUTINE sed_mat_btb 
     461 
     462    SUBROUTINE sed_mat_btbjac( nksed, nvar, NEQ, NROWPD, jacvode, accmask ) 
     463       !!--------------------------------------------------------------------- 
     464       !!                  ***  ROUTINE sed_mat_btb  *** 
     465       !! 
     466       !! ** Purpose :  solves tridiagonal system of linear equations  
     467       !! 
     468       !! ** Method  :  
     469       !!        1 - computes left hand side of linear system of equations 
     470       !!            for dissolution reaction 
     471       !! 
     472       !!         2 - forward/backward substitution.  
     473       !! 
     474       !!   History : 
     475       !!        !  04-10 (N. Emprin, M. Gehlen ) original 
     476       !!        !  06-04 (C. Ethe)  Module Re-organization 
     477       !!---------------------------------------------------------------------- 
     478       !! * Arguments 
     479       INTEGER , INTENT(in) ::  nvar, nksed, NEQ, NROWPD  ! number of variable 
     480       REAL, DIMENSION(jpoce,NROWPD,NEQ), INTENT(inout) :: jacvode 
     481       INTEGER, DIMENSION(jpoce), INTENT(in) :: accmask 
     482 
     483       !---Local declarations 
     484       INTEGER  ::  ji, jk, jn, jnn, jni, jnj ,jnij 
     485 
     486       REAL(wp) ::  & 
     487          aplus,aminus   ,  & 
     488          rplus,rminus   ,  & 
     489          dxplus,dxminus 
     490 
     491       REAL(wp), DIMENSION(nksed)      :: za, zb, zc 
     492 
     493       !---------------------------------------------------------------------- 
     494 
     495      ! Computation left hand side of linear system of  
     496      ! equations for dissolution reaction 
     497      !--------------------------------------------- 
     498      IF( ln_timing )  CALL timing_start('sed_mat_btbjac') 
     499 
     500      ! first sediment level           
     501      jn = nvar 
     502      DO ji = 1, jpoce 
     503         IF (accmask(ji) == 0) THEN 
     504         aplus  = ( por1(2) + por1(3) ) / 2.0 
     505         dxplus = ( dz(2) + dz(3) ) / 2. 
     506         rplus  = ( 1.0 / vols(2) ) * db(ji,2) * aplus / dxplus * rads1sol(2,jn) 
     507 
     508         za(2) = 0. 
     509         zb(2) = rplus 
     510         zc(2) = -rplus 
     511 
     512         DO jk = 3, nksed - 1 
     513            aminus  = ( por1(jk-1) + por1(jk) ) * 0.5 
     514            aminus  = ( ( vols(jk-1) / dz(jk-1) ) + ( vols(jk) / dz(jk) ) ) / 2. 
     515            dxminus = ( dz(jk-1) + dz(jk) ) / 2. 
     516            rminus  = ( 1.0 / vols(jk) ) * db(ji,jk-1) * aminus / dxminus * rads1sol(jk,jn) 
     517            ! 
     518            aplus   = ( por1(jk) + por1(jk+1) ) * 0.5 
     519            dxplus  = ( dz(jk) + dz(jk+1) ) / 2. 
     520            rplus   = ( 1.0 / vols(jk) ) * db(ji,jk) * aplus / dxplus * rads1sol(jk,jn) 
     521            !      
     522            za(jk) = -rminus 
     523            zb(jk) = rminus + rplus 
     524            zc(jk) = -rplus 
     525 
     526         ENDDO 
     527 
     528         aminus = ( por1(nksed-1) + por1(nksed) ) * 0.5 
     529         dxminus = ( dz(nksed-1) + dz(nksed) ) / 2. 
     530         rminus  = ( 1.0 / vols(nksed) ) * db(ji,nksed-1) * aminus / dxminus * rads1sol(nksed,jn) 
     531         ! 
     532         za(nksed) = -rminus 
     533         zb(nksed) = rminus 
     534         zc(nksed) = 0. 
     535 
     536         ! solves tridiagonal system of linear equations  
     537         ! -----------------------------------------------     
     538         jnn = isvode(jn) 
     539         jni = jpvode + jnn 
     540         jnij = jpvode + 1 
     541         jacvode(ji, jnij, jni) = jacvode(ji,jnij,jni) - zb(2) 
     542         jnj = 2 * jpvode + jnn 
     543         jnij = jni - jnj + jpvode + 1 
     544         jacvode(ji, jnij, jnj) = jacvode(ji, jnij, jnj) -zc(2) 
     545         DO jk = 3, nksed-1 
     546            jni = (jk-1) * jpvode + jnn 
     547            jnj = (jk-2) * jpvode + jnn 
     548            jnij = jni - jnj + jpvode + 1 
     549            jacvode(ji, jnij, jnj) = jacvode(ji, jnij, jnj) - za(jk) 
     550            jnj = (jk-1) * jpvode + jnn 
     551            jnij = jni - jnj + jpvode + 1 
     552            jacvode(ji, jnij, jnj) = jacvode(ji, jnij, jnj) - zb(jk) 
     553            jnj = (jk) * jpvode + jnn 
     554            jnij = jni - jnj + jpvode + 1 
     555            jacvode(ji, jnij, jnj)   = jacvode(ji, jnij, jnj) - zc(jk) 
     556         ENDDO 
     557         jni = (nksed-1) * jpvode + jnn 
     558         jnj = (nksed-2) * jpvode + jnn 
     559         jnij = jni - jnj + jpvode + 1 
     560         jacvode(ji, jnij, jnj) = jacvode(ji, jnij, jnj) - za(nksed) 
     561         jnij = jpvode + 1 
     562         jacvode(ji, jnij, jni) = jacvode(ji, jnij, jni) - zb(nksed) 
     563         !  
     564         ENDIF 
     565      END DO 
     566      ! 
     567      IF( ln_timing )  CALL timing_stop('sed_mat_btbjac') 
     568 
     569    END SUBROUTINE sed_mat_btbjac 
     570 
     571 
     572    SUBROUTINE sed_mat_dsri( nksed, nvar, preac, psms, dtsed_in, psol ) 
     573       !!--------------------------------------------------------------------- 
     574       !!                  ***  ROUTINE sed_mat_dsr  *** 
     575       !! 
     576       !! ** Purpose :  solves tridiagonal system of linear equations  
     577       !! 
     578       !! ** Method  :  
     579       !!        1 - computes left hand side of linear system of equations 
     580       !!            for dissolution reaction 
     581       !!         For mass balance in kbot+sediment : 
     582       !!              dz3d  (:,1) = dz(1) = 0.5 cm 
     583       !!              volw3d(:,1) = dzkbot ( see sedini.F90 )  
     584       !!              dz(2)       = 0.3 cm  
     585       !!              dz3d(:,2)   = 0.3 + dzdep   ( see seddsr.F90 )      
     586       !!              volw3d(:,2) and vols3d(l,2) are thickened ( see 
     587       !seddsr.F90 )  
     588       !! 
     589       !!         2 - forward/backward substitution.  
     590       !! 
     591       !!   History : 
     592       !!        !  04-10 (N. Emprin, M. Gehlen ) original 
     593       !!        !  06-04 (C. Ethe)  Module Re-organization 
     594       !!---------------------------------------------------------------------- 
     595       !! * Arguments 
     596       INTEGER , INTENT(in) ::  nksed, nvar  ! number of variable 
     597 
     598       REAL(wp), DIMENSION(jpoce,nksed), INTENT(in   ) :: preac  ! reaction rates 
     599       REAL(wp), DIMENSION(jpoce,nksed), INTENT(in   ) :: psms  ! reaction rates 
     600       REAL(wp), DIMENSION(jpoce,nksed), INTENT(inout) :: psol  ! reaction rates 
     601       REAL(wp), INTENT(in) ::  dtsed_in 
     602 
     603       !---Local declarations 
     604       INTEGER  ::  ji, jk, jn 
     605       REAL(wp), DIMENSION(jpoce,nksed) :: za, zb, zc, zr 
     606       REAL(wp), DIMENSION(jpoce)        :: zbet 
     607       REAL(wp), DIMENSION(jpoce,nksed) :: zgamm 
     608 
     609       REAL(wp) ::  rplus,rminus 
     610       !---------------------------------------------------------------------- 
     611 
     612       IF( ln_timing )  CALL timing_start('sed_mat_dsri') 
     613 
     614       ! Computation left hand side of linear system of  
     615       ! equations for dissolution reaction 
     616       !--------------------------------------------- 
     617       jn = nvar 
     618       ! first sediment level           
     619       DO ji = 1, jpoce 
     620          rplus  = dtsed_in * apluss(ji,1) * diff(ji,1,jn) * radssol(1,jn) 
    86621 
    87622          za(ji,1) = 0. 
     
    89624          zc(ji,1) = -rplus 
    90625       ENDDO 
    91   
    92        DO jk = 2, nlev - 1 
    93           DO ji = 1, ndim 
    94              aminus  = ( ( volw3d(ji,jk-1) / ( dz3d(ji,jk-1) ) ) + & 
    95              &        ( volw3d(ji,jk  ) / ( dz3d(ji,jk  ) ) ) ) / 2. 
    96              dxminus = ( dz3d(ji,jk-1) + dz3d(ji,jk) ) / 2. 
    97  
    98              aplus   = ( ( volw3d(ji,jk  ) / ( dz3d(ji,jk  ) ) ) + & 
    99              &        ( volw3d(ji,jk+1) / ( dz3d(ji,jk+1) ) ) ) / 2. 
    100              dxplus  = ( dz3d(ji,jk) + dz3d(ji,jk+1) ) / 2 
    101                 ! 
    102              rminus  = ( dtsed_in / volw3d(ji,jk) ) * diff(ji,jk-1,jn) * aminus / dxminus 
    103              rplus   = ( dtsed_in / volw3d(ji,jk) ) * diff(ji,jk,jn)   * aplus / dxplus 
     626 
     627       DO jk = 2, nksed - 1 
     628          DO ji = 1, jpoce 
     629             rminus  = dtsed_in * aminuss(ji,jk) * diff(ji,jk-1,jn) * radssol(jk,jn) 
     630             rplus   = dtsed_in * apluss (ji,jk) * diff(ji,jk,jn) * radssol(jk,jn) 
    104631                !      
    105632             za(ji,jk) = -rminus 
    106              zb(ji,jk) = 1. + rminus + rplus  
     633             zb(ji,jk) = 1. + rminus + rplus 
    107634             zc(ji,jk) = -rplus 
    108635          END DO 
    109636       END DO 
    110637 
    111        DO ji = 1, ndim 
    112           aminus  = ( ( volw3d(ji,nlev-1) / dz3d(ji,nlev-1)  ) + & 
    113           &        ( volw3d(ji,nlev)  / dz3d(ji,nlev) ) ) / 2. 
    114           dxminus = ( dz3d(ji,nlev-1) + dz3d(ji,nlev) ) / 2. 
    115           rminus  = ( dtsed_in / volw3d(ji,nlev) ) * diff(ji,nlev-1,jn) * aminus / dxminus 
     638       DO ji = 1, jpoce 
     639          rminus  = dtsed_in * aminuss(ji,nksed) * diff(ji,nksed-1,jn) * radssol(nksed,jn) 
    116640          ! 
    117           za(ji,nlev) = -rminus 
    118           zb(ji,nlev) = 1. + rminus 
    119           zc(ji,nlev) = 0. 
     641          za(ji,nksed) = -rminus 
     642          zb(ji,nksed) = 1. + rminus 
     643          zc(ji,nksed) = 0. 
    120644       END DO 
    121645 
     
    124648       ! ----------------------------------------------- 
    125649 
    126        zr  (:,:) = psol(:,:) + psms(:,:) 
    127        zb  (:,:) = zb(:,:) + preac(:,:) 
     650       zr  (:,:) = psol(:,:) + psms(:,:) * dtsed_in 
     651       zb  (:,:) = zb(:,:) - preac(:,:) * dtsed_in 
    128652       zbet(:  ) = zb(:,1) 
    129653       psol(:,1) = zr(:,1) / zbet(:) 
    130654 
    131655          !  
    132        DO jk = 2, nlev 
    133           DO ji = 1, ndim 
     656       DO jk = 2, nksed 
     657          DO ji = 1, jpoce 
    134658             zgamm(ji,jk) =  zc(ji,jk-1) / zbet(ji) 
    135659             zbet(ji)     =  zb(ji,jk) - za(ji,jk) * zgamm(ji,jk) 
     
    138662       ENDDO 
    139663          !  
    140        DO jk = nlev - 1, 1, -1 
    141           DO ji = 1,ndim 
     664       DO jk = nksed - 1, 1, -1 
     665          DO ji = 1, jpoce 
    142666             psol(ji,jk) = psol(ji,jk) - zgamm(ji,jk+1) * psol(ji,jk+1) 
    143667          END DO 
    144668       ENDDO 
    145669 
    146        IF( ln_timing )  CALL timing_stop('sed_mat_dsr') 
    147  
    148  
    149     END SUBROUTINE sed_mat_dsr 
    150  
    151     SUBROUTINE sed_mat_btb( nvar, ndim, nlev, psol, dtsed_in ) 
    152        !!--------------------------------------------------------------------- 
    153        !!                  ***  ROUTINE sed_mat_btb  *** 
    154        !! 
    155        !! ** Purpose :  solves tridiagonal system of linear equations  
    156        !! 
    157        !! ** Method  :  
    158        !!        1 - computes left hand side of linear system of equations 
    159        !!            for dissolution reaction 
    160        !! 
    161        !!         2 - forward/backward substitution.  
    162        !! 
    163        !!   History : 
    164        !!        !  04-10 (N. Emprin, M. Gehlen ) original 
    165        !!        !  06-04 (C. Ethe)  Module Re-organization 
    166        !!---------------------------------------------------------------------- 
    167        !! * Arguments 
    168        INTEGER , INTENT(in) :: & 
    169           nvar , &  ! number of variables 
    170           ndim , &  ! number of points 
    171           nlev      ! number of sediment levels 
    172  
    173       REAL(wp), DIMENSION(ndim,nlev,nvar), INTENT(inout) :: & 
    174           psol      ! solution 
    175  
    176       REAL(wp), INTENT(in) :: dtsed_in 
    177  
    178        !---Local declarations 
    179        INTEGER  ::  & 
    180           ji, jk, jn 
    181  
    182        REAL(wp) ::  & 
    183           aplus,aminus   ,  &  
    184           rplus,rminus   ,  & 
    185           dxplus,dxminus  
    186  
    187        REAL(wp), DIMENSION(nlev)      :: za, zb, zc 
    188        REAL(wp), DIMENSION(ndim,nlev) :: zr 
    189        REAL(wp), DIMENSION(nmax)      :: zgamm  
    190        REAL(wp) ::  zbet 
    191  
    192   
    193        !---------------------------------------------------------------------- 
    194  
    195        ! Computation left hand side of linear system of  
    196        ! equations for dissolution reaction 
    197        !--------------------------------------------- 
    198  
    199  
    200       IF( ln_timing )  CALL timing_start('sed_mat_btb') 
    201  
    202        ! first sediment level           
    203       DO ji = 1, ndim 
    204          aplus  = ( ( vols(2) / dz(2) ) + ( vols(3) / dz(3) ) ) / 2. 
    205          dxplus = ( dz(2) + dz(3) ) / 2. 
    206          rplus  = ( dtsed_in / vols(2) ) * db(ji,2) * aplus / dxplus  
    207  
    208          za(1) = 0. 
    209          zb(1) = 1. + rplus 
    210          zc(1) = -rplus 
    211  
    212               
    213          DO jk = 2, nlev - 1 
    214             aminus  = ( ( vols(jk) / dz(jk) ) + ( vols(jk+1) / dz(jk+1) ) ) / 2. 
    215             dxminus = ( dz(jk) + dz(jk+1) ) / 2. 
    216             rminus  = ( dtsed_in / vols(jk+1) ) * db(ji,jk) * aminus / dxminus 
    217             ! 
    218             aplus   = ( ( vols(jk+1) / dz(jk+1  ) ) + ( vols(jk+2) / dz(jk+2) ) ) / 2. 
    219             dxplus  = ( dz(jk+1) + dz(jk+2) ) / 2. 
    220             rplus   = ( dtsed_in / vols(jk+1) ) * db(ji,jk+1) * aplus / dxplus 
    221             !      
    222             za(jk) = -rminus 
    223             zb(jk) = 1. + rminus + rplus  
    224             zc(jk) = -rplus 
    225          ENDDO 
    226   
    227          aminus  = ( ( vols(nlev) / dz(nlev) ) + ( vols(nlev+1) / dz(nlev+1) ) ) / 2. 
    228          dxminus = ( dz(nlev) + dz(nlev+1) ) / 2. 
    229          rminus  = ( dtsed_in / vols(nlev+1) ) * db(ji,nlev) * aminus / dxminus 
    230          ! 
    231          za(nlev) = -rminus 
    232          zb(nlev) = 1. + rminus 
    233          zc(nlev) = 0. 
    234  
    235  
    236          ! solves tridiagonal system of linear equations  
    237          ! -----------------------------------------------     
    238          DO jn = 1, nvar 
    239            
    240             DO jk = 1, nlev 
    241                zr  (ji,jk)    = psol(ji,jk,jn) 
    242             END DO 
    243             zbet          = zb(1) 
    244             psol(ji,1,jn) = zr(ji,1) / zbet 
    245             !  
    246             DO jk = 2, nlev 
    247                zgamm(jk) =  zc(jk-1) / zbet 
    248                zbet      =  zb(jk) - za(jk) * zgamm(jk) 
    249                psol(ji,jk,jn) = ( zr(ji,jk) - za(jk) * psol(ji,jk-1,jn) ) / zbet 
    250             ENDDO 
    251             !  
    252             DO jk = nlev - 1, 1, -1 
    253                 psol(ji,jk,jn) = psol(ji,jk,jn) - zgamm(jk+1) * psol(ji,jk+1,jn) 
    254             ENDDO 
    255  
    256          ENDDO 
    257  
    258        END DO 
    259        ! 
    260        IF( ln_timing )  CALL timing_stop('sed_mat_btb') 
    261  
    262         
    263     END SUBROUTINE sed_mat_btb 
     670       IF( ln_timing )  CALL timing_stop('sed_mat_dsri') 
     671 
     672 
     673    END SUBROUTINE sed_mat_dsri 
     674 
    264675 
    265676 END MODULE sedmat 
  • NEMO/trunk/src/TOP/PISCES/SED/sedorg.F90

    r10225 r15450  
    11MODULE sedorg 
    22   !!====================================================================== 
    3    !!              ***  MODULE  seddsr  *** 
    4    !!    Sediment : dissolution and reaction in pore water related  
    5    !!    related to organic matter 
     3   !!              ***  MODULE  sedinorg  *** 
     4   !!    Sediment : dissolution and reaction in pore water of  
     5   !!               inorganic species 
    66   !!===================================================================== 
    77   !! * Modules used 
     
    99   USE sed_oce 
    1010   USE sedini 
    11    USE seddiff 
    12    USE seddsr 
     11   USE sedmat 
    1312   USE lib_mpp         ! distribued memory computing library 
    1413   USE lib_fortran 
     
    1918   PUBLIC sed_org 
    2019 
    21    !! * Module variables 
    22  
    23    REAL(wp) :: zadsnh4 
    24  
    2520   !! $Id: seddsr.F90 5215 2015-04-15 16:11:56Z nicolasmartin $ 
    2621CONTAINS 
    2722    
    28    SUBROUTINE sed_org( kt )  
     23   SUBROUTINE sed_org( kt ) 
    2924      !!---------------------------------------------------------------------- 
    3025      !!                   ***  ROUTINE sed_org  *** 
    3126      !!  
    32       !!  ** Purpose :  computes pore water diffusion and reaction 
     27      !!  ** Purpose :  computes pore water dissolution and reaction 
    3328      !! 
    34       !!  ** Methode :  Computation of the redox reactions in sediment. 
    35       !!                The main redox reactions are solved in sed_dsr whereas 
    36       !!                the secondary reactions are solved in sed_dsr_redoxb. 
    37       !!                A strand spliting approach is being used here (see  
    38       !!                sed_dsr_redoxb for more information).  
    39       !!                Diffusive fluxes are computed in sed_diff 
     29      !!  ** Methode :  implicit simultaneous computation of undersaturation 
     30      !!               resulting from diffusive pore water transport and chemical 
     31      !!               pore water reactions. Solid material is consumed according 
     32      !!               to redissolution and remineralisation 
     33      !! 
     34      !!  ** Remarks : 
     35      !!              - undersaturation : deviation from saturation concentration 
     36      !!              - reaction rate   : sink of undersaturation from dissolution 
     37      !!                                 of solid material  
    4038      !! 
    4139      !!   History : 
     
    4341      !!        !  04-10 (N. Emprin, M. Gehlen ) f90 
    4442      !!        !  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. 
     43      !!        !  19-08 (O. Aumont) Debugging and improvement of the model 
    4944      !!---------------------------------------------------------------------- 
    5045      !! Arguments 
    51       INTEGER, INTENT(in) ::   kt 
     46      INTEGER, INTENT(in)  :: kt   ! time step 
    5247      ! --- local variables 
    53       INTEGER  :: ji, jk, js, jw, jnt   ! dummy looop indices 
    54       REAL(wp) :: zadsnh4 
     48      REAL(wp), DIMENSION(jpoce, jpksed) :: psms, preac 
    5549      !! 
    5650      !!---------------------------------------------------------------------- 
    5751 
    5852      IF( ln_timing )  CALL timing_start('sed_org') 
     53 
     54      IF( kt == nitsed000 ) THEN 
     55         IF (lwp) WRITE(numsed,*) ' sed_org : solute species which do not experience redox reactions ' 
     56         IF (lwp) WRITE(numsed,*) ' ' 
     57      ENDIF 
    5958! 
    60       IF( kt == nitsed000 ) THEN 
    61          IF (lwp) THEN 
    62             WRITE(numsed,*) ' sed_org : Organic degradation related reactions and diffusion' 
    63             WRITE(numsed,*) ' ' 
    64          ENDIF 
    65 !         !  
    66          dens_mol_wgt(1:jpsol) = denssol / mol_wgt(1:jpsol) 
    67          !  
    68       ENDIF 
     59      ! DIC in pore water 
     60      preac(:,:) = 0.0_wp 
     61      psms (:,:) = rearatpom(:,:) 
     62      CALL sed_mat_dsri( jpksed, jwdic, preac, psms, dtsed, pwcp(:,:,jwdic) ) 
     63       
     64      ! Silicate in pore water 
     65      psms (:,:) = 0.0_wp 
     66      CALL sed_mat_dsri( jpksed, jwsil, preac, psms, dtsed, pwcp(:,:,jwsil) ) 
    6967 
    70       dtsed2 = dtsed / REAL( nrseddt, wp ) 
    71  
    72       ! 1. Change of geometry 
    73       !    Increase of dz3d(2) thickness : dz3d(2) = dz3d(2)+dzdep 
    74       !    Warning : no change for dz(2) 
    75       !--------------------------------------------------------- 
    76       dz3d(1:jpoce,2) = dz3d(1:jpoce,2) + dzdep(1:jpoce) 
    77  
    78       ! New values for volw3d(:,2) and vols3d(:,2) 
    79       ! Warning : no change neither for volw(2) nor  vols(2) 
    80       !------------------------------------------------------ 
    81       volw3d(1:jpoce,2) = dz3d(1:jpoce,2) * por(2) 
    82       vols3d(1:jpoce,2) = dz3d(1:jpoce,2) * por1(2) 
    83  
    84       ! 2. Change of previous solid fractions (due to volum changes) for k=2 
    85       !--------------------------------------------------------------------- 
    86  
    87       DO js = 1, jpsol 
    88          DO ji = 1, jpoce 
    89             solcp(ji,2,js) = solcp(ji,2,js) * dz(2) / dz3d(ji,2) 
    90          ENDDO 
    91       END DO 
    92  
    93       ! 3. New solid fractions (including solid rain fractions) for k=2 
    94       !------------------------------------------------------------------ 
    95       DO js = 1, jpsol 
    96          DO ji = 1, jpoce 
    97             IF (raintg(ji) .ne. 0) THEN 
    98                solcp(ji,2,js) = solcp(ji,2,js) + & 
    99                &           ( rainrg(ji,js) / raintg(ji) ) * ( dzdep(ji) / dz3d(ji,2) ) 
    100                ! rainrm are temporary cancel 
    101                rainrm(ji,js) = 0. 
    102             ENDIF 
    103          END DO 
    104       ENDDO 
    105  
    106       ! 4.  Adjustment of bottom water concen.(pwcp(1)): 
    107       !     We impose that pwcp(2) is constant. Including dzdep in dz3d(:,2) we assume 
    108       !     that dzdep has got a porosity of por(2). So pore water volum of jk=2 increase. 
    109       !     To keep pwcp(2) cste we must compensate this "increase" by a slight adjusment 
    110       !     of bottom water concentration. 
    111       !     This adjustment is compensate at the end of routine 
    112       !------------------------------------------------------------- 
    113       DO jw = 1, jpwat 
    114          DO ji = 1, jpoce 
    115             pwcp(ji,1,jw) = pwcp(ji,1,jw) - & 
    116                &            pwcp(ji,2,jw) * dzdep(ji) * por(2) / ( dzkbot(ji) + rtrn ) 
    117          END DO 
    118       ENDDO 
    119  
    120       zadsnh4 = 1.0 / ( 1.0 + adsnh4 ) 
    121  
    122       ! -------------------------------------------------- 
    123       ! Computation of the diffusivities 
    124       ! -------------------------------------------------- 
    125  
    126       DO js = 1, jpwat 
    127          DO jk = 1, jpksed 
    128             DO ji = 1, jpoce 
    129                diff(ji,jk,js) = ( diff1(js) + diff2(js) * temp(ji) ) / ( 1.0 - 2.0 * log( por(jk) ) ) 
    130             END DO 
    131          END DO 
    132       END DO 
    133  
    134       ! Impact of bioirrigation and adsorption on diffusion 
    135       ! --------------------------------------------------- 
    136  
    137       diff(:,:,jwnh4) = diff(:,:,jwnh4) * ( 1.0 + irrig(:,:) ) * zadsnh4 
    138       diff(:,:,jwsil) = diff(:,:,jwsil) * ( 1.0 + irrig(:,:) ) 
    139       diff(:,:,jwoxy) = diff(:,:,jwoxy) * ( 1.0 + irrig(:,:) ) 
    140       diff(:,:,jwdic) = diff(:,:,jwdic) * ( 1.0 + irrig(:,:) ) 
    141       diff(:,:,jwno3) = diff(:,:,jwno3) * ( 1.0 + irrig(:,:) ) 
    142       diff(:,:,jwpo4) = diff(:,:,jwpo4) * ( 1.0 + irrig(:,:) ) 
    143       diff(:,:,jwalk) = diff(:,:,jwalk) * ( 1.0 + irrig(:,:) ) 
    144       diff(:,:,jwh2s) = diff(:,:,jwh2s) * ( 1.0 + irrig(:,:) ) 
    145       diff(:,:,jwso4) = diff(:,:,jwso4) * ( 1.0 + irrig(:,:) ) 
    146       diff(:,:,jwfe2) = diff(:,:,jwfe2) * ( 1.0 + 0.2 * irrig(:,:) ) 
    147  
    148       DO jnt = 1, nrseddt 
    149          CALL sed_diff( kt, jnt )        ! 1st pass in diffusion to get values at t+1/2 
    150          CALL sed_dsr ( kt, jnt )        ! Dissolution reaction 
    151          CALL sed_diff( kt, jnt )        ! 2nd pass in diffusion to get values at t+1 
    152       END DO 
    153  
     68      ! Iron ligands in pore water 
     69      psms (:,:) = ratligc * rearatpom(:,:) 
     70      preac(:,:) = -reac_ligc 
     71      CALL sed_mat_dsri( jpksed, jwlgw, preac, psms, dtsed, pwcp(:,:,jwlgw) ) 
    15472 
    15573      IF( ln_timing )  CALL timing_stop('sed_org') 
  • NEMO/trunk/src/TOP/PISCES/SED/sedrst.F90

    r14239 r15450  
    4242      CHARACTER(LEN=50)   ::   clname   ! trc output restart file name 
    4343      CHARACTER(LEN=256)  ::   clpath   ! full path to ocean output restart file 
    44       CHARACTER(LEN=52)   ::   clpname   ! trc output restart file name including AGRIF 
     44      CHARACTER(LEN=3)    ::   cdcomp 
    4545      !!---------------------------------------------------------------------- 
    4646      ! 
     
    6565 
    6666      IF( .NOT. ln_rst_list .AND. nn_stock == -1 )   RETURN   ! we will never do any restart 
    67  
    6867      ! to get better performances with NetCDF format: 
    69       ! we open and define the tracer restart file one tracer time step before writing the data (-> at nitrst - 1) 
    70       ! except if we write tracer restart files every tracer time step or if a tracer restart file was writen at nitend - 1 
    71       IF( kt == nitrst - 2*nn_dtsed .OR. nn_stock == nn_dtsed .OR. ( kt == nitend - nn_dtsed .AND. .NOT. lrst_sed ) ) THEN 
     68      ! we open and define the tracer restart file one tracer time step before writing the data (-> at nitrst - 2*nn_dttrc + 1) 
     69      ! except if we write tracer restart files every tracer time step or if a tracer restart file was writen at nitend - 2*nn_dttrc + 1 
     70      IF( kt == nitrst - 1 .OR. nn_stock == 1 .OR. ( kt == nitend - 1 .AND. .NOT. lrst_sed ) ) THEN 
    7271         ! beware of the format used to write kt (default is i8.8, that should be large enough) 
    7372         IF( nitrst > 1.0e9 ) THEN   ;   WRITE(clkt,*       ) nitrst 
     
    8180         IF(lwp) WRITE(numsed,*) & 
    8281             '             open sed restart.output NetCDF file: ',TRIM(clpath)//clname 
    83          IF(.NOT.lwxios) THEN 
    84             CALL iom_open( TRIM(clpath)//TRIM(clname), numrsw, ldwrt = .TRUE., kdlev = jpksed, cdcomp = 'SED' ) 
    85          ELSE 
    86 #if defined key_xios 
    87             cw_sedrst_cxt = "rstws_"//TRIM(ADJUSTL(clkt)) 
    88             IF( TRIM(Agrif_CFixed()) == '0' ) THEN 
    89                clpname = clname 
    90             ELSE 
    91                clpname = TRIM(Agrif_CFixed())//"_"//clname 
    92             ENDIF 
    93             numrsw = iom_xios_setid(TRIM(clpath)//TRIM(clpname)) 
    94             CALL iom_init( cw_sedrst_cxt, kdid = numrsw, ld_closedef = .FALSE. ) 
    95 #else 
    96             CALL ctl_stop( 'Can not use XIOS in trc_rst_opn' ) 
    97 #endif 
    98          ENDIF 
    99  
     82         cdcomp ='SED' 
     83         CALL iom_open( TRIM(clpath)//TRIM(clname), numrsw, ldwrt = .TRUE., kdlev = jpksed, cdcomp = cdcomp ) 
    10084         lrst_sed = .TRUE. 
    10185      ENDIF 
     
    204188         &             zdta2(1:jpi,1:jpj,1:jpksed), iarroce(1:jpoce) ) 
    205189 
    206       cltra = "sedligand" 
    207       IF( iom_varid( numrsr, TRIM(cltra) , ldstop = .FALSE. ) > 0 ) THEN 
    208          CALL iom_get( numrsr, jpdom_auto, TRIM(cltra), zdta2(:,:,:) ) 
    209       ELSE 
    210          zdta2(:,:,:) = 0.0 
    211       ENDIF 
    212  
    213       CALL pack_arr( jpoce, sedligand(1:jpoce,1:jpksed), & 
    214          &             zdta2(1:jpi,1:jpj,1:jpksed), iarroce(1:jpoce) ) 
    215190      IF( ln_timing )  CALL timing_stop('sed_rst_read') 
    216191      
     
    247222      IF(lwp) WRITE(numsed,*) '~~~~~~~~~' 
    248223 
    249  
    250       trcsedi(:,:,:,:)   = 0.0 
    251       flxsedi3d(:,:,:,:) = 0.0 
    252224      zdta(:,:)          = 1.0 
    253225      zdta2(:,:,:)       = 0.0 
    254  
    255226          
    256227      !! 1. WRITE in nutwrs 
    257228      !! ------------------ 
    258 !     zinfo(1) = REAL( kt) 
    259       CALL iom_rstput( kt, nitrst, numrsw, 'kt', REAL( kt    , wp) ) 
     229      zinfo(1) = REAL( kt) 
     230      CALL iom_rstput( kt, nitrst, numrsw, 'kt', zinfo ) 
    260231 
    261232      ! Back to 2D geometry 
    262233      DO jn = 1, jpsol 
    263          CALL unpack_arr( jpoce, trcsedi(1:jpi,1:jpj,1:jpksed,jn) , iarroce(1:jpoce), & 
     234         CALL unpack_arr( jpoce, zdta2(1:jpi,1:jpj,1:jpksed) , iarroce(1:jpoce), & 
    264235         &                       solcp(1:jpoce,1:jpksed,jn ) ) 
     236         cltra = TRIM(sedtrcd(jn)) 
     237         CALL iom_rstput( kt, nitrst, numrsw, TRIM(cltra), zdta2(:,:,:) ) 
    265238      END DO 
    266239 
    267240      DO jn = 1, jpwat 
    268          CALL unpack_arr( jpoce, trcsedi(1:jpi,1:jpj,1:jpksed,jpsol+jn) , iarroce(1:jpoce), & 
     241         CALL unpack_arr( jpoce, zdta2(1:jpi,1:jpj,1:jpksed) , iarroce(1:jpoce), & 
    269242         &                       pwcp(1:jpoce,1:jpksed,jn  )  ) 
     243         cltra = TRIM(sedtrcd(jpsol+jn)) 
     244         CALL iom_rstput( kt, nitrst, numrsw, TRIM(cltra), zdta2(:,:,:) ) 
    270245      END DO 
    271246      ! pH 
     
    276251      ENDDO 
    277252 
    278       CALL unpack_arr( jpoce, flxsedi3d(1:jpi,1:jpj,1:jpksed,1)  , iarroce(1:jpoce), & 
     253      CALL unpack_arr( jpoce, zdta2(1:jpi,1:jpj,1:jpksed)  , iarroce(1:jpoce), & 
    279254      &                   zdta(1:jpoce,1:jpksed)  ) 
     255      cltra = TRIM(seddia3d(1)) 
     256      CALL iom_rstput( kt, nitrst, numrsw, TRIM(cltra), zdta2(:,:,:) ) 
    280257          
    281       CALL unpack_arr( jpoce, flxsedi3d(1:jpi,1:jpj,1:jpksed,2)  , iarroce(1:jpoce), & 
     258      CALL unpack_arr( jpoce, zdta2(1:jpi,1:jpj,1:jpksed)  , iarroce(1:jpoce), & 
    282259      &                   co3por(1:jpoce,1:jpksed)  ) 
     260      cltra = TRIM(seddia3d(2)) 
     261      CALL iom_rstput( kt, nitrst, numrsw, TRIM(cltra), zdta2(:,:,:) ) 
    283262 
    284263      ! prognostic variables 
    285264      ! -------------------- 
    286265 
    287       DO jn = 1, jptrased 
    288          cltra = TRIM(sedtrcd(jn)) 
    289          CALL iom_rstput( kt, nitrst, numrsw, TRIM(cltra), trcsedi(:,:,:,jn) ) 
    290       ENDDO 
    291  
    292       DO jn = 1, 2 
    293          cltra = TRIM(seddia3d(jn)) 
    294          CALL iom_rstput( kt, nitrst, numrsw, TRIM(cltra), flxsedi3d(:,:,:,jn) ) 
    295       ENDDO 
    296  
    297266      CALL unpack_arr( jpoce, zdta2(1:jpi,1:jpj,1:jpksed)  , iarroce(1:jpoce), & 
    298267      &                   db(1:jpoce,1:jpksed)  ) 
     
    307276      CALL iom_rstput( kt, nitrst, numrsw, TRIM(cltra), zdta2(:,:,:) ) 
    308277 
    309       CALL unpack_arr( jpoce, zdta2(1:jpi,1:jpj,1:jpksed)  , iarroce(1:jpoce), & 
    310       &                   sedligand(1:jpoce,1:jpksed)  ) 
    311  
    312       cltra = "sedligand" 
    313       CALL iom_rstput( kt, nitrst, numrsw, TRIM(cltra), zdta2(:,:,:) ) 
    314  
    315278      IF( kt == nitrst ) THEN 
    316           IF(.NOT.lwxios) THEN 
    317              CALL iom_close( numrsw )     ! close the restart file (only at last time step) 
    318           ELSE 
    319              CALL iom_context_finalize( cw_sedrst_cxt )  
    320              iom_file(numrsw)%nfid       = 0 
    321              numrsw = 0 
    322           ENDIF 
     279          CALL iom_close( numrsw )     ! close the restart file (only at last time step) 
    323280          IF( l_offline .AND. ln_rst_list ) THEN 
    324281             nrst_lst = nrst_lst + 1 
     
    351308      !!       In both those options, the  exact duration of the experiment 
    352309      !!       since the beginning (cumulated duration of all previous restart runs) 
    353       !!       is not stored in the restart and is assumed to be (nittrc000-1)*rn_Dt. 
     310      !!       is not stored in the restart and is assumed to be (nittrc000-1)*rdt. 
    354311      !!       This is valid is the time step has remained constant. 
    355312      !! 
     
    363320      REAL(wp) ::  zkt, zrdttrc1 
    364321      REAL(wp) ::  zndastp 
    365       CHARACTER(len = 82) :: clpname 
    366322 
    367323      ! Time domain : restart 
     
    375331 
    376332         IF( ln_rst_sed ) THEN 
    377             lxios_sini = .FALSE. 
    378333            CALL iom_open( TRIM(cn_sedrst_indir)//'/'//cn_sedrst_in, numrsr ) 
    379  
    380             IF( lrxios) THEN 
    381                 cr_sedrst_cxt = 'sed_rst' 
    382                 IF(lwp) WRITE(numout,*) 'Enable restart reading by XIOS for SED' 
    383 !               IF( TRIM(Agrif_CFixed()) == '0' ) THEN 
    384 !                  clpname = cn_sedrst_in 
    385 !               ELSE 
    386 !                  clpname = TRIM(Agrif_CFixed())//"_"//cn_sedrst_in    
    387 !               ENDIF 
    388                 CALL iom_init( cr_sedrst_cxt, kdid = numrsr, ld_closedef = .TRUE. ) 
    389             ENDIF 
    390334            CALL iom_get ( numrsr, 'kt', zkt )   ! last time-step of previous run 
     335 
    391336            IF(lwp) THEN 
    392337               WRITE(numsed,*) ' *** Info read in restart : ' 
     
    401346            ENDIF 
    402347            ! Control of date  
    403             IF( nittrc000  - NINT( zkt ) /= nn_dtsed .AND.  nn_rstsed /= 0 )                                  & 
     348            IF( nittrc000  - NINT( zkt ) /= 1 .AND.  nn_rstsed /= 0 )                                  & 
    404349               &   CALL ctl_stop( ' ===>>>> : problem with nittrc000 for the restart',                 & 
    405350               &                  ' verify the restart file or rerun with nn_rsttr = 0 (namelist)' ) 
     
    414359             ELSE 
    415360               ndastp = ndate0 - 1     ! ndate0 read in the namelist in dom_nam 
    416                adatrj = ( REAL( nittrc000-1, wp ) * rn_Dt ) / rday 
     361               adatrj = ( REAL( nittrc000-1, wp ) * rdt ) / rday 
    417362               ! note this is wrong if time step has changed during run 
    418363            ENDIF 
     
    435380            IF(lwp) WRITE(numsed,*) 'trc_wri : write the TOP restart file (NetCDF) at it= ', kt, ' date= ', ndastp 
    436381            IF(lwp) WRITE(numsed,*) '~~~~~~~' 
    437             IF( lwxios ) CALL iom_init_closedef(cw_sedrst_cxt) 
    438382         ENDIF 
    439383         CALL iom_rstput( kt, nitrst, numrsw, 'kt'     , REAL( kt    , wp) )   ! time-step 
    440384         CALL iom_rstput( kt, nitrst, numrsw, 'ndastp' , REAL( ndastp, wp) )   ! date 
    441          CALL iom_rstput( kt, nitrst, numrsw, 'adatrj' , adatrj )   ! number of elapsed days since 
    442          !                                                                                      ! the begining of the run [s] 
     385         CALL iom_rstput( kt, nitrst, numrsw, 'adatrj' , adatrj            )   ! number of elapsed days since 
     386         !                                                                     ! the begining of the run [s] 
    443387      ENDIF 
    444388 
  • NEMO/trunk/src/TOP/PISCES/SED/sedsfc.F90

    r13295 r15450  
    1313   !! * Substitutions 
    1414#  include "do_loop_substitute.h90" 
    15    !! $Id$ 
     15#  include "domzgr_substitute.h90" 
     16 
    1617CONTAINS 
    1718 
     
    2829      !!* Arguments 
    2930      INTEGER, INTENT(in) ::  kt              ! time step 
    30       INTEGER, INTENT(in) ::  Kbb             ! time index 
     31      INTEGER, INTENT(in) ::  Kbb             ! time level indices 
    3132 
    3233      ! * local variables 
     
    4647      CALL unpack_arr ( jpoce, trc_data(1:jpi,1:jpj,7), iarroce(1:jpoce), pwcp(1:jpoce,1,jwnh4) ) 
    4748      CALL unpack_arr ( jpoce, trc_data(1:jpi,1:jpj,8), iarroce(1:jpoce), pwcp(1:jpoce,1,jwfe2) ) 
     49      CALL unpack_arr ( jpoce, trc_data(1:jpi,1:jpj,9), iarroce(1:jpoce), pwcp(1:jpoce,1,jwlgw) ) 
    4850 
    49  
    50       DO_2D( 1, 1, 1, 1 ) 
     51      DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) 
    5152         ikt = mbkt(ji,jj) 
    5253         IF ( tmask(ji,jj,ikt) == 1 ) THEN 
  • NEMO/trunk/src/TOP/PISCES/SED/sedstp.F90

    r13970 r15450  
    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 
    12    USE sedbtb   ! bioturbation 
     10   USE sedsol   ! Organic reactions and diffusion 
    1311   USE sedadv   ! vertical advection 
    14    USE sedmbc   ! mass balance calculation 
    1512   USE sedsfc   ! sediment surface data 
    1613   USE sedrst   ! restart 
    1714   USE sedwri   ! outputs 
     15   USE sedini 
    1816   USE trcdmp_sed 
    1917   USE lib_mpp         ! distribued memory computing library 
     
    2523   !! * Routine accessibility 
    2624   PUBLIC sed_stp  ! called by step.F90 
     25 
     26   !! * Substitutions 
     27#  include "do_loop_substitute.h90" 
     28#  include "domzgr_substitute.h90" 
    2729 
    2830   !! $Id$ 
     
    4446      !!        !  06-04 (C. Ethe)  Re-organization 
    4547      !!---------------------------------------------------------------------- 
    46       INTEGER, INTENT(in) ::   kt                ! number of iteration 
    47       INTEGER, INTENT(in) ::   Kbb, Kmm, Krhs    ! time level indices 
    48       INTEGER :: ji,jk,js,jn,jw 
     48      INTEGER, INTENT(in) ::   kt       ! number of iteration 
     49      INTEGER, INTENT(in) ::   Kbb, Kmm, Krhs  ! time level indices 
     50 
     51      INTEGER :: ji,jk,js,jn,jw,jkmax,jsmax 
    4952      !!---------------------------------------------------------------------- 
    50       IF( ln_timing )      CALL timing_start('sed_stp') 
     53      IF( ln_timing )           CALL timing_start('sed_stp') 
    5154        ! 
    5255                                CALL sed_rst_opn  ( kt )       ! Open tracer restart file  
     
    5659 
    5760      dtsed  = rDt_trc 
    58 !      dtsed2 = dtsed 
    5961      IF (kt /= nitsed000) THEN 
    60          CALL sed_dta( kt, Kbb, Kmm )       ! Load  Data for bot. wat. Chem and fluxes 
     62         CALL sed_dta( kt, Kbb, Kmm )    ! Load  Data for bot. wat. Chem and fluxes 
    6163      ENDIF 
    6264 
     
    6668         ENDIF 
    6769 
    68          CALL sed_btb( kt )         ! 1st pass of bioturbation at t+1/2 
    69          CALL sed_org( kt )         ! Organic related reactions and diffusion 
    70          CALL sed_inorg( kt )       ! Dissolution reaction 
    71          CALL sed_btb( kt )         ! 2nd pass of bioturbation at t+1 
    72          tokbot(:,:) = 0.0 
    73          DO jw = 1, jpwat 
    74             DO ji = 1, jpoce 
    75                tokbot(ji,jw) = pwcp(ji,1,jw) * 1.e-3 * dzkbot(ji) 
    76             END DO 
    77          ENDDO 
     70         CALL sed_sol( kt )        ! Solute diffusion and reactions  
    7871         CALL sed_adv( kt )         ! advection 
    7972         CALL sed_co3( kt )         ! pH actualization for saving 
    80          ! This routine is commented out since it does not work at all 
    81          CALL sed_mbc( kt )         ! cumulation for mass balance calculation 
    8273 
    83          IF (ln_sed_2way) CALL sed_sfc( kt, Kbb )         ! Give back new bottom wat chem to tracer model 
     74         IF (ln_sed_2way) CALL sed_sfc( kt, Kbb )   ! Give back new bottom wat chem to tracer model 
    8475      ENDIF 
    8576      CALL sed_wri( kt )         ! outputs 
    8677      IF( kt == nitsed000 ) THEN 
    8778          CALL iom_close( numrsr )       ! close input tracer restart file 
    88           IF(lrxios) CALL iom_context_finalize(      cr_sedrst_cxt  ) 
    89 !         IF(lwm) CALL FLUSH( numont )   ! flush namelist output 
     79!          IF(lwm) CALL FLUSH( numont )   ! flush namelist output 
    9080      ENDIF 
    9181      IF( lrst_sed )            CALL sed_rst_wri( kt )   ! restart file output 
    9282 
    93       IF( kt == nitsedend )  CLOSE( numsed ) 
     83      IF( kt == nitsedend )     CLOSE( numsed ) 
    9484 
    95       IF( ln_timing )   CALL timing_stop('sed_stp') 
     85      IF( ln_timing )           CALL timing_stop('sed_stp') 
    9686 
    9787   END SUBROUTINE sed_stp 
  • NEMO/trunk/src/TOP/PISCES/SED/sedwri.F90

    r14086 r15450  
    44   !!         Sediment diagnostics :  write sediment output files 
    55   !!====================================================================== 
    6    USE par_sed 
    76   USE sed 
    87   USE sedarr 
    98   USE lib_mpp         ! distribued memory computing library 
    109   USE iom 
     10   USE lbclnk          ! ocean lateral boundary conditions (or mpp link) 
    1111 
    1212   IMPLICIT NONE 
     
    3838      CHARACTER(len = 20)  ::  cltra  
    3939      REAL(wp)  :: zrate 
    40       REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zdta, zflx 
     40      REAL(wp), DIMENSION(jpoce, jpksed)     :: zdta 
     41      REAL(wp), DIMENSION(jpoce, jptrased+1) :: zflx 
     42      REAL(wp), DIMENSION(jpi, jpj, jpksed, jptrased)   :: trcsedi 
     43      REAL(wp), DIMENSION(jpi, jpj, jpksed, jpdia3dsed) :: flxsedi3d 
     44      REAL(wp), DIMENSION(jpi, jpj, jpdia2dsed) :: flxsedi2d 
    4145 
    4246      !!------------------------------------------------------------------- 
     
    5458      IF (lwp) WRITE(numsed,*) ' ' 
    5559       
    56       ALLOCATE( zdta(jpoce,jpksed) )    ;   ALLOCATE( zflx(jpoce,jpwatp1) ) 
    57  
    5860      ! Initialize variables 
    5961      ! -------------------- 
     
    8890      CALL unpack_arr( jpoce, flxsedi3d(1:jpi,1:jpj,1:jpksed,2)  , iarroce(1:jpoce), & 
    8991         &                   co3por(1:jpoce,1:jpksed)  ) 
     92 
     93      CALL unpack_arr( jpoce, flxsedi3d(1:jpi,1:jpj,1:jpksed,3)  , iarroce(1:jpoce), & 
     94         &                   saturco3(1:jpoce,1:jpksed)  ) 
     95 
    9096       
    9197!      flxsedi3d = 0. 
     
    95101         DO ji = 1, jpoce 
    96102            zflx(ji,jw) = ( pwcp(ji,1,jw) - pwcp_dta(ji,jw) ) & 
    97                &         * 1.e3 / 1.e2 * dzkbot(ji) / rDt_trc 
     103               &         * 1.e3 * ( 1.e-2 * dzkbot(ji) ) / 1.E4 / rDt_trc 
     104         ENDDO 
     105      ENDDO 
     106 
     107      ! Calculation of fluxes g/cm2/s 
     108      DO js = 1, jpsol 
     109         zrate =  1.0 / rDt_trc 
     110         DO ji = 1, jpoce 
     111            zflx(ji,jpwat+js) = zflx(ji,jpwat+js) + ( tosed(ji,js) - fromsed(ji,js) ) * zrate 
    98112         ENDDO 
    99113      ENDDO 
     
    101115      ! Calculation of accumulation rate per dt 
    102116      DO js = 1, jpsol 
    103          zrate =  1.0 / ( denssol * por1(jpksed) ) / rDt_trc 
     117         zrate =  1.0 / rDt_trc 
    104118         DO ji = 1, jpoce 
    105             zflx(ji,jpwatp1) = zflx(ji,jpwatp1) + ( tosed(ji,js) - fromsed(ji,js) ) * zrate 
     119            zflx(ji,jptrased+1) = zflx(ji,jptrased+1) + ( tosed(ji,js) - fromsed(ji,js) ) * zrate 
    106120         ENDDO 
    107121      ENDDO 
    108122 
    109       DO jn = 1, jpdia2dsed - 1  
     123      DO jn = 1, jpdia2dsed - 2  
    110124         CALL unpack_arr( jpoce, flxsedi2d(1:jpi,1:jpj,jn), iarroce(1:jpoce), zflx(1:jpoce,jn)  ) 
    111125      END DO 
     126 
    112127      zflx(:,1) = dzdep(:) / dtsed 
    113       CALL unpack_arr( jpoce, flxsedi2d(1:jpi,1:jpj,jpdia2dsed), iarroce(1:jpoce), zflx(1:jpoce,1) ) 
     128      CALL unpack_arr( jpoce, flxsedi2d(1:jpi,1:jpj,jpdia2dsed-1), iarroce(1:jpoce), zflx(1:jpoce,1) ) 
    114129 
    115        ! Start writing data 
    116        ! --------------------- 
    117        DO jn = 1, jptrased 
    118           cltra = sedtrcd(jn) ! short title for 3D diagnostic 
    119           CALL iom_put( cltra, trcsedi(:,:,:,jn) ) 
    120        END DO 
     130      CALL unpack_arr( jpoce, flxsedi2d(1:jpi,1:jpj,jpdia2dsed), iarroce(1:jpoce), rstepros(1:jpoce) ) 
     131      ! 
     132!      CALL lbc_lnk( 'sedwri', trcsedi(:,:,:,:), 'T', 1._wp ) 
     133!      CALL lbc_lnk( 'sedwri', flxsedi3d(:,:,:,:), 'T', 1._wp ) 
     134!      CALL lbc_lnk( 'sedwri', flxsedi2d(:,:,:), 'T', 1._wp ) 
    121135 
    122        DO jn = 1, jpdia3dsed 
    123           cltra = seddia3d(jn) ! short title for 3D diagnostic 
    124           CALL iom_put( cltra, flxsedi3d(:,:,:,jn) ) 
    125        END DO 
     136      ! Start writing data 
     137      ! --------------------- 
     138      DO jn = 1, jptrased 
     139         cltra = sedtrcd(jn) ! short title for 3D diagnostic 
     140         CALL iom_put( cltra, trcsedi(:,:,:,jn) ) 
     141      END DO 
    126142 
    127        DO jn = 1, jpdia2dsed 
    128           cltra = seddia2d(jn) ! short title for 2D diagnostic 
    129           CALL iom_put( cltra, flxsedi2d(:,:,jn) ) 
    130        END DO 
     143      DO jn = 1, jpdia3dsed 
     144         cltra = seddia3d(jn) ! short title for 3D diagnostic 
     145         CALL iom_put( cltra, flxsedi3d(:,:,:,jn) ) 
     146      END DO 
    131147 
    132  
    133       DEALLOCATE( zdta )    ;   DEALLOCATE( zflx ) 
     148      DO jn = 1, jpdia2dsed 
     149         cltra = seddia2d(jn) ! short title for 2D diagnostic 
     150         CALL iom_put( cltra, flxsedi2d(:,:,jn) ) 
     151      END DO 
    134152 
    135153      IF( ln_timing )  CALL timing_stop('sed_wri') 
  • NEMO/trunk/src/TOP/PISCES/SED/trcdmp_sed.F90

    r15023 r15450  
    9393               CALL trc_dta( kt, jl, ztrcdta )   ! read tracer data at nit000 
    9494               ! 
    95                DO_2D( 1, 1, 1, 1 ) 
     95               DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) 
    9696                  ikt = mbkt(ji,jj) 
    9797                  tr(ji,jj,ikt,jn,Kbb) = ztrcdta(ji,jj,ikt) + ( tr(ji,jj,ikt,jn,Kbb) -  ztrcdta(ji,jj,ikt) )     & 
  • NEMO/trunk/src/TOP/PISCES/trcini_pisces.F90

    r14086 r15450  
    280280      ! Initialization of the sediment model 
    281281      IF( ln_sediment)   & 
    282         & CALL sed_init ! Initialization of the sediment model  
     282        & CALL sed_ini ! Initialization of the sediment model  
    283283 
    284284      CALL p4z_sed_init          ! loss of organic matter in the sediments  
Note: See TracChangeset for help on using the changeset viewer.