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 12928 for NEMO/branches/2019/dev_r11078_OSMOSIS_IMMERSE_Nurser/src/TOP/PISCES/P4Z/p4zsed.F90 – NEMO

Ignore:
Timestamp:
2020-05-14T21:46:00+02:00 (4 years ago)
Author:
smueller
Message:

Synchronizing with /NEMO/trunk@12925 (ticket #2170)

Location:
NEMO/branches/2019/dev_r11078_OSMOSIS_IMMERSE_Nurser
Files:
2 edited

Legend:

Unmodified
Added
Removed
  • NEMO/branches/2019/dev_r11078_OSMOSIS_IMMERSE_Nurser

    • Property svn:externals
      •  

        old new  
        66^/vendors/FCM@HEAD            ext/FCM 
        77^/vendors/IOIPSL@HEAD         ext/IOIPSL 
         8 
         9# SETTE 
         10^/utils/CI/sette@HEAD         sette 
  • NEMO/branches/2019/dev_r11078_OSMOSIS_IMMERSE_Nurser/src/TOP/PISCES/P4Z/p4zsed.F90

    r10788 r12928  
    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 .AND. iom_use( "Ironice" ) )   & 
    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 ) THEN 
    146             IF( knt == nrdttrc ) THEN 
    147                 IF( iom_use( "Irondep" ) )   & 
    148                 &  CALL iom_put( "Irondep", zirondep(:,:,1) * 1.e+3 * rfact2r * e3t_n(:,:,1) * tmask(:,:,1) ) ! surface downward dust depo of iron 
    149                 IF( iom_use( "pdust" ) )   & 
    150                 &  CALL iom_put( "pdust"  , dust(:,:) / ( wdust * rday )  * tmask(:,:,1) ) ! dust concentration at surface 
    151             ENDIF 
    152          ENDIF 
    153          DEALLOCATE( zsidep, zpdep, zirondep ) 
    154          !                                               
    155       ENDIF 
    156       
    157       ! Add the external input of nutrients from river 
    158       ! ---------------------------------------------------------- 
    159       IF( ln_river ) THEN 
    160          DO jj = 1, jpj 
    161             DO ji = 1, jpi 
    162                DO jk = 1, nk_rnf(ji,jj) 
    163                   tra(ji,jj,jk,jppo4) = tra(ji,jj,jk,jppo4) +  rivdip(ji,jj) * rfact2 
    164                   tra(ji,jj,jk,jpno3) = tra(ji,jj,jk,jpno3) +  rivdin(ji,jj) * rfact2 
    165                   tra(ji,jj,jk,jpfer) = tra(ji,jj,jk,jpfer) +  rivdic(ji,jj) * 5.e-5 * rfact2 
    166                   tra(ji,jj,jk,jpsil) = tra(ji,jj,jk,jpsil) +  rivdsi(ji,jj) * rfact2 
    167                   tra(ji,jj,jk,jpdic) = tra(ji,jj,jk,jpdic) +  rivdic(ji,jj) * rfact2 
    168                   tra(ji,jj,jk,jptal) = tra(ji,jj,jk,jptal) +  ( rivalk(ji,jj) - rno3 * rivdin(ji,jj) ) * rfact2 
    169                   tra(ji,jj,jk,jpdoc) = tra(ji,jj,jk,jpdoc) +  rivdoc(ji,jj) * rfact2 
    170                ENDDO 
    171             ENDDO 
    172          ENDDO 
    173          IF (ln_ligand) THEN 
    174             DO jj = 1, jpj 
    175                DO ji = 1, jpi 
    176                   DO jk = 1, nk_rnf(ji,jj) 
    177                      tra(ji,jj,jk,jplgw) = tra(ji,jj,jk,jplgw) +  rivdic(ji,jj) * 5.e-5 * rfact2 
    178                   ENDDO 
    179                ENDDO 
    180             ENDDO 
    181          ENDIF 
    182          IF( ln_p5z ) THEN 
    183             DO jj = 1, jpj 
    184                DO ji = 1, jpi 
    185                   DO jk = 1, nk_rnf(ji,jj) 
    186                      tra(ji,jj,jk,jpdop) = tra(ji,jj,jk,jpdop) + rivdop(ji,jj) * rfact2 
    187                      tra(ji,jj,jk,jpdon) = tra(ji,jj,jk,jpdon) + rivdon(ji,jj) * rfact2 
    188                   ENDDO 
    189                ENDDO 
    190             ENDDO 
    191          ENDIF 
    192       ENDIF 
    193        
    194       ! Add the external input of nutrients from nitrogen deposition 
    195       ! ---------------------------------------------------------- 
    196       IF( ln_ndepo ) THEN 
    197          tra(:,:,1,jpno3) = tra(:,:,1,jpno3) + nitdep(:,:) * rfact2 
    198          tra(:,:,1,jptal) = tra(:,:,1,jptal) - rno3 * nitdep(:,:) * rfact2 
    199       ENDIF 
    200  
    201       ! Add the external input of iron from hydrothermal vents 
    202       ! ------------------------------------------------------ 
    203       IF( ln_hydrofe ) THEN 
    204             tra(:,:,:,jpfer) = tra(:,:,:,jpfer) + hydrofe(:,:,:) * rfact2 
    205          IF( ln_ligand ) THEN 
    206             tra(:,:,:,jplgw) = tra(:,:,:,jplgw) + ( hydrofe(:,:,:) * lgw_rath ) * rfact2 
    207          ENDIF 
    208          ! 
    209          IF( lk_iomput .AND. knt == nrdttrc .AND. iom_use( "HYDR" ) )   & 
    210             &   CALL iom_put( "HYDR", hydrofe(:,:,:) * 1.e+3 * tmask(:,:,:) ) ! hydrothermal iron input 
    211       ENDIF 
    212  
    213       ! OA: Warning, the following part is necessary to avoid CFL problems above the sediments 
    214       ! -------------------------------------------------------------------- 
    215       DO jj = 1, jpj 
    216          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 
    21796            ikt  = mbkt(ji,jj) 
    218             zdep = e3t_n(ji,jj,ikt) / xstep 
     97            zdep = e3t(ji,jj,ikt,Kmm) / xstep 
    21998            zwsbio4(ji,jj) = MIN( 0.99 * zdep, wsbio4(ji,jj,ikt) ) 
    22099            zwsbio3(ji,jj) = MIN( 0.99 * zdep, wsbio3(ji,jj,ikt) ) 
    221          END DO 
    222       END DO 
    223       ! 
    224       IF( .NOT.lk_sed ) THEN 
    225 ! 
    226          ! Add the external input of iron from sediment mobilization 
    227          ! ------------------------------------------------------ 
    228          IF( ln_ironsed ) THEN 
    229                             tra(:,:,:,jpfer) = tra(:,:,:,jpfer) + ironsed(:,:,:) * rfact2 
    230             ! 
    231             IF( lk_iomput .AND. knt == nrdttrc .AND. iom_use( "Ironsed" ) )   & 
    232                &   CALL iom_put( "Ironsed", ironsed(:,:,:) * 1.e+3 * tmask(:,:,:) ) ! iron inputs from sediments 
    233          ENDIF 
     100         END_2D 
    234101 
    235102         ! Computation of the sediment denitrification proportion: The metamodel from midlleburg (2006) is being used 
    236103         ! Computation of the fraction of organic matter that is permanently buried from Dunne's model 
    237104         ! ------------------------------------------------------- 
    238          DO jj = 1, jpj 
    239             DO ji = 1, jpi 
    240               IF( tmask(ji,jj,1) == 1 ) THEN 
    241                  ikt = mbkt(ji,jj) 
    242                  zflx = (  trb(ji,jj,ikt,jpgoc) * zwsbio4(ji,jj)   & 
    243                    &     + trb(ji,jj,ikt,jppoc) * zwsbio3(ji,jj) )  * 1E3 * 1E6 / 1E4 
    244                  zflx  = LOG10( MAX( 1E-3, zflx ) ) 
    245                  zo2   = LOG10( MAX( 10. , trb(ji,jj,ikt,jpoxy) * 1E6 ) ) 
    246                  zno3  = LOG10( MAX( 1.  , trb(ji,jj,ikt,jpno3) * 1E6 * rno3 ) ) 
    247                  zdep  = LOG10( gdepw_n(ji,jj,ikt+1) ) 
    248                  zdenit2d(ji,jj) = -2.2567 - 1.185 * zflx - 0.221 * zflx**2 - 0.3995 * zno3 * zo2 + 1.25 * zno3    & 
    249                    &                + 0.4721 * zo2 - 0.0996 * zdep + 0.4256 * zflx * zo2 
    250                  zdenit2d(ji,jj) = 10.0**( zdenit2d(ji,jj) ) 
    251                    ! 
    252                  zflx = (  trb(ji,jj,ikt,jpgoc) * zwsbio4(ji,jj)   & 
    253                    &     + trb(ji,jj,ikt,jppoc) * zwsbio3(ji,jj) ) * 1E6 
    254                  zbureff(ji,jj) = 0.013 + 0.53 * zflx**2 / ( 7.0 + zflx )**2 
    255               ENDIF 
    256             END DO 
    257          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 
    258123         ! 
    259124      ENDIF 
     
    264129      IF( .NOT.lk_sed )  zrivsil = 1._wp - sedsilfrac 
    265130 
    266       DO jj = 1, jpj 
    267          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 
    268144            ikt  = mbkt(ji,jj) 
    269             zdep = xstep / e3t_n(ji,jj,ikt)  
     145            zdep = xstep / e3t(ji,jj,ikt,Kmm)  
    270146            zwsc = zwsbio4(ji,jj) * zdep 
    271             zsiloss = trb(ji,jj,ikt,jpgsi) * zwsc 
    272             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  
    273150            ! 
    274             tra(ji,jj,ikt,jpgsi) = tra(ji,jj,ikt,jpgsi) - zsiloss 
    275             tra(ji,jj,ikt,jpcal) = tra(ji,jj,ikt,jpcal) - zcaloss 
    276          END DO 
    277       END DO 
    278       ! 
    279       IF( .NOT.lk_sed ) THEN 
    280          DO jj = 1, jpj 
    281             DO ji = 1, jpi 
    282                ikt  = mbkt(ji,jj) 
    283                zdep = xstep / e3t_n(ji,jj,ikt)  
    284                zwsc = zwsbio4(ji,jj) * zdep 
    285                zsiloss = trb(ji,jj,ikt,jpgsi) * zwsc 
    286                zcaloss = trb(ji,jj,ikt,jpcal) * zwsc 
    287                tra(ji,jj,ikt,jpsil) = tra(ji,jj,ikt,jpsil) + zsiloss * zrivsil  
    288                ! 
    289                zfactcal = MIN( excess(ji,jj,ikt), 0.2 ) 
    290                zfactcal = MIN( 1., 1.3 * ( 0.2 - zfactcal ) / ( 0.4 - zfactcal ) ) 
    291                zrivalk  = sedcalfrac * zfactcal 
    292                tra(ji,jj,ikt,jptal) =  tra(ji,jj,ikt,jptal) + zcaloss * zrivalk * 2.0 
    293                tra(ji,jj,ikt,jpdic) =  tra(ji,jj,ikt,jpdic) + zcaloss * zrivalk 
    294                zsedcal(ji,jj) = (1.0 - zrivalk) * zcaloss * e3t_n(ji,jj,ikt)  
    295                zsedsi (ji,jj) = (1.0 - zrivsil) * zsiloss * e3t_n(ji,jj,ikt)  
    296             END DO 
    297          END DO 
    298       ENDIF 
    299       ! 
    300       DO jj = 1, jpj 
    301          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 
    302174            ikt  = mbkt(ji,jj) 
    303             zdep = xstep / e3t_n(ji,jj,ikt)  
     175            zdep = xstep / e3t(ji,jj,ikt,Kmm)  
    304176            zws4 = zwsbio4(ji,jj) * zdep 
    305177            zws3 = zwsbio3(ji,jj) * zdep 
    306             tra(ji,jj,ikt,jpgoc) = tra(ji,jj,ikt,jpgoc) - trb(ji,jj,ikt,jpgoc) * zws4  
    307             tra(ji,jj,ikt,jppoc) = tra(ji,jj,ikt,jppoc) - trb(ji,jj,ikt,jppoc) * zws3 
    308             tra(ji,jj,ikt,jpbfe) = tra(ji,jj,ikt,jpbfe) - trb(ji,jj,ikt,jpbfe) * zws4 
    309             tra(ji,jj,ikt,jpsfe) = tra(ji,jj,ikt,jpsfe) - trb(ji,jj,ikt,jpsfe) * zws3 
    310          END DO 
    311       END DO 
    312       ! 
    313       IF( ln_p5z ) THEN 
    314          DO jj = 1, jpj 
    315             DO ji = 1, jpi 
    316                ikt  = mbkt(ji,jj) 
    317                zdep = xstep / e3t_n(ji,jj,ikt)  
    318                zws4 = zwsbio4(ji,jj) * zdep 
    319                zws3 = zwsbio3(ji,jj) * zdep 
    320                tra(ji,jj,ikt,jpgon) = tra(ji,jj,ikt,jpgon) - trb(ji,jj,ikt,jpgon) * zws4 
    321                tra(ji,jj,ikt,jppon) = tra(ji,jj,ikt,jppon) - trb(ji,jj,ikt,jppon) * zws3 
    322                tra(ji,jj,ikt,jpgop) = tra(ji,jj,ikt,jpgop) - trb(ji,jj,ikt,jpgop) * zws4 
    323                tra(ji,jj,ikt,jppop) = tra(ji,jj,ikt,jppop) - trb(ji,jj,ikt,jppop) * zws3 
    324             END DO 
    325          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 
    326183      ENDIF 
    327184 
     
    329186         ! The 0.5 factor in zpdenit is to avoid negative NO3 concentration after 
    330187         ! denitrification in the sediments. Not very clever, but simpliest option. 
    331          DO jj = 1, jpj 
    332             DO ji = 1, jpi 
    333                ikt  = mbkt(ji,jj) 
    334                zdep = xstep / e3t_n(ji,jj,ikt)  
    335                zws4 = zwsbio4(ji,jj) * zdep 
    336                zws3 = zwsbio3(ji,jj) * zdep 
    337                zrivno3 = 1. - zbureff(ji,jj) 
    338                zwstpoc = trb(ji,jj,ikt,jpgoc) * zws4 + trb(ji,jj,ikt,jppoc) * zws3 
    339                zpdenit  = MIN( 0.5 * ( trb(ji,jj,ikt,jpno3) - rtrn ) / rdenit, zdenit2d(ji,jj) * zwstpoc * zrivno3 ) 
    340                z1pdenit = zwstpoc * zrivno3 - zpdenit 
    341                zolimit = MIN( ( trb(ji,jj,ikt,jpoxy) - rtrn ) / o2ut, z1pdenit * ( 1.- nitrfac(ji,jj,ikt) ) ) 
    342                tra(ji,jj,ikt,jpdoc) = tra(ji,jj,ikt,jpdoc) + z1pdenit - zolimit 
    343                tra(ji,jj,ikt,jppo4) = tra(ji,jj,ikt,jppo4) + zpdenit + zolimit 
    344                tra(ji,jj,ikt,jpnh4) = tra(ji,jj,ikt,jpnh4) + zpdenit + zolimit 
    345                tra(ji,jj,ikt,jpno3) = tra(ji,jj,ikt,jpno3) - rdenit * zpdenit 
    346                tra(ji,jj,ikt,jpoxy) = tra(ji,jj,ikt,jpoxy) - zolimit * o2ut 
    347                tra(ji,jj,ikt,jptal) = tra(ji,jj,ikt,jptal) + rno3 * (zolimit + (1.+rdenit) * zpdenit ) 
    348                tra(ji,jj,ikt,jpdic) = tra(ji,jj,ikt,jpdic) + zpdenit + zolimit  
    349                sdenit(ji,jj) = rdenit * zpdenit * e3t_n(ji,jj,ikt) 
    350                zsedc(ji,jj)   = (1. - zrivno3) * zwstpoc * e3t_n(ji,jj,ikt) 
    351                IF( ln_p5z ) THEN 
    352                   zwstpop              = trb(ji,jj,ikt,jpgop) * zws4 + trb(ji,jj,ikt,jppop) * zws3 
    353                   zwstpon              = trb(ji,jj,ikt,jpgon) * zws4 + trb(ji,jj,ikt,jppon) * zws3 
    354                   tra(ji,jj,ikt,jpdon) = tra(ji,jj,ikt,jpdon) + ( z1pdenit - zolimit ) * zwstpon / (zwstpoc + rtrn) 
    355                   tra(ji,jj,ikt,jpdop) = tra(ji,jj,ikt,jpdop) + ( z1pdenit - zolimit ) * zwstpop / (zwstpoc + rtrn) 
    356                ENDIF 
    357             END DO 
    358          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 
    359214       ENDIF 
    360215 
     
    368223      ENDDO 
    369224      IF( ln_p4z ) THEN 
    370          DO jk = 1, jpkm1 
    371             DO jj = 1, jpj 
    372                DO ji = 1, jpi 
    373                   !                      ! Potential nitrogen fixation dependant on temperature and iron 
    374                   ztemp = tsn(ji,jj,jk,jp_tem) 
    375                   zmudia = MAX( 0.,-0.001096*ztemp**2 + 0.057*ztemp -0.637 ) * 7.625 
    376                   !       Potential nitrogen fixation dependant on temperature and iron 
    377                   xdianh4 = trb(ji,jj,jk,jpnh4) / ( concnnh4 + trb(ji,jj,jk,jpnh4) ) 
    378                   xdiano3 = trb(ji,jj,jk,jpno3) / ( concnno3 + trb(ji,jj,jk,jpno3) ) * (1. - xdianh4) 
    379                   zlim = ( 1.- xdiano3 - xdianh4 ) 
    380                   IF( zlim <= 0.1 )   zlim = 0.01 
    381                   zfact = zlim * rfact2 
    382                   ztrfer = biron(ji,jj,jk) / ( concfediaz + biron(ji,jj,jk) ) 
    383                   ztrpo4(ji,jj,jk) = trb(ji,jj,jk,jppo4) / ( 1E-6 + trb(ji,jj,jk,jppo4) ) 
    384                   ztrdp = ztrpo4(ji,jj,jk) 
    385                   nitrpot(ji,jj,jk) =  zmudia * r1_rday * zfact * MIN( ztrfer, ztrdp ) * zlight(ji,jj,jk) 
    386                END DO 
    387             END DO 
    388          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 
    389240      ELSE       ! p5z 
    390          DO jk = 1, jpkm1 
    391             DO jj = 1, jpj 
    392                DO ji = 1, jpi 
    393                   !                      ! Potential nitrogen fixation dependant on temperature and iron 
    394                   ztemp = tsn(ji,jj,jk,jp_tem) 
    395                   zmudia = MAX( 0.,-0.001096*ztemp**2 + 0.057*ztemp -0.637 ) * 7.625 
    396                   !       Potential nitrogen fixation dependant on temperature and iron 
    397                   xdianh4 = trb(ji,jj,jk,jpnh4) / ( concnnh4 + trb(ji,jj,jk,jpnh4) ) 
    398                   xdiano3 = trb(ji,jj,jk,jpno3) / ( concnno3 + trb(ji,jj,jk,jpno3) ) * (1. - xdianh4) 
    399                   zlim = ( 1.- xdiano3 - xdianh4 ) 
    400                   IF( zlim <= 0.1 )   zlim = 0.01 
    401                   zfact = zlim * rfact2 
    402                   ztrfer = biron(ji,jj,jk) / ( concfediaz + biron(ji,jj,jk) ) 
    403                   ztrpo4(ji,jj,jk) = trb(ji,jj,jk,jppo4) / ( 1E-6 + trb(ji,jj,jk,jppo4) ) 
    404                   ztrdop(ji,jj,jk) = trb(ji,jj,jk,jpdop) / ( 1E-6 + trb(ji,jj,jk,jpdop) ) * (1. - ztrpo4(ji,jj,jk)) 
    405                   ztrdp = ztrpo4(ji,jj,jk) + ztrdop(ji,jj,jk) 
    406                   nitrpot(ji,jj,jk) =  zmudia * r1_rday * zfact * MIN( ztrfer, ztrdp ) * zlight(ji,jj,jk) 
    407                END DO 
    408             END DO 
    409          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 
    410257      ENDIF 
    411258 
     
    413260      ! ---------------------------------------- 
    414261      IF( ln_p4z ) THEN 
    415          DO jk = 1, jpkm1 
    416             DO jj = 1, jpj 
    417                DO ji = 1, jpi 
    418                   zfact = nitrpot(ji,jj,jk) * nitrfix 
    419                   tra(ji,jj,jk,jpnh4) = tra(ji,jj,jk,jpnh4) + zfact / 3.0 
    420                   tra(ji,jj,jk,jptal) = tra(ji,jj,jk,jptal) + rno3 * zfact / 3.0 
    421                   tra(ji,jj,jk,jppo4) = tra(ji,jj,jk,jppo4) - zfact * 2.0 / 3.0 
    422                   tra(ji,jj,jk,jpdoc) = tra(ji,jj,jk,jpdoc) + zfact * 1.0 / 3.0 
    423                   tra(ji,jj,jk,jppoc) = tra(ji,jj,jk,jppoc) + zfact * 1.0 / 3.0 * 2.0 / 3.0 
    424                   tra(ji,jj,jk,jpgoc) = tra(ji,jj,jk,jpgoc) + zfact * 1.0 / 3.0 * 1.0 / 3.0 
    425                   tra(ji,jj,jk,jpoxy) = tra(ji,jj,jk,jpoxy) + ( o2ut + o2nit ) * zfact * 2.0 / 3.0 + o2nit * zfact / 3.0 
    426                   tra(ji,jj,jk,jpfer) = tra(ji,jj,jk,jpfer) - 30E-6 * zfact * 1.0 / 3.0 
    427                   tra(ji,jj,jk,jpsfe) = tra(ji,jj,jk,jpsfe) + 30E-6 * zfact * 1.0 / 3.0 * 2.0 / 3.0 
    428                   tra(ji,jj,jk,jpbfe) = tra(ji,jj,jk,jpbfe) + 30E-6 * zfact * 1.0 / 3.0 * 1.0 / 3.0 
    429                   tra(ji,jj,jk,jpfer) = tra(ji,jj,jk,jpfer) + 0.002 * 4E-10 * zsoufer(ji,jj,jk) * rfact2 / rday 
    430                   tra(ji,jj,jk,jppo4) = tra(ji,jj,jk,jppo4) + concdnh4 / ( concdnh4 + trb(ji,jj,jk,jppo4) ) & 
    431                   &                     * 0.001 * trb(ji,jj,jk,jpdoc) * xstep 
    432               END DO 
    433             END DO  
    434          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 
    435278      ELSE    ! p5z 
    436          DO jk = 1, jpkm1 
    437             DO jj = 1, jpj 
    438                DO ji = 1, jpi 
    439                   zfact = nitrpot(ji,jj,jk) * nitrfix 
    440                   tra(ji,jj,jk,jpnh4) = tra(ji,jj,jk,jpnh4) + zfact / 3.0 
    441                   tra(ji,jj,jk,jptal) = tra(ji,jj,jk,jptal) + rno3 * zfact / 3.0 
    442                   tra(ji,jj,jk,jppo4) = tra(ji,jj,jk,jppo4) - 16.0 / 46.0 * zfact * ( 1.0 - 1.0 / 3.0 ) & 
    443                   &                     * ztrpo4(ji,jj,jk) / (ztrpo4(ji,jj,jk) + ztrdop(ji,jj,jk) + rtrn) 
    444                   tra(ji,jj,jk,jpdon) = tra(ji,jj,jk,jpdon) + zfact * 1.0 / 3.0 
    445                   tra(ji,jj,jk,jpdoc) = tra(ji,jj,jk,jpdoc) + zfact * 1.0 / 3.0 
    446                   tra(ji,jj,jk,jpdop) = tra(ji,jj,jk,jpdop) + 16.0 / 46.0 * zfact / 3.0  & 
    447                   &                     - 16.0 / 46.0 * zfact * ztrdop(ji,jj,jk)   & 
    448                   &                     / (ztrpo4(ji,jj,jk) + ztrdop(ji,jj,jk) + rtrn) 
    449                   tra(ji,jj,jk,jppoc) = tra(ji,jj,jk,jppoc) + zfact * 1.0 / 3.0 * 2.0 / 3.0 
    450                   tra(ji,jj,jk,jppon) = tra(ji,jj,jk,jppon) + zfact * 1.0 / 3.0 * 2.0 /3.0 
    451                   tra(ji,jj,jk,jppop) = tra(ji,jj,jk,jppop) + 16.0 / 46.0 * zfact * 1.0 / 3.0 * 2.0 /3.0 
    452                   tra(ji,jj,jk,jpgoc) = tra(ji,jj,jk,jpgoc) + zfact * 1.0 / 3.0 * 1.0 / 3.0 
    453                   tra(ji,jj,jk,jpgon) = tra(ji,jj,jk,jpgon) + zfact * 1.0 / 3.0 * 1.0 /3.0 
    454                   tra(ji,jj,jk,jpgop) = tra(ji,jj,jk,jpgop) + 16.0 / 46.0 * zfact * 1.0 / 3.0 * 1.0 /3.0 
    455                   tra(ji,jj,jk,jpoxy) = tra(ji,jj,jk,jpoxy) + ( o2ut + o2nit ) * zfact * 2.0 / 3.0 + o2nit * zfact / 3.0 
    456                   tra(ji,jj,jk,jpfer) = tra(ji,jj,jk,jpfer) - 30E-6 * zfact * 1.0 / 3.0  
    457                   tra(ji,jj,jk,jpsfe) = tra(ji,jj,jk,jpsfe) + 30E-6 * zfact * 1.0 / 3.0 * 2.0 / 3.0 
    458                   tra(ji,jj,jk,jpbfe) = tra(ji,jj,jk,jpbfe) + 30E-6 * zfact * 1.0 / 3.0 * 1.0 / 3.0 
    459                   tra(ji,jj,jk,jpfer) = tra(ji,jj,jk,jpfer) + 0.002 * 4E-10 * zsoufer(ji,jj,jk) * rfact2 / rday 
    460               END DO 
    461             END DO  
    462          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 
    463302         ! 
    464303      ENDIF 
    465304 
    466       IF( lk_iomput ) THEN 
    467          IF( knt == nrdttrc ) THEN 
    468             zfact = 1.e+3 * rfact2r !  conversion from molC/l/kt  to molN/m3/s 
    469             IF( iom_use("Nfix"   ) ) CALL iom_put( "Nfix", nitrpot(:,:,:) * nitrfix * rno3 * zfact * tmask(:,:,:) )  ! nitrogen fixation  
    470             IF( iom_use("INTNFIX") ) THEN   ! nitrogen fixation rate in ocean ( vertically integrated ) 
    471                zwork(:,:) = 0. 
    472                DO jk = 1, jpkm1 
    473                  zwork(:,:) = zwork(:,:) + nitrpot(:,:,jk) * nitrfix * rno3 * zfact * e3t_n(:,:,jk) * tmask(:,:,jk) 
    474                ENDDO 
    475                CALL iom_put( "INTNFIX" , zwork )  
    476             ENDIF 
    477             IF( iom_use("SedCal" ) ) CALL iom_put( "SedCal", zsedcal(:,:) * zfact ) 
    478             IF( iom_use("SedSi" ) )  CALL iom_put( "SedSi",  zsedsi (:,:) * zfact ) 
    479             IF( iom_use("SedC" ) )   CALL iom_put( "SedC",   zsedc  (:,:) * zfact ) 
    480             IF( iom_use("Sdenit" ) ) CALL iom_put( "Sdenit", sdenit (:,:) * zfact * rno3 ) 
    481          ENDIF 
    482       ENDIF 
    483       ! 
    484       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) 
    485315         WRITE(charout, fmt="('sed ')") 
    486316         CALL prt_ctl_trc_info(charout) 
    487          CALL prt_ctl_trc(tab4d=tra, mask=tmask, clinfo=ctrcnm) 
     317         CALL prt_ctl_trc(tab4d=tr(:,:,:,:,Krhs), mask=tmask, clinfo=ctrcnm) 
    488318      ENDIF 
    489319      ! 
     
    494324   END SUBROUTINE p4z_sed 
    495325 
     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 
    496368 
    497369   INTEGER FUNCTION p4z_sed_alloc() 
Note: See TracChangeset for help on using the changeset viewer.