Index: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/TOP/PISCES/P4Z/p4zsms.F90
===================================================================
--- NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/TOP/PISCES/P4Z/p4zsms.F90	(revision 10966)
+++ NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/TOP/PISCES/P4Z/p4zsms.F90	(revision 10975)
@@ -46,5 +46,5 @@
 CONTAINS
 
-   SUBROUTINE p4z_sms( kt, Kbb, Kmm )
+   SUBROUTINE p4z_sms( kt, Kbb, Kmm, Krhs )
       !!---------------------------------------------------------------------
       !!                     ***  ROUTINE p4z_sms  ***
@@ -58,6 +58,6 @@
       !!---------------------------------------------------------------------
       !
-      INTEGER, INTENT( in ) ::   kt           ! ocean time-step index      
-      INTEGER, INTENT( in ) ::   Kbb, Kmm     ! time level index
+      INTEGER, INTENT( in ) ::   kt              ! ocean time-step index      
+      INTEGER, INTENT( in ) ::   Kbb, Kmm, Krhs  ! time level index
       !!
       INTEGER ::   ji, jj, jk, jnt, jn, jl
@@ -73,14 +73,14 @@
         !
         IF( .NOT. ln_rsttr ) THEN
-            CALL p4z_che                              ! initialize the chemical constants
-            CALL ahini_for_at(hi)   !  set PH at kt=nit000
+            CALL p4z_che( Kbb, Kmm )                  ! initialize the chemical constants
+            CALL ahini_for_at( hi, Kbb )              !  set PH at kt=nit000
             t_oce_co2_flx_cum = 0._wp
         ELSE
-            CALL p4z_rst( nittrc000, 'READ' )  !* read or initialize all required fields
+            CALL p4z_rst( nittrc000, Kbb, Kmm,  'READ' )  !* read or initialize all required fields
         ENDIF
         !
       ENDIF
       !
-      IF( ln_pisdmp .AND. MOD( kt - nn_dttrc, nn_pisdmp ) == 0 )   CALL p4z_dmp( kt )      ! Relaxation of some tracers
+      IF( ln_pisdmp .AND. MOD( kt - nn_dttrc, nn_pisdmp ) == 0 )   CALL p4z_dmp( kt, Kbb, Kmm )      ! Relaxation of some tracers
       !
       rfact = r2dttrc
@@ -99,20 +99,20 @@
       IF( ( neuler == 0 .AND. kt == nittrc000 ) .OR. ln_top_euler ) THEN
          DO jn = jp_pcs0, jp_pcs1              !   SMS on tracer without Asselin time-filter
-            trb(:,:,:,jn) = trn(:,:,:,jn)
+            tr(:,:,:,jn,Kbb) = tr(:,:,:,jn,Kmm)
          END DO
       ENDIF
       !
-      IF( ll_sbc ) CALL p4z_sbc( kt )   ! external sources of nutrients 
+      IF( ll_sbc ) CALL p4z_sbc( kt, Kmm )   ! external sources of nutrients 
       !
 #if ! defined key_sed_off
-      CALL p4z_che              ! computation of chemical constants
-      CALL p4z_int( kt )        ! computation of various rates for biogeochemistry
+      CALL p4z_che(     Kbb, Kmm       ) ! computation of chemical constants
+      CALL p4z_int( kt, Kbb, Kmm       ) ! computation of various rates for biogeochemistry
       !
       DO jnt = 1, nrdttrc          ! Potential time splitting if requested
          !
