Index: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/TOP/PISCES/P4Z/p4zmort.F90
===================================================================
--- NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/TOP/PISCES/P4Z/p4zmort.F90	(revision 10740)
+++ NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/TOP/PISCES/P4Z/p4zmort.F90	(revision 10975)
@@ -36,5 +36,5 @@
 CONTAINS
 
-   SUBROUTINE p4z_mort( kt )
+   SUBROUTINE p4z_mort( kt, Kbb, Krhs )
       !!---------------------------------------------------------------------
       !!                     ***  ROUTINE p4z_mort  ***
@@ -46,14 +46,15 @@
       !!---------------------------------------------------------------------
       INTEGER, INTENT(in) ::   kt ! ocean time step
-      !!---------------------------------------------------------------------
-      !
-      CALL p4z_nano            ! nanophytoplankton
-      !
-      CALL p4z_diat            ! diatoms
+      INTEGER, INTENT(in) ::   Kbb, Krhs  ! time level indices
+      !!---------------------------------------------------------------------
+      !
+      CALL p4z_nano( Kbb, Krhs )            ! nanophytoplankton
+      !
+      CALL p4z_diat( Kbb, Krhs )            ! diatoms
       !
    END SUBROUTINE p4z_mort
 
 
-   SUBROUTINE p4z_nano
+   SUBROUTINE p4z_nano( Kbb, Krhs )
       !!---------------------------------------------------------------------
       !!                     ***  ROUTINE p4z_nano  ***
@@ -63,4 +64,5 @@
       !! ** Method  : - ???
       !!---------------------------------------------------------------------
+      INTEGER, INTENT(in) ::   Kbb, Krhs  ! time level indices
       INTEGER  ::   ji, jj, jk
       REAL(wp) ::   zsizerat, zcompaph
@@ -76,10 +78,10 @@
          DO jj = 1, jpj
             DO ji = 1, jpi
-               zcompaph = MAX( ( trb(ji,jj,jk,jpphy) - 1e-8 ), 0.e0 )
+               zcompaph = MAX( ( tr(ji,jj,jk,jpphy,Kbb) - 1e-8 ), 0.e0 )
                !     When highly limited by macronutrients, very small cells 
                !     dominate the community. As a consequence, aggregation
                !     due to turbulence is negligible. Mortality is also set
                !     to 0
-               zsizerat = MIN(1., MAX( 0., (quotan(ji,jj,jk) - 0.2) / 0.3) ) * trb(ji,jj,jk,jpphy)
+               zsizerat = MIN(1., MAX( 0., (quotan(ji,jj,jk) - 0.2) / 0.3) ) * tr(ji,jj,jk,jpphy,Kbb)
                !     Squared mortality of Phyto similar to a sedimentation term during
                !     blooms (Doney et al. 1996)
@@ -89,5 +91,5 @@
                !     increased when nutrients are limiting phytoplankton growth
                !     as observed for instance in case of iron limitation.
-               ztortp = mprat * xstep * zcompaph / ( xkmort + trb(ji,jj,jk,jpphy) ) * zsizerat
+               ztortp = mprat * xstep * zcompaph / ( xkmort + tr(ji,jj,jk,jpphy,Kbb) ) * zsizerat
 
                zmortp = zrespp + ztortp
@@ -95,9 +97,9 @@
                !   Update the arrays TRA which contains the biological sources and sinks
 
-               zfactfe = trb(ji,jj,jk,jpnfe)/(trb(ji,jj,jk,jpphy)+rtrn)
-               zfactch = trb(ji,jj,jk,jpnch)/(trb(ji,jj,jk,jpphy)+rtrn)
-               tra(ji,jj,jk,jpphy) = tra(ji,jj,jk,jpphy) - zmortp
-               tra(ji,jj,jk,jpnch) = tra(ji,jj,jk,jpnch) - zmortp * zfactch
-               tra(ji,jj,jk,jpnfe) = tra(ji,jj,jk,jpnfe) - zmortp * zfactfe
+               zfactfe = tr(ji,jj,jk,jpnfe,Kbb)/(tr(ji,jj,jk,jpphy,Kbb)+rtrn)
+               zfactch = tr(ji,jj,jk,jpnch,Kbb)/(tr(ji,jj,jk,jpphy,Kbb)+rtrn)
+               tr(ji,jj,jk,jpphy,Krhs) = tr(ji,jj,jk,jpphy,Krhs) - zmortp
+               tr(ji,jj,jk,jpnch,Krhs) = tr(ji,jj,jk,jpnch,Krhs) - zmortp * zfactch
+               tr(ji,jj,jk,jpnfe,Krhs) = tr(ji,jj,jk,jpnfe,Krhs) - zmortp * zfactfe
                zprcaca = xfracal(ji,jj,jk) * zmortp
                !
