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 12377 for NEMO/trunk/src/TOP/PISCES/P4Z/p4zsed.F90 – NEMO

Ignore:
Timestamp:
2020-02-12T15:39:06+01:00 (4 years ago)
Author:
acc
Message:

The big one. Merging all 2019 developments from the option 1 branch back onto the trunk.

This changeset reproduces 2019/dev_r11943_MERGE_2019 on the trunk using a 2-URL merge
onto a working copy of the trunk. I.e.:

svn merge --ignore-ancestry \

​svn+ssh://acc@forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/NEMO/trunk \
​svn+ssh://acc@forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/NEMO/branches/2019/dev_r11943_MERGE_2019 ./

The --ignore-ancestry flag avoids problems that may otherwise arise from the fact that
the merge history been trunk and branch may have been applied in a different order but
care has been taken before this step to ensure that all applicable fixes and updates
are present in the merge branch.

The trunk state just before this step has been branched to releases/release-4.0-HEAD
and that branch has been immediately tagged as releases/release-4.0.2. Any fixes
or additions in response to tickets on 4.0, 4.0.1 or 4.0.2 should be done on
releases/release-4.0-HEAD. From now on future 'point' releases (e.g. 4.0.2) will
remain unchanged with periodic releases as needs demand. Note release-4.0-HEAD is a
transitional naming convention. Future full releases, say 4.2, will have a release-4.2
branch which fulfills this role and the first point release (e.g. 4.2.0) will be made
immediately following the release branch creation.

2020 developments can be started from any trunk revision later than this one.