-         CALL p4z_bio( kt, jnt, Kbb, Kmm )   ! Biology
-         CALL p4z_lys( kt, jnt )   ! Compute CaCO3 saturation
-         CALL p4z_sed( kt, jnt )   ! Surface and Bottom boundary conditions
-         CALL p4z_flx( kt, jnt )   ! Compute surface fluxes
+         CALL p4z_bio( kt, jnt, Kbb, Kmm, Krhs )   ! Biology
+         CALL p4z_lys( kt, jnt, Kbb,      Krhs )   ! Compute CaCO3 saturation
+         CALL p4z_sed( kt, jnt, Kbb, Kmm, Krhs )   ! Surface and Bottom boundary conditions
+         CALL p4z_flx( kt, jnt, Kbb, Kmm, Krhs )   ! Compute surface fluxes
          !
          xnegtr(:,:,:) = 1.e0
@@ -121,6 +121,6 @@
                DO jj = 1, jpj
                   DO ji = 1, jpi
-                     IF( ( trb(ji,jj,jk,jn) + tra(ji,jj,jk,jn) ) < 0.e0 ) THEN
-                        ztra             = ABS( trb(ji,jj,jk,jn) ) / ( ABS( tra(ji,jj,jk,jn) ) + rtrn )
+                     IF( ( tr(ji,jj,jk,jn,Kbb) + tr(ji,jj,jk,jn,Krhs) ) < 0.e0 ) THEN
+                        ztra             = ABS( tr(ji,jj,jk,jn,Kbb) ) / ( ABS( tr(ji,jj,jk,jn,Krhs) ) + rtrn )
                         xnegtr(ji,jj,jk) = MIN( xnegtr(ji,jj,jk),  ztra )
                      ENDIF
@@ -132,14 +132,14 @@
          !                                ! 
          DO jn = jp_pcs0, jp_pcs1
-           trb(:,:,:,jn) = trb(:,:,:,jn) + xnegtr(:,:,:) * tra(:,:,:,jn)
+           tr(:,:,:,jn,Kbb) = tr(:,:,:,jn,Kbb) + xnegtr(:,:,:) * tr(:,:,:,jn,Krhs)
          END DO
         !
          DO jn = jp_pcs0, jp_pcs1
-            tra(:,:,:,jn) = 0._wp
+            tr(:,:,:,jn,Krhs) = 0._wp
          END DO
          !
          IF( ln_top_euler ) THEN
             DO jn = jp_pcs0, jp_pcs1
-               trn(:,:,:,jn) = trb(:,:,:,jn)
+               tr(:,:,:,jn,Kmm) = tr(:,:,:,jn,Kbb)
             END DO
          ENDIF
@@ -149,5 +149,5 @@
       IF( l_trdtrc ) THEN
          DO jn = jp_pcs0, jp_pcs1
-           CALL trd_trc( tra(:,:,:,jn), jn, jptra_sms, kt, Kmm )   ! save trends
+           CALL trd_trc( tr(:,:,:,jn,Krhs), jn, jptra_sms, kt, Kmm )   ! save trends
          END DO
       END IF
@@ -156,9 +156,9 @@
       IF( ln_sediment ) THEN 
          !
-         CALL sed_model( kt, Kmm )     !  Main program of Sediment model
+         CALL sed_model( kt, Kbb, Kmm, Krhs )     !  Main program of Sediment model
          !
          IF( ln_top_euler ) THEN
             DO jn = jp_pcs0, jp_pcs1
-               trn(:,:,:,jn) = trb(:,:,:,jn)
+               tr(:,:,:,jn,Kmm) = tr(:,:,:,jn,Kbb)
             END DO
          ENDIF
@@ -166,10 +166,10 @@
       ENDIF
       !
-      IF( lrst_trc )  CALL p4z_rst( kt, 'WRITE' )  !* Write PISCES informations in restart file 
-      !
-
-      IF( lk_iomput .OR. ln_check_mass )  CALL p4z_chk_mass( kt )    ! Mass conservation checking
-
-      IF( lwm .AND. kt == nittrc000    )  CALL FLUSH( numonp )       ! flush output namelist PISCES
+      IF( lrst_trc )  CALL p4z_rst( kt, Kbb, Kmm,  'WRITE' )           !* Write PISCES informations in restart file 
+      !
+
+      IF( lk_iomput .OR. ln_check_mass )  CALL p4z_chk_mass( kt, Kmm ) ! Mass conservation checking
+
+      IF( lwm .AND. kt == nittrc000    )  CALL FLUSH( numonp )         ! flush output namelist PISCES
       !
       IF( ln_timing )  CALL timing_stop('p4z_sms')
