Index: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DOM/istate.F90
===================================================================
--- NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DOM/istate.F90	(revision 10922)
+++ NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DOM/istate.F90	(revision 10928)
@@ -51,5 +51,5 @@
 CONTAINS
 
-   SUBROUTINE istate_init( Kbb, Kmm )
+   SUBROUTINE istate_init( Kbb, Kmm, Kaa )
       !!----------------------------------------------------------------------
       !!                   ***  ROUTINE istate_init  ***
@@ -57,5 +57,6 @@
       !! ** Purpose :   Initialization of the dynamics and tracer fields.
       !!----------------------------------------------------------------------
-      INTEGER, INTENT(in) ::   Kbb, Kmm   ! ocean time level indices
+      INTEGER, INTENT( in )  ::  Kbb, Kmm, Kaa   ! ocean time level indices
+      !
       INTEGER ::   ji, jj, jk   ! dummy loop indices
 !!gm see comment further down
@@ -77,9 +78,9 @@
       rhd  (:,:,:  ) = 0._wp   ;   rhop (:,:,:  ) = 0._wp      ! set one for all to 0 at level jpk
       rn2b (:,:,:  ) = 0._wp   ;   rn2  (:,:,:  ) = 0._wp      ! set one for all to 0 at levels 1 and jpk
-      tsa  (:,:,:,:) = 0._wp                                   ! set one for all to 0 at level jpk
+      ts  (:,:,:,:,Kaa) = 0._wp                                   ! set one for all to 0 at level jpk
       rab_b(:,:,:,:) = 0._wp   ;   rab_n(:,:,:,:) = 0._wp      ! set one for all to 0 at level jpk
 #if defined key_agrif
-      ua   (:,:,:  ) = 0._wp   ! used in agrif_oce_sponge at initialization
-      va   (:,:,:  ) = 0._wp   ! used in agrif_oce_sponge at initialization    
+      uu   (:,:,:  ,Kaa) = 0._wp   ! used in agrif_oce_sponge at initialization
+      vv   (:,:,:  ,Kaa) = 0._wp   ! used in agrif_oce_sponge at initialization    
 #endif
 
@@ -98,9 +99,9 @@
          !
          IF( ln_tsd_init ) THEN               
-            CALL dta_tsd( nit000, tsb )       ! read 3D T and S data at nit000
+            CALL dta_tsd( nit000, ts(:,:,:,:,Kbb) )       ! read 3D T and S data at nit000
             !
-            sshb(:,:)   = 0._wp               ! set the ocean at rest
+            ssh(:,:,Kbb)   = 0._wp               ! set the ocean at rest
             IF( ll_wd ) THEN
-               sshb(:,:) =  -ssh_ref  ! Added in 30 here for bathy that adds 30 as Iterative test CEOD 
+               ssh(:,:,Kbb) =  -ssh_ref  ! Added in 30 here for bathy that adds 30 as Iterative test CEOD 
                !
                ! Apply minimum wetdepth criterion
@@ -108,26 +109,26 @@
                DO jj = 1,jpj
                   DO ji = 1,jpi
-                     IF( ht_0(ji,jj) + sshb(ji,jj)  < rn_wdmin1 ) THEN
-                        sshb(ji,jj) = tmask(ji,jj,1)*( rn_wdmin1 - (ht_0(ji,jj)) )
+                     IF( ht_0(ji,jj) + ssh(ji,jj,Kbb)  < rn_wdmin1 ) THEN
+                        ssh(ji,jj,Kbb) = tmask(ji,jj,1)*( rn_wdmin1 - (ht_0(ji,jj)) )
                      ENDIF
                   END DO
                END DO 
             ENDIF 
-            ub  (:,:,:) = 0._wp
-            vb  (:,:,:) = 0._wp  
+            uu  (:,:,:,Kbb) = 0._wp
+            vv  (:,:,:,Kbb) = 0._wp  
             !
          ELSE                                 ! user defined initial T and S
-            CALL usr_def_istate( gdept_b, tmask, tsb, ub, vb, sshb  )         
+            CALL usr_def_istate( gdept(:,:,:,Kbb), tmask, ts(:,:,:,:,Kbb), uu(:,:,:,Kbb), vv(:,:,:,Kbb), ssh(:,:,Kbb)  )         
          ENDIF
-         tsn  (:,:,:,:) = tsb (:,:,:,:)       ! set now values from to before ones
-         sshn (:,:)     = sshb(:,:)   
-         un   (:,:,:)   = ub  (:,:,:)
-         vn   (:,:,:)   = vb  (:,:,:)
-         hdivn(:,:,jpk) = 0._wp               ! bottom divergence set one for 0 to zero at jpk level
-         CALL div_hor( 0, Kmm )               ! compute interior hdivn value  
-!!gm                                    hdivn(:,:,:) = 0._wp
+         ts  (:,:,:,:,Kmm) = ts (:,:,:,:,Kbb)       ! set now values from to before ones
+         ssh (:,:,Kmm)     = ssh(:,:,Kbb)   
+         uu   (:,:,:,Kmm)   = uu  (:,:,:,Kbb)
+         vv   (:,:,:,Kmm)   = vv  (:,:,:,Kbb)
+         hdiv(:,:,jpk) = 0._wp               ! bottom divergence set one for 0 to zero at jpk level
+         CALL div_hor( 0, Kmm )              ! compute interior hdiv value  
+!!gm                                    hdiv(:,:,:) = 0._wp
 
 !!gm POTENTIAL BUG :