Location:
NEMO/trunk
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • NEMO/trunk

    • Property svn:externals

      •  

        old new  
        33^/utils/build/mk@HEAD         mk 
        44^/utils/tools@HEAD            tools 
        5 ^/vendors/AGRIF/dev@HEAD      ext/AGRIF 
         5^/vendors/AGRIF/dev_r11615_ENHANCE-04_namelists_as_internalfiles_agrif@HEAD      ext/AGRIF 
        66^/vendors/FCM@HEAD            ext/FCM 
        77^/vendors/IOIPSL@HEAD         ext/IOIPSL 

  • NEMO/trunk/src/TOP/PISCES/P4Z/p4zsed.F90

    r12276 r12377  
    1515   USE sms_pisces      !  PISCES Source Minus Sink variables 
    1616   USE p4zlim          !  Co-limitations of differents nutrients 
    17    USE p4zsbc          !  External source of nutrients  
    1817   USE p4zint          !  interpolation and computation of various fields 
    1918   USE sed             !  Sediment module 
    … …  
    2524 
    2625   PUBLIC   p4z_sed   
     26   PUBLIC   p4z_sed_init 
    2727   PUBLIC   p4z_sed_alloc 
    2828  
     29   REAL(wp), PUBLIC ::   nitrfix      !: Nitrogen fixation rate 
     30   REAL(wp), PUBLIC ::   diazolight   !: Nitrogen fixation sensitivty to light 
     31   REAL(wp), PUBLIC ::   concfediaz   !: Fe half-saturation Cste for diazotrophs 
     32 
    2933   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: nitrpot    !: Nitrogen fixation  
    3034   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:  ) :: sdenit     !: Nitrate reduction in the sediments 
    31    REAL(wp) :: r1_rday                  !: inverse of rday 
    32    LOGICAL, SAVE :: lk_sed 
    33  
     35   ! 
     36   REAL(wp), SAVE :: r1_rday           
     37   REAL(wp), SAVE :: sedsilfrac, sedcalfrac 
     38 
     39   !! * Substitutions 
     40#  include "do_loop_substitute.h90" 
    3441   !!---------------------------------------------------------------------- 
    3542   !! NEMO/TOP 4.0 , NEMO Consortium (2018) 
    … …  
    3946CONTAINS 
    4047 
    41    SUBROUTINE p4z_sed( kt, knt ) 
     48   SUBROUTINE p4z_sed( kt, knt, Kbb, Kmm, Krhs ) 
    4249      !!--------------------------------------------------------------------- 
    4350      !!                     ***  ROUTINE p4z_sed  *** 
    … …  
    5158      ! 
    5259      INTEGER, INTENT(in) ::   kt, knt ! ocean time step 
     60      INTEGER, INTENT(in) ::   Kbb, Kmm, Krhs  ! time level indices 
    5361      INTEGER  ::  ji, jj, jk, ikt 
    5462      REAL(wp) ::  zrivalk, zrivsil, zrivno3 
    55       REAL(wp) ::  zwflux, zlim, zfact, zfactcal 
     63      REAL(wp) ::  zlim, zfact, zfactcal 
    5664      REAL(wp) ::  zo2, zno3, zflx, zpdenit, z1pdenit, zolimit 
    5765      REAL(wp) ::  zsiloss, zcaloss, zws3, zws4, zwsc, zdep 
    … …  
    6674      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zsoufer, zlight 
    6775      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: ztrpo4, ztrdop, zirondep, zpdep 
    68       REAL(wp), ALLOCATABLE, DIMENSION(:,:  ) :: zsidep, zironice 
    6976      !!--------------------------------------------------------------------- 
    7077      ! 
    7178      IF( ln_timing )  CALL timing_start('p4z_sed') 
    7279      ! 
    73       IF( kt == nittrc000 .AND. knt == 1 )   THEN 
    74           r1_rday  = 1. / rday 
    75           IF (ln_sediment .AND. ln_sed_2way) THEN 
    76              lk_sed = .TRUE. 
    77           ELSE 
    78              lk_sed = .FALSE. 
    79           ENDIF 
    80       ENDIF 
    81       ! 
    82       IF( kt == nittrc000 .AND. knt == 1 )   r1_rday  = 1. / rday 
    83       ! 
     80 
    8481      ! Allocate temporary workspace 
    8582      ALLOCATE( ztrpo4(jpi,jpj,jpk) ) 
    … …  
    9390      zsedc   (:,:) = 0.e0 
    9491 
    95       ! Iron input/uptake due to sea ice : Crude parameterization based on Lancelot et al. 
    96       ! ---------------------------------------------------- 
    97       IF( ln_ironice ) THEN   
    98          !                                               
    99          ALLOCATE( zironice(jpi,jpj) ) 
    100          !                                               
    101          DO jj = 1, jpj 
    102             DO ji = 1, jpi 
    103                zdep    = rfact2 / e3t_n(ji,jj,1) 
    104                zwflux  = fmmflx(ji,jj) / 1000._wp 
    105                zironice(ji,jj) =  MAX( -0.99 * trb(ji,jj,1,jpfer), -zwflux * icefeinput * zdep ) 
    106             END DO 
    107          END DO 
    108          ! 
    109          tra(:,:,1,jpfer) = tra(:,:,1,jpfer) + zironice(:,:)  
    110          !  
    111          IF( lk_iomput .AND. knt == nrdttrc )   & 
    112             &   CALL iom_put( "Ironice", zironice(:,:) * 1.e+3 * rfact2r * e3t_n(:,:,1) * tmask(:,:,1) ) ! iron flux from ice 
    113          ! 
    114          DEALLOCATE( zironice ) 
    115          !                                               
    116       ENDIF 
    117  
    118       ! Add the external input of nutrients from dust deposition 
    119       ! ---------------------------------------------------------- 
    120       IF( ln_dust ) THEN 
    121          !                                               
    122          ALLOCATE( zsidep(jpi,jpj), zpdep(jpi,jpj,jpk), zirondep(jpi,jpj,jpk) ) 
    123          !                                              ! Iron and Si deposition at the surface 
    124          IF( ln_solub ) THEN 
    125             zirondep(:,:,1) = solub(:,:) * dust(:,:) * mfrac * rfact2 / e3t_n(:,:,1) / 55.85 + 3.e-10 * r1_ryyss  
    126          ELSE 
    127             zirondep(:,:,1) = dustsolub  * dust(:,:) * mfrac * rfact2 / e3t_n(:,:,1) / 55.85 + 3.e-10 * r1_ryyss  
    128          ENDIF 
    129          zsidep(:,:)   = 8.8 * 0.075 * dust(:,:) * mfrac * rfact2 / e3t_n(:,:,1) / 28.1  
    130          zpdep (:,:,1) = 0.1 * 0.021 * dust(:,:) * mfrac * rfact2 / e3t_n(:,:,1) / 31. / po4r  
    131          !                                              ! Iron solubilization of particles in the water column 
    132          !                                              ! dust in kg/m2/s ---> 1/55.85 to put in mol/Fe ;  wdust in m/j 
    133          zwdust = 0.03 * rday / ( wdust * 55.85 ) / ( 270. * rday ) 
    134          DO jk = 2, jpkm1 
    135             zirondep(:,:,jk) = dust(:,:) * mfrac * zwdust * rfact2 * EXP( -gdept_n(:,:,jk) / 540. ) 
    136             zpdep   (:,:,jk) = zirondep(:,:,jk) * 0.023 
    137          END DO 
    138          !                                              ! Iron solubilization of particles in the water column 
    139          tra(:,:,1,jpsil) = tra(:,:,1,jpsil) + zsidep  (:,:) 
    140          DO jk = 1, jpkm1 
    141             tra(:,:,jk,jppo4) = tra(:,:,jk,jppo4) + zpdep   (:,:,jk) 
    142             tra(:,:,jk,jpfer) = tra(:,:,jk,jpfer) + zirondep(:,:,jk)  
    143          ENDDO 
    144          !  
    145          IF( lk_iomput .AND. knt == nrdttrc ) THEN 
    146              CALL iom_put( "Irondep", zirondep(:,:,1) * 1.e+3 * rfact2r * e3t_n(:,:,1) * tmask(:,:,1) ) ! surface downward dust depo of iron 
    147              CALL iom_put( "pdust"  , dust(:,:) / ( wdust * rday )  * tmask(:,:,1) ) ! dust concentration at surface 
    148          ENDIF 
    149          DEALLOCATE( zsidep, zpdep, zirondep ) 
    150          !                                               
    151       ENDIF 
    152       
    153       ! Add the external input of nutrients from river 
    154       ! ---------------------------------------------------------- 
    155       IF( ln_river ) THEN 
    156          DO jj = 1, jpj 
    157             DO ji = 1, jpi 
    158                DO jk = 1, nk_rnf(ji,jj) 
    159                   tra(ji,jj,jk,jppo4) = tra(ji,jj,jk,jppo4) +  rivdip(ji,jj) * rfact2 
    160                   tra(ji,jj,jk,jpno3) = tra(ji,jj,jk,jpno3) +  rivdin(ji,jj) * rfact2 
    161                   tra(ji,jj,jk,jpfer) = tra(ji,jj,jk,jpfer) +  rivdic(ji,jj) * 5.e-5 * rfact2 
    162                   tra(ji,jj,jk,jpsil) = tra(ji,jj,jk,jpsil) +  rivdsi(ji,jj) * rfact2 
    163                   tra(ji,jj,jk,jpdic) = tra(ji,jj,jk,jpdic) +  rivdic(ji,jj) * rfact2 
    164                   tra(ji,jj,jk,jptal) = tra(ji,jj,jk,jptal) +  ( rivalk(ji,jj) - rno3 * rivdin(ji,jj) ) * rfact2 
    165                   tra(ji,jj,jk,jpdoc) = tra(ji,jj,jk,jpdoc) +  rivdoc(ji,jj) * rfact2 
    166                ENDDO 
    167             ENDDO 
    168          ENDDO 
    169          IF (ln_ligand) THEN 
    170             DO jj = 1, jpj 
    171                DO ji = 1, jpi 
    172                   DO jk = 1, nk_rnf(ji,jj) 
    173                      tra(ji,jj,jk,jplgw) = tra(ji,jj,jk,jplgw) +  rivdic(ji,jj) * 5.e-5 * rfact2 
    174                   ENDDO 
    175                ENDDO 
    176             ENDDO 
    177          ENDIF 
    178          IF( ln_p5z ) THEN 
    179             DO jj = 1, jpj 
    180                DO ji = 1, jpi 
    181                   DO jk = 1, nk_rnf(ji,jj) 
    182                      tra(ji,jj,jk,jpdop) = tra(ji,jj,jk,jpdop) + rivdop(ji,jj) * rfact2 
    183                      tra(ji,jj,jk,jpdon) = tra(ji,jj,jk,jpdon) + rivdon(ji,jj) * rfact2 
    184                   ENDDO 
    185                ENDDO 
    186             ENDDO 
    187          ENDIF 
    188       ENDIF 
    189        
    190       ! Add the external input of nutrients from nitrogen deposition 
    191       ! ---------------------------------------------------------- 
    192       IF( ln_ndepo ) THEN 
    193          tra(:,:,1,jpno3) = tra(:,:,1,jpno3) + nitdep(:,:) * rfact2 
    194          tra(:,:,1,jptal) = tra(:,:,1,jptal) - rno3 * nitdep(:,:) * rfact2 
    195       ENDIF 
    196  
    197       ! Add the external input of iron from hydrothermal vents 
    198       ! ------------------------------------------------------ 
    199       IF( ln_hydrofe ) THEN 
    200             tra(:,:,:,jpfer) = tra(:,:,:,jpfer) + hydrofe(:,:,:) * rfact2 
    201          IF( ln_ligand ) THEN 
    202             tra(:,:,:,jplgw) = tra(:,:,:,jplgw) + ( hydrofe(:,:,:) * lgw_rath ) * rfact2 
    203          ENDIF 
    204          ! 
    205          IF( lk_iomput .AND. knt == nrdttrc )   & 
    206             &   CALL iom_put( "HYDR", hydrofe(:,:,:) * 1.e+3 * tmask(:,:,:) ) ! hydrothermal iron input 
    207       ENDIF 
    208  
    209       ! OA: Warning, the following part is necessary to avoid CFL problems above the sediments 
    210       ! -------------------------------------------------------------------- 
    211       DO jj = 1, jpj 
    212          DO ji = 1, jpi 
     92      IF( .NOT.lk_sed ) THEN 
     93         ! OA: Warning, the following part is necessary to avoid CFL problems above the sediments 
     94         ! -------------------------------------------------------------------- 
     95         DO_2D_11_11 
    21396            ikt  = mbkt(ji,jj) 
    214             zdep = e3t_n(ji,jj,ikt) / xstep 
     97            zdep = e3t(ji,jj,ikt,Kmm) / xstep 
    21598            zwsbio4(ji,jj) = MIN( 0.99 * zdep, wsbio4(ji,jj,ikt) ) 
    21699            zwsbio3(ji,jj) = MIN( 0.99 * zdep, wsbio3(ji,jj,ikt) ) 
    217          END DO 
    218       END DO 
    219       ! 
    220       IF( .NOT.lk_sed ) THEN 
    221 ! 
    222          ! Add the external input of iron from sediment mobilization 
    223          ! ------------------------------------------------------ 
    224          IF( ln_ironsed ) THEN 
    225                             tra(:,:,:,jpfer) = tra(:,:,:,jpfer) + ironsed(:,:,:) * rfact2 
    226             ! 
    227             IF( lk_iomput .AND. knt == nrdttrc )   & 
    228                &   CALL iom_put( "Ironsed", ironsed(:,:,:) * 1.e+3 * tmask(:,:,:) ) ! iron inputs from sediments 
    229          ENDIF 
     100         END_2D 
    230101 
    231102         ! Computation of the sediment denitrification proportion: The metamodel from midlleburg (2006) is being used 
    232103         ! Computation of the fraction of organic matter that is permanently buried from Dunne's model 
    233104         ! ------------------------------------------------------- 
    234          DO jj = 1, jpj 
    235             DO ji = 1, jpi 
    236               IF( tmask(ji,jj,1) == 1 ) THEN 
    237                  ikt = mbkt(ji,jj) 
    238                  zflx = (  trb(ji,jj,ikt,jpgoc) * zwsbio4(ji,jj)   & 
    239                    &     + trb(ji,jj,ikt,jppoc) * zwsbio3(ji,jj) )  * 1E3 * 1E6 / 1E4 
    240                  zflx  = LOG10( MAX( 1E-3, zflx ) ) 
    241                  zo2   = LOG10( MAX( 10. , trb(ji,jj,ikt,jpoxy) * 1E6 ) ) 
    242                  zno3  = LOG10( MAX( 1.  , trb(ji,jj,ikt,jpno3) * 1E6 * rno3 ) ) 
    243                  zdep  = LOG10( gdepw_n(ji,jj,ikt+1) ) 
    244                  zdenit2d(ji,jj) = -2.2567 - 1.185 * zflx - 0.221 * zflx**2 - 0.3995 * zno3 * zo2 + 1.25 * zno3    & 
    245                    &                + 0.4721 * zo2 - 0.0996 * zdep + 0.4256 * zflx * zo2 
    246                  zdenit2d(ji,jj) = 10.0**( zdenit2d(ji,jj) ) 
    247                    ! 
    248                  zflx = (  trb(ji,jj,ikt,jpgoc) * zwsbio4(ji,jj)   & 
    249                    &     + trb(ji,jj,ikt,jppoc) * zwsbio3(ji,jj) ) * 1E6 
    250                  zbureff(ji,jj) = 0.013 + 0.53 * zflx**2 / ( 7.0 + zflx )**2 
    251               ENDIF 
    252             END DO 
    253          END DO  
     105         DO_2D_11_11 
     106           IF( tmask(ji,jj,1) == 1 ) THEN 
     107              ikt = mbkt(ji,jj) 
     108              zflx = (  tr(ji,jj,ikt,jpgoc,Kbb) * zwsbio4(ji,jj)   & 
     109                &     + tr(ji,jj,ikt,jppoc,Kbb) * zwsbio3(ji,jj) )  * 1E3 * 1E6 / 1E4 
     110              zflx  = LOG10( MAX( 1E-3, zflx ) ) 
     111              zo2   = LOG10( MAX( 10. , tr(ji,jj,ikt,jpoxy,Kbb) * 1E6 ) ) 
     112              zno3  = LOG10( MAX( 1.  , tr(ji,jj,ikt,jpno3,Kbb) * 1E6 * rno3 ) ) 
     113              zdep  = LOG10( gdepw(ji,jj,ikt+1,Kmm) ) 
     114              zdenit2d(ji,jj) = -2.2567 - 1.185 * zflx - 0.221 * zflx**2 - 0.3995 * zno3 * zo2 + 1.25 * zno3    & 
     115                &                + 0.4721 * zo2 - 0.0996 * zdep + 0.4256 * zflx * zo2 
     116              zdenit2d(ji,jj) = 10.0**( zdenit2d(ji,jj) ) 
     117                ! 
     118              zflx = (  tr(ji,jj,ikt,jpgoc,Kbb) * zwsbio4(ji,jj)   & 
     119                &     + tr(ji,jj,ikt,jppoc,Kbb) * zwsbio3(ji,jj) ) * 1E6 
     120              zbureff(ji,jj) = 0.013 + 0.53 * zflx**2 / ( 7.0 + zflx )**2 
     121           ENDIF 
     122         END_2D 
    254123         ! 
    255124      ENDIF 
    … …  
    260129      IF( .NOT.lk_sed )  zrivsil = 1._wp - sedsilfrac 
    261130 
    262       DO jj = 1, jpj 
    263          DO ji = 1, jpi 
     131      DO_2D_11_11 
     132         ikt  = mbkt(ji,jj) 
     133         zdep = xstep / e3t(ji,jj,ikt,Kmm)  
     134         zwsc = zwsbio4(ji,jj) * zdep 
     135         zsiloss = tr(ji,jj,ikt,jpgsi,Kbb) * zwsc 
     136         zcaloss = tr(ji,jj,ikt,jpcal,Kbb) * zwsc 
     137         ! 
     138         tr(ji,jj,ikt,jpgsi,Krhs) = tr(ji,jj,ikt,jpgsi,Krhs) - zsiloss 
     139         tr(ji,jj,ikt,jpcal,Krhs) = tr(ji,jj,ikt,jpcal,Krhs) - zcaloss 
     140      END_2D 
     141      ! 
     142      IF( .NOT.lk_sed ) THEN 
     143         DO_2D_11_11 
    264144            ikt  = mbkt(ji,jj) 
    265             zdep = xstep / e3t_n(ji,jj,ikt)  
     145            zdep = xstep / e3t(ji,jj,ikt,Kmm)  
    266146            zwsc = zwsbio4(ji,jj) * zdep 
    267             zsiloss = trb(ji,jj,ikt,jpgsi) * zwsc 
    268             zcaloss = trb(ji,jj,ikt,jpcal) * zwsc 
     147            zsiloss = tr(ji,jj,ikt,jpgsi,Kbb) * zwsc 
     148            zcaloss = tr(ji,jj,ikt,jpcal,Kbb) * zwsc 
     149            tr(ji,jj,ikt,jpsil,Krhs) = tr(ji,jj,ikt,jpsil,Krhs) + zsiloss * zrivsil  
    269150            ! 
    270             tra(ji,jj,ikt,jpgsi) = tra(ji,jj,ikt,jpgsi) - zsiloss 
    271             tra(ji,jj,ikt,jpcal) = tra(ji,jj,ikt,jpcal) - zcaloss 
    272          END DO 
    273       END DO 
    274       ! 
    275       IF( .NOT.lk_sed ) THEN 
    276          DO jj = 1, jpj 
    277             DO ji = 1, jpi 
    278                ikt  = mbkt(ji,jj) 
    279                zdep = xstep / e3t_n(ji,jj,ikt)  
    280                zwsc = zwsbio4(ji,jj) * zdep 
    281                zsiloss = trb(ji,jj,ikt,jpgsi) * zwsc 
    282                zcaloss = trb(ji,jj,ikt,jpcal) * zwsc 
    283                tra(ji,jj,ikt,jpsil) = tra(ji,jj,ikt,jpsil) + zsiloss * zrivsil  
    284                ! 
    285                zfactcal = MIN( excess(ji,jj,ikt), 0.2 ) 
    286                zfactcal = MIN( 1., 1.3 * ( 0.2 - zfactcal ) / ( 0.4 - zfactcal ) ) 
    287                zrivalk  = sedcalfrac * zfactcal 
    288                tra(ji,jj,ikt,jptal) =  tra(ji,jj,ikt,jptal) + zcaloss * zrivalk * 2.0 
    289                tra(ji,jj,ikt,jpdic) =  tra(ji,jj,ikt,jpdic) + zcaloss * zrivalk 
    290                zsedcal(ji,jj) = (1.0 - zrivalk) * zcaloss * e3t_n(ji,jj,ikt)  
    291                zsedsi (ji,jj) = (1.0 - zrivsil) * zsiloss * e3t_n(ji,jj,ikt)  
    292             END DO 
    293          END DO 
    294       ENDIF 
    295       ! 
    296       DO jj = 1, jpj 
    297          DO ji = 1, jpi 
     151            zfactcal = MIN( excess(ji,jj,ikt), 0.2 ) 
     152            zfactcal = MIN( 1., 1.3 * ( 0.2 - zfactcal ) / ( 0.4 - zfactcal ) ) 
     153            zrivalk  = sedcalfrac * zfactcal 
     154            tr(ji,jj,ikt,jptal,Krhs) =  tr(ji,jj,ikt,jptal,Krhs) + zcaloss * zrivalk * 2.0 
     155            tr(ji,jj,ikt,jpdic,Krhs) =  tr(ji,jj,ikt,jpdic,Krhs) + zcaloss * zrivalk 
     156            zsedcal(ji,jj) = (1.0 - zrivalk) * zcaloss * e3t(ji,jj,ikt,Kmm)  
     157            zsedsi (ji,jj) = (1.0 - zrivsil) * zsiloss * e3t(ji,jj,ikt,Kmm)  
     158         END_2D 
     159      ENDIF 
     160      ! 
     161      DO_2D_11_11 
     162         ikt  = mbkt(ji,jj) 
     163         zdep = xstep / e3t(ji,jj,ikt,Kmm)  
     164         zws4 = zwsbio4(ji,jj) * zdep 
     165         zws3 = zwsbio3(ji,jj) * zdep 
     166         tr(ji,jj,ikt,jpgoc,Krhs) = tr(ji,jj,ikt,jpgoc,Krhs) - tr(ji,jj,ikt,jpgoc,Kbb) * zws4  
     167         tr(ji,jj,ikt,jppoc,Krhs) = tr(ji,jj,ikt,jppoc,Krhs) - tr(ji,jj,ikt,jppoc,Kbb) * zws3 
     168         tr(ji,jj,ikt,jpbfe,Krhs) = tr(ji,jj,ikt,jpbfe,Krhs) - tr(ji,jj,ikt,jpbfe,Kbb) * zws4 
     169         tr(ji,jj,ikt,jpsfe,Krhs) = tr(ji,jj,ikt,jpsfe,Krhs) - tr(ji,jj,ikt,jpsfe,Kbb) * zws3 
     170      END_2D 
     171      ! 
     172      IF( ln_p5z ) THEN 
     173         DO_2D_11_11 
    298174            ikt  = mbkt(ji,jj) 
    299             zdep = xstep / e3t_n(ji,jj,ikt)  
     175            zdep = xstep / e3t(ji,jj,ikt,Kmm)  
    300176            zws4 = zwsbio4(ji,jj) * zdep 
    301177            zws3 = zwsbio3(ji,jj) * zdep 
    302             tra(ji,jj,ikt,jpgoc) = tra(ji,jj,ikt,jpgoc) - trb(ji,jj,ikt,jpgoc) * zws4  
    303             tra(ji,jj,ikt,jppoc) = tra(ji,jj,ikt,jppoc) - trb(ji,jj,ikt,jppoc) * zws3 
    304             tra(ji,jj,ikt,jpbfe) = tra(ji,jj,ikt,jpbfe) - trb(ji,jj,ikt,jpbfe) * zws4 
    305             tra(ji,jj,ikt,jpsfe) = tra(ji,jj,ikt,jpsfe) - trb(ji,jj,ikt,jpsfe) * zws3 
    306          END DO 
    307       END DO 
    308       ! 
    309       IF( ln_p5z ) THEN 
    310          DO jj = 1, jpj 
    311             DO ji = 1, jpi 
    312                ikt  = mbkt(ji,jj) 
    313                zdep = xstep / e3t_n(ji,jj,ikt)  
    314                zws4 = zwsbio4(ji,jj) * zdep 
    315                zws3 = zwsbio3(ji,jj) * zdep 
    316                tra(ji,jj,ikt,jpgon) = tra(ji,jj,ikt,jpgon) - trb(ji,jj,ikt,jpgon) * zws4 
    317                tra(ji,jj,ikt,jppon) = tra(ji,jj,ikt,jppon) - trb(ji,jj,ikt,jppon) * zws3 
    318                tra(ji,jj,ikt,jpgop) = tra(ji,jj,ikt,jpgop) - trb(ji,jj,ikt,jpgop) * zws4 
    319                tra(ji,jj,ikt,jppop) = tra(ji,jj,ikt,jppop) - trb(ji,jj,ikt,jppop) * zws3 
    320             END DO 
    321          END DO 
     178            tr(ji,jj,ikt,jpgon,Krhs) = tr(ji,jj,ikt,jpgon,Krhs) - tr(ji,jj,ikt,jpgon,Kbb) * zws4 
     179            tr(ji,jj,ikt,jppon,Krhs) = tr(ji,jj,ikt,jppon,Krhs) - tr(ji,jj,ikt,jppon,Kbb) * zws3 
     180            tr(ji,jj,ikt,jpgop,Krhs) = tr(ji,jj,ikt,jpgop,Krhs) - tr(ji,jj,ikt,jpgop,Kbb) * zws4 
     181            tr(ji,jj,ikt,jppop,Krhs) = tr(ji,jj,ikt,jppop,Krhs) - tr(ji,jj,ikt,jppop,Kbb) * zws3 
     182         END_2D 
    322183      ENDIF 
    323184 
    … …  
    325186         ! The 0.5 factor in zpdenit is to avoid negative NO3 concentration after 
    326187         ! denitrification in the sediments. Not very clever, but simpliest option. 
    327          DO jj = 1, jpj 
    328             DO ji = 1, jpi 
    329                ikt  = mbkt(ji,jj) 
    330                zdep = xstep / e3t_n(ji,jj,ikt)  
    331                zws4 = zwsbio4(ji,jj) * zdep 
    332                zws3 = zwsbio3(ji,jj) * zdep 
    333                zrivno3 = 1. - zbureff(ji,jj) 
    334                zwstpoc = trb(ji,jj,ikt,jpgoc) * zws4 + trb(ji,jj,ikt,jppoc) * zws3 
    335                zpdenit  = MIN( 0.5 * ( trb(ji,jj,ikt,jpno3) - rtrn ) / rdenit, zdenit2d(ji,jj) * zwstpoc * zrivno3 ) 
    336                z1pdenit = zwstpoc * zrivno3 - zpdenit 
    337                zolimit = MIN( ( trb(ji,jj,ikt,jpoxy) - rtrn ) / o2ut, z1pdenit * ( 1.- nitrfac(ji,jj,ikt) ) ) 
    338                tra(ji,jj,ikt,jpdoc) = tra(ji,jj,ikt,jpdoc) + z1pdenit - zolimit 
    339                tra(ji,jj,ikt,jppo4) = tra(ji,jj,ikt,jppo4) + zpdenit + zolimit 
    340                tra(ji,jj,ikt,jpnh4) = tra(ji,jj,ikt,jpnh4) + zpdenit + zolimit 
    341                tra(ji,jj,ikt,jpno3) = tra(ji,jj,ikt,jpno3) - rdenit * zpdenit 
    342                tra(ji,jj,ikt,jpoxy) = tra(ji,jj,ikt,jpoxy) - zolimit * o2ut 
    343                tra(ji,jj,ikt,jptal) = tra(ji,jj,ikt,jptal) + rno3 * (zolimit + (1.+rdenit) * zpdenit ) 
    344                tra(ji,jj,ikt,jpdic) = tra(ji,jj,ikt,jpdic) + zpdenit + zolimit  
    345                sdenit(ji,jj) = rdenit * zpdenit * e3t_n(ji,jj,ikt) 
    346                zsedc(ji,jj)   = (1. - zrivno3) * zwstpoc * e3t_n(ji,jj,ikt) 
    347                IF( ln_p5z ) THEN 
    348                   zwstpop              = trb(ji,jj,ikt,jpgop) * zws4 + trb(ji,jj,ikt,jppop) * zws3 
    349                   zwstpon              = trb(ji,jj,ikt,jpgon) * zws4 + trb(ji,jj,ikt,jppon) * zws3 
    350                   tra(ji,jj,ikt,jpdon) = tra(ji,jj,ikt,jpdon) + ( z1pdenit - zolimit ) * zwstpon / (zwstpoc + rtrn) 
    351                   tra(ji,jj,ikt,jpdop) = tra(ji,jj,ikt,jpdop) + ( z1pdenit - zolimit ) * zwstpop / (zwstpoc + rtrn) 
    352                ENDIF 
    353             END DO 
    354          END DO 
     188         DO_2D_11_11 
     189            ikt  = mbkt(ji,jj) 
     190            zdep = xstep / e3t(ji,jj,ikt,Kmm)  
     191            zws4 = zwsbio4(ji,jj) * zdep 
     192            zws3 = zwsbio3(ji,jj) * zdep 
     193            zrivno3 = 1. - zbureff(ji,jj) 
     194            zwstpoc = tr(ji,jj,ikt,jpgoc,Kbb) * zws4 + tr(ji,jj,ikt,jppoc,Kbb) * zws3 
     195            zpdenit  = MIN( 0.5 * ( tr(ji,jj,ikt,jpno3,Kbb) - rtrn ) / rdenit, zdenit2d(ji,jj) * zwstpoc * zrivno3 ) 
     196            z1pdenit = zwstpoc * zrivno3 - zpdenit 
     197            zolimit = MIN( ( tr(ji,jj,ikt,jpoxy,Kbb) - rtrn ) / o2ut, z1pdenit * ( 1.- nitrfac(ji,jj,ikt) ) ) 
     198            tr(ji,jj,ikt,jpdoc,Krhs) = tr(ji,jj,ikt,jpdoc,Krhs) + z1pdenit - zolimit 
     199            tr(ji,jj,ikt,jppo4,Krhs) = tr(ji,jj,ikt,jppo4,Krhs) + zpdenit + zolimit 
     200            tr(ji,jj,ikt,jpnh4,Krhs) = tr(ji,jj,ikt,jpnh4,Krhs) + zpdenit + zolimit 
     201            tr(ji,jj,ikt,jpno3,Krhs) = tr(ji,jj,ikt,jpno3,Krhs) - rdenit * zpdenit 
     202            tr(ji,jj,ikt,jpoxy,Krhs) = tr(ji,jj,ikt,jpoxy,Krhs) - zolimit * o2ut 
     203            tr(ji,jj,ikt,jptal,Krhs) = tr(ji,jj,ikt,jptal,Krhs) + rno3 * (zolimit + (1.+rdenit) * zpdenit ) 
     204            tr(ji,jj,ikt,jpdic,Krhs) = tr(ji,jj,ikt,jpdic,Krhs) + zpdenit + zolimit  
     205            sdenit(ji,jj) = rdenit * zpdenit * e3t(ji,jj,ikt,Kmm) 
     206            zsedc(ji,jj)   = (1. - zrivno3) * zwstpoc * e3t(ji,jj,ikt,Kmm) 
     207            IF( ln_p5z ) THEN 
     208               zwstpop              = tr(ji,jj,ikt,jpgop,Kbb) * zws4 + tr(ji,jj,ikt,jppop,Kbb) * zws3 
     209               zwstpon              = tr(ji,jj,ikt,jpgon,Kbb) * zws4 + tr(ji,jj,ikt,jppon,Kbb) * zws3 
     210               tr(ji,jj,ikt,jpdon,Krhs) = tr(ji,jj,ikt,jpdon,Krhs) + ( z1pdenit - zolimit ) * zwstpon / (zwstpoc + rtrn) 
     211               tr(ji,jj,ikt,jpdop,Krhs) = tr(ji,jj,ikt,jpdop,Krhs) + ( z1pdenit - zolimit ) * zwstpop / (zwstpoc + rtrn) 
     212            ENDIF 
     213         END_2D 
    355214       ENDIF 
    356215 
    … …  
    364223      ENDDO 
    365224      IF( ln_p4z ) THEN 
    366          DO jk = 1, jpkm1 
    367             DO jj = 1, jpj 
    368                DO ji = 1, jpi 
    369                   !                      ! Potential nitrogen fixation dependant on temperature and iron 
    370                   ztemp = tsn(ji,jj,jk,jp_tem) 
    371                   zmudia = MAX( 0.,-0.001096*ztemp**2 + 0.057*ztemp -0.637 ) * 7.625 
    372                   !       Potential nitrogen fixation dependant on temperature and iron 
    373                   xdianh4 = trb(ji,jj,jk,jpnh4) / ( concnnh4 + trb(ji,jj,jk,jpnh4) ) 
    374                   xdiano3 = trb(ji,jj,jk,jpno3) / ( concnno3 + trb(ji,jj,jk,jpno3) ) * (1. - xdianh4) 
    375                   zlim = ( 1.- xdiano3 - xdianh4 ) 
    376                   IF( zlim <= 0.1 )   zlim = 0.01 
    377                   zfact = zlim * rfact2 
    378                   ztrfer = biron(ji,jj,jk) / ( concfediaz + biron(ji,jj,jk) ) 
    379                   ztrpo4(ji,jj,jk) = trb(ji,jj,jk,jppo4) / ( 1E-6 + trb(ji,jj,jk,jppo4) ) 
    380                   ztrdp = ztrpo4(ji,jj,jk) 
    381                   nitrpot(ji,jj,jk) =  zmudia * r1_rday * zfact * MIN( ztrfer, ztrdp ) * zlight(ji,jj,jk) 
    382                END DO 
    383             END DO 
    384          END DO 
     225         DO_3D_11_11( 1, jpkm1 ) 
     226            !                      ! Potential nitrogen fixation dependant on temperature and iron 
     227            ztemp = ts(ji,jj,jk,jp_tem,Kmm) 
     228            zmudia = MAX( 0.,-0.001096*ztemp**2 + 0.057*ztemp -0.637 ) * 7.625 
     229            !       Potential nitrogen fixation dependant on temperature and iron 
     230            xdianh4 = tr(ji,jj,jk,jpnh4,Kbb) / ( concnnh4 + tr(ji,jj,jk,jpnh4,Kbb) ) 
     231            xdiano3 = tr(ji,jj,jk,jpno3,Kbb) / ( concnno3 + tr(ji,jj,jk,jpno3,Kbb) ) * (1. - xdianh4) 
     232            zlim = ( 1.- xdiano3 - xdianh4 ) 
     233            IF( zlim <= 0.1 )   zlim = 0.01 
     234            zfact = zlim * rfact2 
     235            ztrfer = biron(ji,jj,jk) / ( concfediaz + biron(ji,jj,jk) ) 
     236            ztrpo4(ji,jj,jk) = tr(ji,jj,jk,jppo4,Kbb) / ( 1E-6 + tr(ji,jj,jk,jppo4,Kbb) ) 
     237            ztrdp = ztrpo4(ji,jj,jk) 
     238            nitrpot(ji,jj,jk) =  zmudia * r1_rday * zfact * MIN( ztrfer, ztrdp ) * zlight(ji,jj,jk) 
     239         END_3D 
    385240      ELSE       ! p5z 
    386          DO jk = 1, jpkm1 
    387             DO jj = 1, jpj 
    388                DO ji = 1, jpi 
    389                   !                      ! Potential nitrogen fixation dependant on temperature and iron 
    390                   ztemp = tsn(ji,jj,jk,jp_tem) 
    391                   zmudia = MAX( 0.,-0.001096*ztemp**2 + 0.057*ztemp -0.637 ) * 7.625 
    392                   !       Potential nitrogen fixation dependant on temperature and iron 
    393                   xdianh4 = trb(ji,jj,jk,jpnh4) / ( concnnh4 + trb(ji,jj,jk,jpnh4) ) 
    394                   xdiano3 = trb(ji,jj,jk,jpno3) / ( concnno3 + trb(ji,jj,jk,jpno3) ) * (1. - xdianh4) 
    395                   zlim = ( 1.- xdiano3 - xdianh4 ) 
    396                   IF( zlim <= 0.1 )   zlim = 0.01 
    397                   zfact = zlim * rfact2 
    398                   ztrfer = biron(ji,jj,jk) / ( concfediaz + biron(ji,jj,jk) ) 
    399                   ztrpo4(ji,jj,jk) = trb(ji,jj,jk,jppo4) / ( 1E-6 + trb(ji,jj,jk,jppo4) ) 
    400                   ztrdop(ji,jj,jk) = trb(ji,jj,jk,jpdop) / ( 1E-6 + trb(ji,jj,jk,jpdop) ) * (1. - ztrpo4(ji,jj,jk)) 
    401                   ztrdp = ztrpo4(ji,jj,jk) + ztrdop(ji,jj,jk) 
    402                   nitrpot(ji,jj,jk) =  zmudia * r1_rday * zfact * MIN( ztrfer, ztrdp ) * zlight(ji,jj,jk) 
    403                END DO 
    404             END DO 
    405          END DO 
     241         DO_3D_11_11( 1, jpkm1 ) 
     242            !                      ! Potential nitrogen fixation dependant on temperature and iron 
     243            ztemp = ts(ji,jj,jk,jp_tem,Kmm) 
     244            zmudia = MAX( 0.,-0.001096*ztemp**2 + 0.057*ztemp -0.637 ) * 7.625 
     245            !       Potential nitrogen fixation dependant on temperature and iron 
     246            xdianh4 = tr(ji,jj,jk,jpnh4,Kbb) / ( concnnh4 + tr(ji,jj,jk,jpnh4,Kbb) ) 
     247            xdiano3 = tr(ji,jj,jk,jpno3,Kbb) / ( concnno3 + tr(ji,jj,jk,jpno3,Kbb) ) * (1. - xdianh4) 
     248            zlim = ( 1.- xdiano3 - xdianh4 ) 
     249            IF( zlim <= 0.1 )   zlim = 0.01 
     250            zfact = zlim * rfact2 
     251            ztrfer = biron(ji,jj,jk) / ( concfediaz + biron(ji,jj,jk) ) 
     252            ztrpo4(ji,jj,jk) = tr(ji,jj,jk,jppo4,Kbb) / ( 1E-6 + tr(ji,jj,jk,jppo4,Kbb) ) 
     253            ztrdop(ji,jj,jk) = tr(ji,jj,jk,jpdop,Kbb) / ( 1E-6 + tr(ji,jj,jk,jpdop,Kbb) ) * (1. - ztrpo4(ji,jj,jk)) 
     254            ztrdp = ztrpo4(ji,jj,jk) + ztrdop(ji,jj,jk) 
     255            nitrpot(ji,jj,jk) =  zmudia * r1_rday * zfact * MIN( ztrfer, ztrdp ) * zlight(ji,jj,jk) 
     256         END_3D 
    406257      ENDIF 
    407258 
    … …  
    409260      ! ---------------------------------------- 
    410261      IF( ln_p4z ) THEN 
    411          DO jk = 1, jpkm1 
    412             DO jj = 1, jpj 
    413                DO ji = 1, jpi 
    414                   zfact = nitrpot(ji,jj,jk) * nitrfix 
    415                   tra(ji,jj,jk,jpnh4) = tra(ji,jj,jk,jpnh4) + zfact / 3.0 
    416                   tra(ji,jj,jk,jptal) = tra(ji,jj,jk,jptal) + rno3 * zfact / 3.0 
    417                   tra(ji,jj,jk,jppo4) = tra(ji,jj,jk,jppo4) - zfact * 2.0 / 3.0 
    418                   tra(ji,jj,jk,jpdoc) = tra(ji,jj,jk,jpdoc) + zfact * 1.0 / 3.0 
    419                   tra(ji,jj,jk,jppoc) = tra(ji,jj,jk,jppoc) + zfact * 1.0 / 3.0 * 2.0 / 3.0 
    420                   tra(ji,jj,jk,jpgoc) = tra(ji,jj,jk,jpgoc) + zfact * 1.0 / 3.0 * 1.0 / 3.0 
    421                   tra(ji,jj,jk,jpoxy) = tra(ji,jj,jk,jpoxy) + ( o2ut + o2nit ) * zfact * 2.0 / 3.0 + o2nit * zfact / 3.0 
    422                   tra(ji,jj,jk,jpfer) = tra(ji,jj,jk,jpfer) - 30E-6 * zfact * 1.0 / 3.0 
    423                   tra(ji,jj,jk,jpsfe) = tra(ji,jj,jk,jpsfe) + 30E-6 * zfact * 1.0 / 3.0 * 2.0 / 3.0 
    424                   tra(ji,jj,jk,jpbfe) = tra(ji,jj,jk,jpbfe) + 30E-6 * zfact * 1.0 / 3.0 * 1.0 / 3.0 
    425                   tra(ji,jj,jk,jpfer) = tra(ji,jj,jk,jpfer) + 0.002 * 4E-10 * zsoufer(ji,jj,jk) * rfact2 / rday 
    426                   tra(ji,jj,jk,jppo4) = tra(ji,jj,jk,jppo4) + concdnh4 / ( concdnh4 + trb(ji,jj,jk,jppo4) ) & 
    427                   &                     * 0.001 * trb(ji,jj,jk,jpdoc) * xstep 
    428               END DO 
    429             END DO  
    430          END DO 
     262         DO_3D_11_11( 1, jpkm1 ) 
     263            zfact = nitrpot(ji,jj,jk) * nitrfix 
     264            tr(ji,jj,jk,jpnh4,Krhs) = tr(ji,jj,jk,jpnh4,Krhs) + zfact / 3.0 
     265            tr(ji,jj,jk,jptal,Krhs) = tr(ji,jj,jk,jptal,Krhs) + rno3 * zfact / 3.0 
     266            tr(ji,jj,jk,jppo4,Krhs) = tr(ji,jj,jk,jppo4,Krhs) - zfact * 2.0 / 3.0 
     267            tr(ji,jj,jk,jpdoc,Krhs) = tr(ji,jj,jk,jpdoc,Krhs) + zfact * 1.0 / 3.0 
     268            tr(ji,jj,jk,jppoc,Krhs) = tr(ji,jj,jk,jppoc,Krhs) + zfact * 1.0 / 3.0 * 2.0 / 3.0 
     269            tr(ji,jj,jk,jpgoc,Krhs) = tr(ji,jj,jk,jpgoc,Krhs) + zfact * 1.0 / 3.0 * 1.0 / 3.0 
     270            tr(ji,jj,jk,jpoxy,Krhs) = tr(ji,jj,jk,jpoxy,Krhs) + ( o2ut + o2nit ) * zfact * 2.0 / 3.0 + o2nit * zfact / 3.0 
     271            tr(ji,jj,jk,jpfer,Krhs) = tr(ji,jj,jk,jpfer,Krhs) - 30E-6 * zfact * 1.0 / 3.0 
     272            tr(ji,jj,jk,jpsfe,Krhs) = tr(ji,jj,jk,jpsfe,Krhs) + 30E-6 * zfact * 1.0 / 3.0 * 2.0 / 3.0 
     273            tr(ji,jj,jk,jpbfe,Krhs) = tr(ji,jj,jk,jpbfe,Krhs) + 30E-6 * zfact * 1.0 / 3.0 * 1.0 / 3.0 
     274            tr(ji,jj,jk,jpfer,Krhs) = tr(ji,jj,jk,jpfer,Krhs) + 0.002 * 4E-10 * zsoufer(ji,jj,jk) * rfact2 / rday 
     275            tr(ji,jj,jk,jppo4,Krhs) = tr(ji,jj,jk,jppo4,Krhs) + concdnh4 / ( concdnh4 + tr(ji,jj,jk,jppo4,Kbb) ) & 
     276            &                     * 0.001 * tr(ji,jj,jk,jpdoc,Kbb) * xstep 
     277         END_3D 
    431278      ELSE    ! p5z 
    432          DO jk = 1, jpkm1 
    433             DO jj = 1, jpj 
    434                DO ji = 1, jpi 
    435                   zfact = nitrpot(ji,jj,jk) * nitrfix 
    436                   tra(ji,jj,jk,jpnh4) = tra(ji,jj,jk,jpnh4) + zfact / 3.0 
    437                   tra(ji,jj,jk,jptal) = tra(ji,jj,jk,jptal) + rno3 * zfact / 3.0 
    438                   tra(ji,jj,jk,jppo4) = tra(ji,jj,jk,jppo4) - 16.0 / 46.0 * zfact * ( 1.0 - 1.0 / 3.0 ) & 
    439                   &                     * ztrpo4(ji,jj,jk) / (ztrpo4(ji,jj,jk) + ztrdop(ji,jj,jk) + rtrn) 
    440                   tra(ji,jj,jk,jpdon) = tra(ji,jj,jk,jpdon) + zfact * 1.0 / 3.0 
    441                   tra(ji,jj,jk,jpdoc) = tra(ji,jj,jk,jpdoc) + zfact * 1.0 / 3.0 
    442                   tra(ji,jj,jk,jpdop) = tra(ji,jj,jk,jpdop) + 16.0 / 46.0 * zfact / 3.0  & 
    443                   &                     - 16.0 / 46.0 * zfact * ztrdop(ji,jj,jk)   & 
    444                   &                     / (ztrpo4(ji,jj,jk) + ztrdop(ji,jj,jk) + rtrn) 
    445                   tra(ji,jj,jk,jppoc) = tra(ji,jj,jk,jppoc) + zfact * 1.0 / 3.0 * 2.0 / 3.0 
    446                   tra(ji,jj,jk,jppon) = tra(ji,jj,jk,jppon) + zfact * 1.0 / 3.0 * 2.0 /3.0 
    447                   tra(ji,jj,jk,jppop) = tra(ji,jj,jk,jppop) + 16.0 / 46.0 * zfact * 1.0 / 3.0 * 2.0 /3.0 
    448                   tra(ji,jj,jk,jpgoc) = tra(ji,jj,jk,jpgoc) + zfact * 1.0 / 3.0 * 1.0 / 3.0 
    449                   tra(ji,jj,jk,jpgon) = tra(ji,jj,jk,jpgon) + zfact * 1.0 / 3.0 * 1.0 /3.0 
    450                   tra(ji,jj,jk,jpgop) = tra(ji,jj,jk,jpgop) + 16.0 / 46.0 * zfact * 1.0 / 3.0 * 1.0 /3.0 
    451                   tra(ji,jj,jk,jpoxy) = tra(ji,jj,jk,jpoxy) + ( o2ut + o2nit ) * zfact * 2.0 / 3.0 + o2nit * zfact / 3.0 
    452                   tra(ji,jj,jk,jpfer) = tra(ji,jj,jk,jpfer) - 30E-6 * zfact * 1.0 / 3.0  
    453                   tra(ji,jj,jk,jpsfe) = tra(ji,jj,jk,jpsfe) + 30E-6 * zfact * 1.0 / 3.0 * 2.0 / 3.0 
    454                   tra(ji,jj,jk,jpbfe) = tra(ji,jj,jk,jpbfe) + 30E-6 * zfact * 1.0 / 3.0 * 1.0 / 3.0 
    455                   tra(ji,jj,jk,jpfer) = tra(ji,jj,jk,jpfer) + 0.002 * 4E-10 * zsoufer(ji,jj,jk) * rfact2 / rday 
    456               END DO 
    457             END DO  
    458          END DO 
     279         DO_3D_11_11( 1, jpkm1 ) 
     280            zfact = nitrpot(ji,jj,jk) * nitrfix 
     281            tr(ji,jj,jk,jpnh4,Krhs) = tr(ji,jj,jk,jpnh4,Krhs) + zfact / 3.0 
     282            tr(ji,jj,jk,jptal,Krhs) = tr(ji,jj,jk,jptal,Krhs) + rno3 * zfact / 3.0 
     283            tr(ji,jj,jk,jppo4,Krhs) = tr(ji,jj,jk,jppo4,Krhs) - 16.0 / 46.0 * zfact * ( 1.0 - 1.0 / 3.0 ) & 
     284            &                     * ztrpo4(ji,jj,jk) / (ztrpo4(ji,jj,jk) + ztrdop(ji,jj,jk) + rtrn) 
     285            tr(ji,jj,jk,jpdon,Krhs) = tr(ji,jj,jk,jpdon,Krhs) + zfact * 1.0 / 3.0 
     286            tr(ji,jj,jk,jpdoc,Krhs) = tr(ji,jj,jk,jpdoc,Krhs) + zfact * 1.0 / 3.0 
     287            tr(ji,jj,jk,jpdop,Krhs) = tr(ji,jj,jk,jpdop,Krhs) + 16.0 / 46.0 * zfact / 3.0  & 
     288            &                     - 16.0 / 46.0 * zfact * ztrdop(ji,jj,jk)   & 
     289            &                     / (ztrpo4(ji,jj,jk) + ztrdop(ji,jj,jk) + rtrn) 
     290            tr(ji,jj,jk,jppoc,Krhs) = tr(ji,jj,jk,jppoc,Krhs) + zfact * 1.0 / 3.0 * 2.0 / 3.0 
     291            tr(ji,jj,jk,jppon,Krhs) = tr(ji,jj,jk,jppon,Krhs) + zfact * 1.0 / 3.0 * 2.0 /3.0 
     292            tr(ji,jj,jk,jppop,Krhs) = tr(ji,jj,jk,jppop,Krhs) + 16.0 / 46.0 * zfact * 1.0 / 3.0 * 2.0 /3.0 
     293            tr(ji,jj,jk,jpgoc,Krhs) = tr(ji,jj,jk,jpgoc,Krhs) + zfact * 1.0 / 3.0 * 1.0 / 3.0 
     294            tr(ji,jj,jk,jpgon,Krhs) = tr(ji,jj,jk,jpgon,Krhs) + zfact * 1.0 / 3.0 * 1.0 /3.0 
     295            tr(ji,jj,jk,jpgop,Krhs) = tr(ji,jj,jk,jpgop,Krhs) + 16.0 / 46.0 * zfact * 1.0 / 3.0 * 1.0 /3.0 
     296            tr(ji,jj,jk,jpoxy,Krhs) = tr(ji,jj,jk,jpoxy,Krhs) + ( o2ut + o2nit ) * zfact * 2.0 / 3.0 + o2nit * zfact / 3.0 
     297            tr(ji,jj,jk,jpfer,Krhs) = tr(ji,jj,jk,jpfer,Krhs) - 30E-6 * zfact * 1.0 / 3.0  
     298            tr(ji,jj,jk,jpsfe,Krhs) = tr(ji,jj,jk,jpsfe,Krhs) + 30E-6 * zfact * 1.0 / 3.0 * 2.0 / 3.0 
     299            tr(ji,jj,jk,jpbfe,Krhs) = tr(ji,jj,jk,jpbfe,Krhs) + 30E-6 * zfact * 1.0 / 3.0 * 1.0 / 3.0 
     300            tr(ji,jj,jk,jpfer,Krhs) = tr(ji,jj,jk,jpfer,Krhs) + 0.002 * 4E-10 * zsoufer(ji,jj,jk) * rfact2 / rday 
     301         END_3D 
    459302         ! 
    460303      ENDIF 
    461304 
    462       IF( lk_iomput ) THEN 
    463          IF( knt == nrdttrc ) THEN 
    464             zfact = 1.e+3 * rfact2r !  conversion from molC/l/kt  to molN/m3/s 
    465             CALL iom_put( "Nfix", nitrpot(:,:,:) * nitrfix * rno3 * zfact * tmask(:,:,:) )  ! nitrogen fixation  
    466             CALL iom_put( "SedCal", zsedcal(:,:) * zfact ) 
    467             CALL iom_put( "SedSi",  zsedsi (:,:) * zfact ) 
    468             CALL iom_put( "SedC",   zsedc  (:,:) * zfact ) 
    469             CALL iom_put( "Sdenit", sdenit (:,:) * zfact * rno3 ) 
    470          ENDIF 
    471       ENDIF 
    472       ! 
    473       IF(ln_ctl) THEN  ! print mean trends (USEd for debugging) 
     305      IF( lk_iomput .AND. knt == nrdttrc ) THEN 
     306         zfact = 1.e+3 * rfact2r !  conversion from molC/l/kt  to molN/m3/s 
     307         CALL iom_put( "Nfix", nitrpot(:,:,:) * nitrfix * rno3 * zfact * tmask(:,:,:) )  ! nitrogen fixation  
     308         CALL iom_put( "SedCal", zsedcal(:,:) * zfact ) 
     309         CALL iom_put( "SedSi" , zsedsi (:,:) * zfact ) 
     310         CALL iom_put( "SedC"  , zsedc  (:,:) * zfact ) 
     311         CALL iom_put( "Sdenit", sdenit (:,:) * zfact * rno3 ) 
     312      ENDIF 
     313      ! 
     314      IF(sn_cfctl%l_prttrc) THEN  ! print mean trends (USEd for debugging) 
    474315         WRITE(charout, fmt="('sed ')") 
    475316         CALL prt_ctl_trc_info(charout) 
    476          CALL prt_ctl_trc(tab4d=tra, mask=tmask, clinfo=ctrcnm) 
     317         CALL prt_ctl_trc(tab4d=tr(:,:,:,:,Krhs), mask=tmask, clinfo=ctrcnm) 
    477318      ENDIF 
    478319      ! 
    … …  
    483324   END SUBROUTINE p4z_sed 
    484325 
     326   SUBROUTINE p4z_sed_init 
     327      !!---------------------------------------------------------------------- 
     328      !!                  ***  routine p4z_sed_init  *** 
     329      !! 
     330      !! ** purpose :   initialization of some parameters 
     331      !! 
     332      !!---------------------------------------------------------------------- 
     333      !!---------------------------------------------------------------------- 
     334      INTEGER  :: ji, jj, jk, jm 
     335      INTEGER  :: ios                 ! Local integer output status for namelist read 
     336      ! 
     337      !! 
     338      NAMELIST/nampissed/ nitrfix, diazolight, concfediaz 
     339      !!---------------------------------------------------------------------- 
     340      ! 
     341      IF(lwp) THEN 
     342         WRITE(numout,*) 
     343         WRITE(numout,*) 'p4z_sed_init : initialization of sediment mobilisation ' 
     344         WRITE(numout,*) '~~~~~~~~~~~~ ' 
     345      ENDIF 
     346      !                            !* set file information 
     347      READ  ( numnatp_ref, nampissed, IOSTAT = ios, ERR = 901) 
     348901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'nampissed in reference namelist' ) 
     349      READ  ( numnatp_cfg, nampissed, IOSTAT = ios, ERR = 902 ) 
     350902   IF( ios >  0 )   CALL ctl_nam ( ios , 'nampissed in configuration namelist' ) 
     351      IF(lwm) WRITE ( numonp, nampissed ) 
     352 
     353      IF(lwp) THEN 
     354         WRITE(numout,*) '   Namelist : nampissed ' 
     355         WRITE(numout,*) '      nitrogen fixation rate                       nitrfix = ', nitrfix 
     356         WRITE(numout,*) '      nitrogen fixation sensitivty to light    diazolight  = ', diazolight 
     357         WRITE(numout,*) '      Fe half-saturation cste for diazotrophs  concfediaz  = ', concfediaz 
     358      ENDIF 
     359      ! 
     360      r1_rday  = 1. / rday 
     361      ! 
     362      sedsilfrac = 0.03     ! percentage of silica loss in the sediments 
     363      sedcalfrac = 0.6      ! percentage of calcite loss in the sediments 
     364      ! 
     365      lk_sed = ln_sediment .AND. ln_sed_2way  
     366      ! 
     367   END SUBROUTINE p4z_sed_init 
    485368 
    486369   INTEGER FUNCTION p4z_sed_alloc() 
Note: See TracChangeset for help on using the changeset viewer.