@@ -265,5 +265,5 @@
 
 
-   SUBROUTINE p4z_rst( kt, cdrw )
+   SUBROUTINE p4z_rst( kt, Kbb, Kmm, cdrw )
       !!---------------------------------------------------------------------
       !!                   ***  ROUTINE p4z_rst  ***
@@ -276,4 +276,5 @@
       !!---------------------------------------------------------------------
       INTEGER         , INTENT(in) ::   kt         ! ocean time-step
+      INTEGER         , INTENT(in) ::   Kbb, Kmm   ! time level indices
       CHARACTER(len=*), INTENT(in) ::   cdrw       ! "READ"/"WRITE" flag
       !!---------------------------------------------------------------------
@@ -288,6 +289,6 @@
             CALL iom_get( numrtr, jpdom_autoglo, 'PH' , hi(:,:,:)  )
          ELSE
-            CALL p4z_che                              ! initialize the chemical constants
-            CALL ahini_for_at(hi)
+            CALL p4z_che( Kbb, Kmm )                  ! initialize the chemical constants
+            CALL ahini_for_at( hi, Kbb )
          ENDIF
          CALL iom_get( numrtr, jpdom_autoglo, 'Silicalim', xksi(:,:) )
@@ -336,5 +337,5 @@
 
 
-   SUBROUTINE p4z_dmp( kt )
+   SUBROUTINE p4z_dmp( kt, Kbb, Kmm )
       !!----------------------------------------------------------------------
       !!                    ***  p4z_dmp  ***
@@ -343,5 +344,6 @@
       !!----------------------------------------------------------------------
       !
-      INTEGER, INTENT( in )  ::     kt ! time step
+      INTEGER, INTENT( in )  ::     kt            ! time step
+      INTEGER, INTENT( in )  ::     Kbb, Kmm      ! time level indices
       !
       REAL(wp) ::  alkmean = 2426.     ! mean value of alkalinity ( Glodap ; for Goyet 2391. )
@@ -364,40 +366,40 @@
             zarea          = 1._wp / glob_sum( 'p4zsms', cvol(:,:,:) ) * 1e6              
 
-            zalksumn = glob_sum( 'p4zsms', trn(:,:,:,jptal) * cvol(:,:,:)  ) * zarea
-            zpo4sumn = glob_sum( 'p4zsms', trn(:,:,:,jppo4) * cvol(:,:,:)  ) * zarea * po4r
-            zno3sumn = glob_sum( 'p4zsms', trn(:,:,:,jpno3) * cvol(:,:,:)  ) * zarea * rno3
-            zsilsumn = glob_sum( 'p4zsms', trn(:,:,:,jpsil) * cvol(:,:,:)  ) * zarea
+            zalksumn = glob_sum( 'p4zsms', tr(:,:,:,jptal,Kmm) * cvol(:,:,:)  ) * zarea
+            zpo4sumn = glob_sum( 'p4zsms', tr(:,:,:,jppo4,Kmm) * cvol(:,:,:)  ) * zarea * po4r
+            zno3sumn = glob_sum( 'p4zsms', tr(:,:,:,jpno3,Kmm) * cvol(:,:,:)  ) * zarea * rno3
+            zsilsumn = glob_sum( 'p4zsms', tr(:,:,:,jpsil,Kmm) * cvol(:,:,:)  ) * zarea
  
             IF(lwp) WRITE(numout,*) '       TALKN mean : ', zalksumn
