Index: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/TOP/PISCES/P4Z/p4zche.F90
===================================================================
--- NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/TOP/PISCES/P4Z/p4zche.F90	(revision 10740)
+++ NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/TOP/PISCES/P4Z/p4zche.F90	(revision 10975)
@@ -137,5 +137,5 @@
 CONTAINS
 
-   SUBROUTINE p4z_che
+   SUBROUTINE p4z_che( Kbb, Kmm )
       !!---------------------------------------------------------------------
       !!                     ***  ROUTINE p4z_che  ***
@@ -145,4 +145,5 @@
       !! ** Method  : - ...
       !!---------------------------------------------------------------------
+      INTEGER, INTENT(in) ::   Kbb, Kmm  ! time level indices
       INTEGER  ::   ji, jj, jk
       REAL(wp) ::   ztkel, ztkel1, zt , zsal  , zsal2 , zbuf1 , zbuf2
@@ -164,7 +165,7 @@
       ! -------------------------------------------------------------
       IF (neos == -1) THEN
-         salinprac(:,:,:) = tsn(:,:,:,jp_sal) * 35.0 / 35.16504
+         salinprac(:,:,:) = ts(:,:,:,jp_sal,Kmm) * 35.0 / 35.16504
       ELSE
-         salinprac(:,:,:) = tsn(:,:,:,jp_sal)
+         salinprac(:,:,:) = ts(:,:,:,jp_sal,Kmm)
       ENDIF
 
@@ -178,8 +179,8 @@
          DO jj = 1, jpj
             DO ji = 1, jpi
-               zpres = gdept_n(ji,jj,jk) / 1000.
-               za1 = 0.04 * ( 1.0 + 0.185 * tsn(ji,jj,jk,jp_tem) + 0.035 * (salinprac(ji,jj,jk) - 35.0) )
-               za2 = 0.0075 * ( 1.0 - tsn(ji,jj,jk,jp_tem) / 30.0 )
-               tempis(ji,jj,jk) = tsn(ji,jj,jk,jp_tem) - za1 * zpres + za2 * zpres**2
+               zpres = gdept(ji,jj,jk,Kmm) / 1000.
+               za1 = 0.04 * ( 1.0 + 0.185 * ts(ji,jj,jk,jp_tem,Kmm) + 0.035 * (salinprac(ji,jj,jk) - 35.0) )
+               za2 = 0.0075 * ( 1.0 - ts(ji,jj,jk,jp_tem,Kmm) / 30.0 )
+               tempis(ji,jj,jk) = ts(ji,jj,jk,jp_tem,Kmm) - za1 * zpres + za2 * zpres**2
             END DO
          END DO
@@ -245,5 +246,5 @@
                zplat   = SIN ( ABS(gphit(ji,jj)*3.141592654/180.) )
                zc1 = 5.92E-3 + zplat**2 * 5.25E-3
-               zpres = ((1-zc1)-SQRT(((1-zc1)**2)-(8.84E-6*gdept_n(ji,jj,jk)))) / 4.42E-6
+               zpres = ((1-zc1)-SQRT(((1-zc1)**2)-(8.84E-6*gdept(ji,jj,jk,Kmm)))) / 4.42E-6
                zpres = zpres / 10.0
 
@@ -448,5 +449,5 @@
    END SUBROUTINE p4z_che
 
-   SUBROUTINE ahini_for_at(p_hini)
+   SUBROUTINE ahini_for_at(p_hini, Kbb )
       !!---------------------------------------------------------------------
       !!                     ***  ROUTINE ahini_for_at  ***
@@ -462,4 +463,5 @@
       !!---------------------------------------------------------------------
       REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(OUT)  ::  p_hini
+      INTEGER,                          INTENT(in)   ::  Kbb      ! time level indices
       INTEGER  ::   ji, jj, jk
       REAL(wp)  ::  zca1, zba1
@@ -474,6 +476,6 @@
         DO jj = 1, jpj
           DO ji = 1, jpi