-!!gm  ISSUE :  if sshb /= 0  then, in non linear free surface, the e3._n, e3._b should be recomputed
+!!gm  ISSUE :  if ssh(:,:,Kbb) /= 0  then, in non linear free surface, the e3._n, e3._b should be recomputed
 !!             as well as gdept and gdepw....   !!!!! 
 !!      ===>>>>   probably a call to domvvl initialisation here....
@@ -139,14 +140,14 @@
 !            ALLOCATE( zuvd(jpi,jpj,jpk,2) )
 !            CALL dta_uvd( nit000, zuvd )
-!            ub(:,:,:) = zuvd(:,:,:,1) ;  un(:,:,:) = ub(:,:,:)
-!            vb(:,:,:) = zuvd(:,:,:,2) ;  vn(:,:,:) = vb(:,:,:)
+!            uu(:,:,:,Kbb) = zuvd(:,:,:,1) ;  uu(:,:,:,Kmm) = uu(:,:,:,Kbb)
+!            vv(:,:,:,Kbb) = zuvd(:,:,:,2) ;  vv(:,:,:,Kmm) = vv(:,:,:,Kbb)
 !            DEALLOCATE( zuvd )
 !         ENDIF
          !
 !!gm This is to be changed !!!!
-!         ! - ML - sshn could be modified by istate_eel, so that initialization of e3t_b is done here
+!         ! - ML - ssh(:,:,Kmm) could be modified by istate_eel, so that initialization of e3t(:,:,:,Kbb) is done here
 !         IF( .NOT.ln_linssh ) THEN
 !            DO jk = 1, jpk
-!               e3t_b(:,:,jk) = e3t_n(:,:,jk)
+!               e3t(:,:,jk,Kbb) = e3t(:,:,jk,Kmm)
 !            END DO
 !         ENDIF
@@ -158,25 +159,25 @@
       ! Do it whatever the free surface method, these arrays being eventually used
       !
-      un_b(:,:) = 0._wp   ;   vn_b(:,:) = 0._wp
-      ub_b(:,:) = 0._wp   ;   vb_b(:,:) = 0._wp
+      uu_b(:,:,Kmm) = 0._wp   ;   vv_b(:,:,Kmm) = 0._wp
+      uu_b(:,:,Kbb) = 0._wp   ;   vv_b(:,:,Kbb) = 0._wp
       !
-!!gm  the use of umsak & vmask is not necessary below as un, vn, ub, vb are always masked
+!!gm  the use of umsak & vmask is not necessary below as uu(:,:,:,Kmm), vv(:,:,:,Kmm), uu(:,:,:,Kbb), vv(:,:,:,Kbb) are always masked
       DO jk = 1, jpkm1
          DO jj = 1, jpj
             DO ji = 1, jpi
-               un_b(ji,jj) = un_b(ji,jj) + e3u_n(ji,jj,jk) * un(ji,jj,jk) * umask(ji,jj,jk)
-               vn_b(ji,jj) = vn_b(ji,jj) + e3v_n(ji,jj,jk) * vn(ji,jj,jk) * vmask(ji,jj,jk)
+               uu_b(ji,jj,Kmm) = uu_b(ji,jj,Kmm) + e3u(ji,jj,jk,Kmm) * uu(ji,jj,jk,Kmm) * umask(ji,jj,jk)
+               vv_b(ji,jj,Kmm) = vv_b(ji,jj,Kmm) + e3v(ji,jj,jk,Kmm) * vv(ji,jj,jk,Kmm) * vmask(ji,jj,jk)
                !
-               ub_b(ji,jj) = ub_b(ji,jj) + e3u_b(ji,jj,jk) * ub(ji,jj,jk) * umask(ji,jj,jk)
-               vb_b(ji,jj) = vb_b(ji,jj) + e3v_b(ji,jj,jk) * vb(ji,jj,jk) * vmask(ji,jj,jk)
+               uu_b(ji,jj,Kbb) = uu_b(ji,jj,Kbb) + e3u(ji,jj,jk,Kbb) * uu(ji,jj,jk,Kbb) * umask(ji,jj,jk)
+               vv_b(ji,jj,Kbb) = vv_b(ji,jj,Kbb) + e3v(ji,jj,jk,Kbb) * vv(ji,jj,jk,Kbb) * vmask(ji,jj,jk)
             END DO
          END DO
       END DO
       !
-      un_b(:,:) = un_b(:,:) * r1_hu_n(:,:)
-      vn_b(:,:) = vn_b(:,:) * r1_hv_n(:,:)
+      uu_b(:,:,Kmm) = uu_b(:,:,Kmm) * r1_hu_n(:,:)
+      vv_b(:,:,Kmm) = vv_b(:,:,Kmm) * r1_hv_n(:,:)
       !
-      ub_b(:,:) = ub_b(:,:) * r1_hu_b(:,:)
-      vb_b(:,:) = vb_b(:,:) * r1_hv_b(:,:)
+      uu_b(:,:,Kbb) = uu_b(:,:,Kbb) * r1_hu_b(:,:)
+      vv_b(:,:,Kbb) = vv_b(:,:,Kbb) * r1_hv_b(:,:)
       !
    END SUBROUTINE istate_init