-            trn(:,:,:,jptal) = trn(:,:,:,jptal) * alkmean / zalksumn
+            tr(:,:,:,jptal,Kmm) = tr(:,:,:,jptal,Kmm) * alkmean / zalksumn
 
             IF(lwp) WRITE(numout,*) '       PO4N  mean : ', zpo4sumn
-            trn(:,:,:,jppo4) = trn(:,:,:,jppo4) * po4mean / zpo4sumn
+            tr(:,:,:,jppo4,Kmm) = tr(:,:,:,jppo4,Kmm) * po4mean / zpo4sumn
 
             IF(lwp) WRITE(numout,*) '       NO3N  mean : ', zno3sumn
-            trn(:,:,:,jpno3) = trn(:,:,:,jpno3) * no3mean / zno3sumn
+            tr(:,:,:,jpno3,Kmm) = tr(:,:,:,jpno3,Kmm) * no3mean / zno3sumn
 
             IF(lwp) WRITE(numout,*) '       SiO3N mean : ', zsilsumn
-            trn(:,:,:,jpsil) = MIN( 400.e-6,trn(:,:,:,jpsil) * silmean / zsilsumn )
+            tr(:,:,:,jpsil,Kmm) = MIN( 400.e-6,tr(:,:,:,jpsil,Kmm) * silmean / zsilsumn )
             !
             !
             IF( .NOT. ln_top_euler ) THEN
-               zalksumb = glob_sum( 'p4zsms', trb(:,:,:,jptal) * cvol(:,:,:)  ) * zarea
-               zpo4sumb = glob_sum( 'p4zsms', trb(:,:,:,jppo4) * cvol(:,:,:)  ) * zarea * po4r
-               zno3sumb = glob_sum( 'p4zsms', trb(:,:,:,jpno3) * cvol(:,:,:)  ) * zarea * rno3
-               zsilsumb = glob_sum( 'p4zsms', trb(:,:,:,jpsil) * cvol(:,:,:)  ) * zarea
+               zalksumb = glob_sum( 'p4zsms', tr(:,:,:,jptal,Kbb) * cvol(:,:,:)  ) * zarea
+               zpo4sumb = glob_sum( 'p4zsms', tr(:,:,:,jppo4,Kbb) * cvol(:,:,:)  ) * zarea * po4r
+               zno3sumb = glob_sum( 'p4zsms', tr(:,:,:,jpno3,Kbb) * cvol(:,:,:)  ) * zarea * rno3
+               zsilsumb = glob_sum( 'p4zsms', tr(:,:,:,jpsil,Kbb) * cvol(:,:,:)  ) * zarea
  
                IF(lwp) WRITE(numout,*) ' '
                IF(lwp) WRITE(numout,*) '       TALKB mean : ', zalksumb
-               trb(:,:,:,jptal) = trb(:,:,:,jptal) * alkmean / zalksumb
+               tr(:,:,:,jptal,Kbb) = tr(:,:,:,jptal,Kbb) * alkmean / zalksumb
 
                IF(lwp) WRITE(numout,*) '       PO4B  mean : ', zpo4sumb
-               trb(:,:,:,jppo4) = trb(:,:,:,jppo4) * po4mean / zpo4sumb
+               tr(:,:,:,jppo4,Kbb) = tr(:,:,:,jppo4,Kbb) * po4mean / zpo4sumb
 
                IF(lwp) WRITE(numout,*) '       NO3B  mean : ', zno3sumb
-               trb(:,:,:,jpno3) = trb(:,:,:,jpno3) * no3mean / zno3sumb
+               tr(:,:,:,jpno3,Kbb) = tr(:,:,:,jpno3,Kbb) * no3mean / zno3sumb
 
                IF(lwp) WRITE(numout,*) '       SiO3B mean : ', zsilsumb
-               trb(:,:,:,jpsil) = MIN( 400.e-6,trb(:,:,:,jpsil) * silmean / zsilsumb )
+               tr(:,:,:,jpsil,Kbb) = MIN( 400.e-6,tr(:,:,:,jpsil,Kbb) * silmean / zsilsumb )
            ENDIF
         ENDIF