-            p_alkcb  = trb(ji,jj,jk,jptal) * 1000. / (rhop(ji,jj,jk) + rtrn)
-            p_dictot = trb(ji,jj,jk,jpdic) * 1000. / (rhop(ji,jj,jk) + rtrn)
+            p_alkcb  = tr(ji,jj,jk,jptal,Kbb) * 1000. / (rhop(ji,jj,jk) + rtrn)
+            p_dictot = tr(ji,jj,jk,jpdic,Kbb) * 1000. / (rhop(ji,jj,jk) + rtrn)
             p_bortot = borat(ji,jj,jk)
             IF (p_alkcb <= 0.) THEN
@@ -516,5 +518,5 @@
    !===============================================================================
 
-   SUBROUTINE anw_infsup( p_alknw_inf, p_alknw_sup )
+   SUBROUTINE anw_infsup( p_alknw_inf, p_alknw_sup, Kbb )
 
    ! Subroutine returns the lower and upper bounds of "non-water-selfionization"
@@ -525,8 +527,9 @@
    REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(OUT) :: p_alknw_inf
    REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(OUT) :: p_alknw_sup
-
-   p_alknw_inf(:,:,:) =  -trb(:,:,:,jppo4) * 1000. / (rhop(:,:,:) + rtrn) - sulfat(:,:,:)  &
+   INTEGER,                          INTENT(in)  ::  Kbb      ! time level indices
+
+   p_alknw_inf(:,:,:) =  -tr(:,:,:,jppo4,Kbb) * 1000. / (rhop(:,:,:) + rtrn) - sulfat(:,:,:)  &
    &              - fluorid(:,:,:)
-   p_alknw_sup(:,:,:) =   (2. * trb(:,:,:,jpdic) + 2. * trb(:,:,:,jppo4) + trb(:,:,:,jpsil) )    &
+   p_alknw_sup(:,:,:) =   (2. * tr(:,:,:,jpdic,Kbb) + 2. * tr(:,:,:,jppo4,Kbb) + tr(:,:,:,jpsil,Kbb) )    &
    &               * 1000. / (rhop(:,:,:) + rtrn) + borat(:,:,:) 
 
@@ -534,5 +537,5 @@
 
 
-   SUBROUTINE solve_at_general( p_hini, zhi )
+   SUBROUTINE solve_at_general( p_hini, zhi, Kbb )
 
    ! Universal pH solver that converges from any given initial value,
@@ -543,4 +546,5 @@
    REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(IN)   :: p_hini
    REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(OUT)  :: zhi
+   INTEGER,                          INTENT(in)   :: Kbb  ! time level indices
 
    ! Local variables
@@ -565,5 +569,5 @@
    IF( ln_timing )  CALL timing_start('solve_at_general')
 
-   CALL anw_infsup( zalknw_inf, zalknw_sup )
+   CALL anw_infsup( zalknw_inf, zalknw_sup, Kbb )
 
    rmask(:,:,:) = tmask(:,:,:)
@@ -575,5 +579,5 @@
          DO ji = 1, jpi
             IF (rmask(ji,jj,jk) == 1.) THEN
-               p_alktot = trb(ji,jj,jk,jptal) * 1000. / (rhop(ji,jj,jk) + rtrn)
+               p_alktot = tr(ji,jj,jk,jptal,Kbb) * 1000. / (rhop(ji,jj,jk) + rtrn)
                aphscale = 1. + sulfat(ji,jj,jk)/aks3(ji,jj,jk)
                zh_ini = p_hini(ji,jj,jk)
@@ -609,9 +613,9 @@
             IF (rmask(ji,jj,jk) == 1.) THEN
                zfact = rhop(ji,jj,jk) / 1000. + rtrn
-               p_alktot = trb(ji,jj,jk,jptal) / zfact
-               zdic  = trb(ji,jj,jk,jpdic) / zfact
+               p_alktot = tr(ji,jj,jk,jptal,Kbb) / zfact
+               zdic  = tr(ji,jj,jk,jpdic,Kbb) / zfact
                zbot  = borat(ji,jj,jk)
-               zpt = trb(ji,jj,jk,jppo4) / zfact * po4r
-               zsit = trb(ji,jj,jk,jpsil) / zfact
+               zpt = tr(ji,jj,jk,jppo4,Kbb) / zfact * po4r
+               zsit = tr(ji,jj,jk,jpsil,Kbb) / zfact
                zst = sulfat (ji,jj,jk)
                zft = fluorid(ji,jj,jk)