@@ -105,13 +107,13 @@
                !
                zfracal = 0.5 * xfracal(ji,jj,jk)
-               tra(ji,jj,jk,jpdic) = tra(ji,jj,jk,jpdic) - zprcaca
-               tra(ji,jj,jk,jptal) = tra(ji,jj,jk,jptal) - 2. * zprcaca
-               tra(ji,jj,jk,jpcal) = tra(ji,jj,jk,jpcal) + zprcaca
-               tra(ji,jj,jk,jpgoc) = tra(ji,jj,jk,jpgoc) + zfracal * zmortp
-               tra(ji,jj,jk,jppoc) = tra(ji,jj,jk,jppoc) + ( 1. - zfracal ) * zmortp
+               tr(ji,jj,jk,jpdic,Krhs) = tr(ji,jj,jk,jpdic,Krhs) - zprcaca
+               tr(ji,jj,jk,jptal,Krhs) = tr(ji,jj,jk,jptal,Krhs) - 2. * zprcaca
+               tr(ji,jj,jk,jpcal,Krhs) = tr(ji,jj,jk,jpcal,Krhs) + zprcaca
+               tr(ji,jj,jk,jpgoc,Krhs) = tr(ji,jj,jk,jpgoc,Krhs) + zfracal * zmortp
+               tr(ji,jj,jk,jppoc,Krhs) = tr(ji,jj,jk,jppoc,Krhs) + ( 1. - zfracal ) * zmortp
                prodpoc(ji,jj,jk) = prodpoc(ji,jj,jk) + ( 1. - zfracal ) * zmortp
                prodgoc(ji,jj,jk) = prodgoc(ji,jj,jk) + zfracal * zmortp
-               tra(ji,jj,jk,jpsfe) = tra(ji,jj,jk,jpsfe) + ( 1. - zfracal ) * zmortp * zfactfe
-               tra(ji,jj,jk,jpbfe) = tra(ji,jj,jk,jpbfe) + zfracal * zmortp * zfactfe
+               tr(ji,jj,jk,jpsfe,Krhs) = tr(ji,jj,jk,jpsfe,Krhs) + ( 1. - zfracal ) * zmortp * zfactfe
+               tr(ji,jj,jk,jpbfe,Krhs) = tr(ji,jj,jk,jpbfe,Krhs) + zfracal * zmortp * zfactfe
             END DO
          END DO
@@ -121,5 +123,5 @@
          WRITE(charout, FMT="('nano')")
          CALL prt_ctl_trc_info(charout)
-         CALL prt_ctl_trc(tab4d=tra, mask=tmask, clinfo=ctrcnm)
+         CALL prt_ctl_trc(tab4d=tr(:,:,:,:,Krhs), mask=tmask, clinfo=ctrcnm)
        ENDIF
       !
@@ -129,5 +131,5 @@
 
 
-   SUBROUTINE p4z_diat
+   SUBROUTINE p4z_diat( Kbb, Krhs )
       !!---------------------------------------------------------------------
       !!                     ***  ROUTINE p4z_diat  ***
@@ -137,4 +139,5 @@
       !! ** Method  : - ???
       !!---------------------------------------------------------------------
+      INTEGER, INTENT(in) ::   Kbb, Krhs  ! time level indices
       INTEGER  ::   ji, jj, jk
       REAL(wp) ::   zfactfe,zfactsi,zfactch, zcompadi
@@ -155,5 +158,5 @@
             DO ji = 1, jpi
 
-               zcompadi = MAX( ( trb(ji,jj,jk,jpdia) - 1e-9), 0. )
+               zcompadi = MAX( ( tr(ji,jj,jk,jpdia,Kbb) - 1e-9), 0. )
 
                !    Aggregation term for diatoms is increased in case of nutrient
@@ -165,28 +168,28 @@
                zlim2   = xlimdia(ji,jj,jk) * xlimdia(ji,jj,jk)
                zlim1   = 0.25 * ( 1. - zlim2 ) / ( 0.25 + zlim2 ) 
-               zrespp2 = 1.e6 * xstep * (  wchld + wchldm * zlim1 ) * xdiss(ji,jj,jk) * zcompadi * trb(ji,jj,jk,jpdia)
+               zrespp2 = 1.e6 * xstep * (  wchld + wchldm * zlim1 ) * xdiss(ji,jj,jk) * zcompadi * tr(ji,jj,jk,jpdia,Kbb)
 
                !     Phytoplankton mortality. 
                !     ------------------------