@@ -408,5 +410,5 @@
 
 
-   SUBROUTINE p4z_chk_mass( kt )
+   SUBROUTINE p4z_chk_mass( kt, Kmm )
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE p4z_chk_mass  ***
@@ -416,4 +418,5 @@
       !!---------------------------------------------------------------------
       INTEGER, INTENT( in ) ::   kt      ! ocean time-step index      
+      INTEGER, INTENT( in ) ::   Kmm     ! time level indices
       REAL(wp)             ::  zrdenittot, zsdenittot, znitrpottot
       CHARACTER(LEN=100)   ::   cltxt
@@ -439,13 +442,13 @@
          !   Compute the budget of NO3, ALK, Si, Fer
          IF( ln_p4z ) THEN
-            zwork(:,:,:) =    trn(:,:,:,jpno3) + trn(:,:,:,jpnh4)                      &
-               &          +   trn(:,:,:,jpphy) + trn(:,:,:,jpdia)                      &
-               &          +   trn(:,:,:,jppoc) + trn(:,:,:,jpgoc)  + trn(:,:,:,jpdoc)  &        
-               &          +   trn(:,:,:,jpzoo) + trn(:,:,:,jpmes) 
+            zwork(:,:,:) =    tr(:,:,:,jpno3,Kmm) + tr(:,:,:,jpnh4,Kmm)                      &
+               &          +   tr(:,:,:,jpphy,Kmm) + tr(:,:,:,jpdia,Kmm)                      &
+               &          +   tr(:,:,:,jppoc,Kmm) + tr(:,:,:,jpgoc,Kmm)  + tr(:,:,:,jpdoc,Kmm)  &        
+               &          +   tr(:,:,:,jpzoo,Kmm) + tr(:,:,:,jpmes,Kmm) 
         ELSE
-            zwork(:,:,:) =    trn(:,:,:,jpno3) + trn(:,:,:,jpnh4) + trn(:,:,:,jpnph)   &
-               &          +   trn(:,:,:,jpndi) + trn(:,:,:,jpnpi)                      & 
-               &          +   trn(:,:,:,jppon) + trn(:,:,:,jpgon) + trn(:,:,:,jpdon)   &
-               &          + ( trn(:,:,:,jpzoo) + trn(:,:,:,jpmes) ) * no3rat3 
+            zwork(:,:,:) =    tr(:,:,:,jpno3,Kmm) + tr(:,:,:,jpnh4,Kmm) + tr(:,:,:,jpnph,Kmm)   &
+               &          +   tr(:,:,:,jpndi,Kmm) + tr(:,:,:,jpnpi,Kmm)                      & 
+               &          +   tr(:,:,:,jppon,Kmm) + tr(:,:,:,jpgon,Kmm) + tr(:,:,:,jpdon,Kmm)   &
+               &          + ( tr(:,:,:,jpzoo,Kmm) + tr(:,:,:,jpmes,Kmm) ) * no3rat3 
         ENDIF
         !
@@ -457,13 +460,13 @@
       IF( iom_use( "ppo4tot" ) .OR. ( ln_check_mass .AND. kt == nitend )  ) THEN
          IF( ln_p4z ) THEN
-            zwork(:,:,:) =    trn(:,:,:,jppo4)                                         &
-               &          +   trn(:,:,:,jpphy) + trn(:,:,:,jpdia)                      &
-               &          +   trn(:,:,:,jppoc) + trn(:,:,:,jpgoc)  + trn(:,:,:,jpdoc)  &        
-               &          +   trn(:,:,:,jpzoo) + trn(:,:,:,jpmes) 
+            zwork(:,:,:) =    tr(:,:,:,jppo4,Kmm)                                         &
+               &          +   tr(:,:,:,jpphy,Kmm) + tr(:,:,:,jpdia,Kmm)                      &
+               &          +   tr(:,:,:,jppoc,Kmm) + tr(:,:,:,jpgoc,Kmm)  + tr(:,:,:,jpdoc,Kmm)  &        
+               &          +   tr(:,:,:,jpzoo,Kmm) + tr(:,:,:,jpmes,Kmm) 
         ELSE