-               ztortp2 = mprat2 * xstep * trb(ji,jj,jk,jpdia)  / ( xkmort + trb(ji,jj,jk,jpdia) ) * zcompadi 
+               ztortp2 = mprat2 * xstep * tr(ji,jj,jk,jpdia,Kbb)  / ( xkmort + tr(ji,jj,jk,jpdia,Kbb) ) * zcompadi 
 
                zmortp2 = zrespp2 + ztortp2
 
-               !   Update the arrays tra which contains the biological sources and sinks
+               !   Update the arrays tr(:,:,:,:,Krhs) which contains the biological sources and sinks
                !   ---------------------------------------------------------------------
-               zfactch = trb(ji,jj,jk,jpdch) / ( trb(ji,jj,jk,jpdia) + rtrn )
-               zfactfe = trb(ji,jj,jk,jpdfe) / ( trb(ji,jj,jk,jpdia) + rtrn )
-               zfactsi = trb(ji,jj,jk,jpdsi) / ( trb(ji,jj,jk,jpdia) + rtrn )
-               tra(ji,jj,jk,jpdia) = tra(ji,jj,jk,jpdia) - zmortp2 
-               tra(ji,jj,jk,jpdch) = tra(ji,jj,jk,jpdch) - zmortp2 * zfactch
-               tra(ji,jj,jk,jpdfe) = tra(ji,jj,jk,jpdfe) - zmortp2 * zfactfe
-               tra(ji,jj,jk,jpdsi) = tra(ji,jj,jk,jpdsi) - zmortp2 * zfactsi
-               tra(ji,jj,jk,jpgsi) = tra(ji,jj,jk,jpgsi) + zmortp2 * zfactsi
-               tra(ji,jj,jk,jpgoc) = tra(ji,jj,jk,jpgoc) + zrespp2 + 0.5 * ztortp2
-               tra(ji,jj,jk,jppoc) = tra(ji,jj,jk,jppoc) + 0.5 * ztortp2
+               zfactch = tr(ji,jj,jk,jpdch,Kbb) / ( tr(ji,jj,jk,jpdia,Kbb) + rtrn )
+               zfactfe = tr(ji,jj,jk,jpdfe,Kbb) / ( tr(ji,jj,jk,jpdia,Kbb) + rtrn )
+               zfactsi = tr(ji,jj,jk,jpdsi,Kbb) / ( tr(ji,jj,jk,jpdia,Kbb) + rtrn )
+               tr(ji,jj,jk,jpdia,Krhs) = tr(ji,jj,jk,jpdia,Krhs) - zmortp2 
+               tr(ji,jj,jk,jpdch,Krhs) = tr(ji,jj,jk,jpdch,Krhs) - zmortp2 * zfactch
+               tr(ji,jj,jk,jpdfe,Krhs) = tr(ji,jj,jk,jpdfe,Krhs) - zmortp2 * zfactfe
+               tr(ji,jj,jk,jpdsi,Krhs) = tr(ji,jj,jk,jpdsi,Krhs) - zmortp2 * zfactsi
+               tr(ji,jj,jk,jpgsi,Krhs) = tr(ji,jj,jk,jpgsi,Krhs) + zmortp2 * zfactsi
+               tr(ji,jj,jk,jpgoc,Krhs) = tr(ji,jj,jk,jpgoc,Krhs) + zrespp2 + 0.5 * ztortp2
+               tr(ji,jj,jk,jppoc,Krhs) = tr(ji,jj,jk,jppoc,Krhs) + 0.5 * ztortp2
                prodpoc(ji,jj,jk) = prodpoc(ji,jj,jk) + 0.5 * ztortp2
                prodgoc(ji,jj,jk) = prodgoc(ji,jj,jk) + zrespp2 + 0.5 * ztortp2
-               tra(ji,jj,jk,jpsfe) = tra(ji,jj,jk,jpsfe) + 0.5 * ztortp2 * zfactfe
-               tra(ji,jj,jk,jpbfe) = tra(ji,jj,jk,jpbfe) + ( zrespp2 + 0.5 * ztortp2 ) * zfactfe
+               tr(ji,jj,jk,jpsfe,Krhs) = tr(ji,jj,jk,jpsfe,Krhs) + 0.5 * ztortp2 * zfactfe
+               tr(ji,jj,jk,jpbfe,Krhs) = tr(ji,jj,jk,jpbfe,Krhs) + ( zrespp2 + 0.5 * ztortp2 ) * zfactfe
             END DO
          END DO
@@ -196,5 +199,5 @@
          WRITE(charout, FMT="('diat')")
          CALL prt_ctl_trc_info(charout)
-         CALL prt_ctl_trc(tab4d=tra, mask=tmask, clinfo=ctrcnm)
+         CALL prt_ctl_trc(tab4d=tr(:,:,:,:,Krhs), mask=tmask, clinfo=ctrcnm)
       ENDIF
       !