-            zwork(:,:,:) =    trn(:,:,:,jppo4) + trn(:,:,:,jppph)                      &
-               &          +   trn(:,:,:,jppdi) + trn(:,:,:,jpppi)                      & 
-               &          +   trn(:,:,:,jppop) + trn(:,:,:,jpgop) + trn(:,:,:,jpdop)   &
-               &          + ( trn(:,:,:,jpzoo) + trn(:,:,:,jpmes) ) * po4rat3 
+            zwork(:,:,:) =    tr(:,:,:,jppo4,Kmm) + tr(:,:,:,jppph,Kmm)                      &
+               &          +   tr(:,:,:,jppdi,Kmm) + tr(:,:,:,jpppi,Kmm)                      & 
+               &          +   tr(:,:,:,jppop,Kmm) + tr(:,:,:,jpgop,Kmm) + tr(:,:,:,jpdop,Kmm)   &
+               &          + ( tr(:,:,:,jpzoo,Kmm) + tr(:,:,:,jpmes,Kmm) ) * po4rat3 
         ENDIF
         !
@@ -474,5 +477,5 @@
       !
       IF( iom_use( "psiltot" ) .OR. ( ln_check_mass .AND. kt == nitend )  ) THEN
-         zwork(:,:,:) =  trn(:,:,:,jpsil) + trn(:,:,:,jpgsi) + trn(:,:,:,jpdsi) 
+         zwork(:,:,:) =  tr(:,:,:,jpsil,Kmm) + tr(:,:,:,jpgsi,Kmm) + tr(:,:,:,jpdsi,Kmm) 
          !
          silbudget = glob_sum( 'p4zsms', zwork(:,:,:) * cvol(:,:,:)  )  
@@ -482,5 +485,5 @@
       !
       IF( iom_use( "palktot" ) .OR. ( ln_check_mass .AND. kt == nitend )  ) THEN
-         zwork(:,:,:) =  trn(:,:,:,jpno3) * rno3 + trn(:,:,:,jptal) + trn(:,:,:,jpcal) * 2.              
+         zwork(:,:,:) =  tr(:,:,:,jpno3,Kmm) * rno3 + tr(:,:,:,jptal,Kmm) + tr(:,:,:,jpcal,Kmm) * 2.              
          !
          alkbudget = glob_sum( 'p4zsms', zwork(:,:,:) * cvol(:,:,:)  )         !
@@ -490,7 +493,7 @@
       !
       IF( iom_use( "pfertot" ) .OR. ( ln_check_mass .AND. kt == nitend )  ) THEN
-         zwork(:,:,:) =   trn(:,:,:,jpfer) + trn(:,:,:,jpnfe) + trn(:,:,:,jpdfe)   &
-            &         +   trn(:,:,:,jpbfe) + trn(:,:,:,jpsfe)                      &
-            &         + ( trn(:,:,:,jpzoo) + trn(:,:,:,jpmes) )  * ferat3    
+         zwork(:,:,:) =   tr(:,:,:,jpfer,Kmm) + tr(:,:,:,jpnfe,Kmm) + tr(:,:,:,jpdfe,Kmm)   &
+            &         +   tr(:,:,:,jpbfe,Kmm) + tr(:,:,:,jpsfe,Kmm)                      &
+            &         + ( tr(:,:,:,jpzoo,Kmm) + tr(:,:,:,jpmes,Kmm) )  * ferat3    
          !
          ferbudget = glob_sum( 'p4zsms', zwork(:,:,:) * cvol(:,:,:)  )  
