Index: branches/NERC/dev_r5549_BDY_ZEROGRAD/NEMOGCM/NEMO/OPA_SRC/TRA/tranxt.F90
===================================================================
--- branches/NERC/dev_r5549_BDY_ZEROGRAD/NEMOGCM/NEMO/OPA_SRC/TRA/tranxt.F90	(revision 6805)
+++ branches/NERC/dev_r5549_BDY_ZEROGRAD/NEMOGCM/NEMO/OPA_SRC/TRA/tranxt.F90	(revision 6808)
@@ -28,13 +28,13 @@
    USE sbc_oce         ! surface boundary condition: ocean
    USE sbcrnf          ! river runoffs
+   USE sbcisf          ! ice shelf melting
    USE zdf_oce         ! ocean vertical mixing
    USE domvvl          ! variable volume
-   USE dynspg_oce      ! surface     pressure gradient variables
-   USE dynhpg          ! hydrostatic pressure gradient 
    USE trd_oce         ! trends: ocean variables
    USE trdtra          ! trends manager: tracers 
    USE traqsr          ! penetrative solar radiation (needed for nksr)
    USE phycst          ! physical constant
-   USE ldftra_oce      ! lateral physics on tracers
+   USE ldftra          ! lateral physics on tracers
+   USE ldfslp
    USE bdy_oce         ! BDY open boundary condition variables
    USE bdytra          ! open boundary condition (bdy_tra routine)
@@ -46,5 +46,4 @@
    USE timing          ! Timing
 #if defined key_agrif
-   USE agrif_opa_update
    USE agrif_opa_interp
 #endif
@@ -57,8 +56,6 @@
    PUBLIC   tra_nxt_vvl   ! to be used in trcnxt
 
-   REAL(wp) ::   rbcp   ! Brown & Campana parameters for semi-implicit hpg
-
    !! * Substitutions
-#  include "domzgr_substitute.h90"
+#  include "vectopt_loop_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/OPA 3.3 , NEMO-Consortium (2010) 
@@ -88,11 +85,10 @@
       !!             domains (lk_agrif=T)
       !!
-      !! ** Action  : - (tb,sb) and (tn,sn) ready for the next time step
-      !!              - (ta,sa) time averaged (t,s)   (ln_dynhpg_imp = T)
+      !! ** Action  : - tsb & tsn ready for the next time step
       !!----------------------------------------------------------------------
       INTEGER, INTENT(in) ::   kt    ! ocean time-step index
       !!
-      INTEGER  ::   jk, jn    ! dummy loop indices
-      REAL(wp) ::   zfact     ! local scalars
+      INTEGER  ::   ji, jj, jk, jn   ! dummy loop indices
+      REAL(wp) ::   zfact            ! local scalars
       REAL(wp), POINTER, DIMENSION(:,:,:) ::  ztrdt, ztrds
       !!----------------------------------------------------------------------
@@ -104,10 +100,12 @@
          IF(lwp) WRITE(numout,*) 'tra_nxt : achieve the time stepping by Asselin filter and array swap'
          IF(lwp) WRITE(numout,*) '~~~~~~~'
-         !
-         rbcp = 0.25_wp * (1._wp + atfp) * (1._wp + atfp) * ( 1._wp - atfp)      ! Brown & Campana parameter for semi-implicit hpg
       ENDIF
 
       ! Update after tracer on domain lateral boundaries
       ! 
+#if defined key_agrif
+      CALL Agrif_tra                     ! AGRIF zoom boundaries
+#endif
+      !
       CALL lbc_lnk( tsa(:,:,:,jp_tem), 'T', 1._wp )      ! local domain boundaries  (T-point, unchanged sign)
       CALL lbc_lnk( tsa(:,:,:,jp_sal), 'T', 1._wp )
@@ -116,11 +114,8 @@
       IF( lk_bdy )   CALL bdy_tra( kt )  ! BDY open boundaries
 #endif
-#if defined key_agrif
-      CALL Agrif_tra                     ! AGRIF zoom boundaries
-#endif
  
       ! set time step size (Euler/Leapfrog)
-      IF( neuler == 0 .AND. kt == nit000 ) THEN   ;   r2dtra(:) =     rdttra(:)      ! at nit000             (Euler)
-      ELSEIF( kt <= nit000 + 1 )           THEN   ;   r2dtra(:) = 2._wp* rdttra(:)      ! at nit000 or nit000+1 (Leapfrog)
+      IF( neuler == 0 .AND. kt == nit000 ) THEN   ;   r2dt =     rdt      ! at nit000             (Euler)
+      ELSEIF( kt <= nit000 + 1 )           THEN   ;   r2dt = 2._wp* rdt   ! at nit000 or nit000+1 (Leapfrog)
       ENDIF
 
@@ -142,21 +137,22 @@
             END DO
          END DO
+         !
       ELSE                                            ! Leap-Frog + Asselin filter time stepping
          !
-         IF( lk_vvl )  THEN   ;   CALL tra_nxt_vvl( kt, nit000, rdttra, 'TRA', tsb, tsn, tsa,   &
-           &                                                              sbc_tsc, sbc_tsc_b, jpts )  ! variable volume level (vvl) 
-         ELSE                 ;   CALL tra_nxt_fix( kt, nit000,         'TRA', tsb, tsn, tsa, jpts )  ! fixed    volume level 
+         IF( ln_linssh ) THEN   ;   CALL tra_nxt_fix( kt, nit000,      'TRA', tsb, tsn, tsa, jpts )  ! linear free surface 
+         ELSE                   ;   CALL tra_nxt_vvl( kt, nit000, rdt, 'TRA', tsb, tsn, tsa,   &
+           &                                                                sbc_tsc, sbc_tsc_b, jpts )  ! non-linear free surface
          ENDIF
-      ENDIF 
-      !
-#if defined key_agrif
-      ! Update tracer at AGRIF zoom boundaries
-      IF( .NOT.Agrif_Root() )    CALL Agrif_Update_Tra( kt )      ! children only
-#endif      
-      !
-      ! trends computation
+         !
+         DO jn = 1, jpts
+            CALL lbc_lnk( tsb(:,:,:,jn), 'T', 1._wp ) 
+            CALL lbc_lnk( tsn(:,:,:,jn), 'T', 1._wp )
+            CALL lbc_lnk( tsa(:,:,:,jn), 'T', 1._wp )
+         END DO
+      ENDIF     
+      !
       IF( l_trdtra ) THEN      ! trend of the Asselin filter (tb filtered - tb)/dt     
          DO jk = 1, jpkm1
-            zfact = 1._wp / r2dtra(jk)             
+            zfact = 1._wp / r2dt             
             ztrdt(:,:,jk) = ( tsb(:,:,jk,jp_tem) - ztrdt(:,:,jk) ) * zfact
             ztrds(:,:,jk) = ( tsb(:,:,jk,jp_sal) - ztrds(:,:,jk) ) * zfact
@@ -184,32 +180,20 @@
       !! 
       !! ** Method  : - Apply a Asselin time filter on now fields.
-      !!              - save in (ta,sa) an average over the three time levels 
-      !!             which will be used to compute rdn and thus the semi-implicit
-      !!             hydrostatic pressure gradient (ln_dynhpg_imp = T)
       !!              - swap tracer fields to prepare the next time_step.
-      !!                This can be summurized for tempearture as:
-      !!             ztm = tn + rbcp * [ta -2 tn + tb ]       ln_dynhpg_imp = T
-      !!             ztm = 0                                   otherwise
-      !!                   with rbcp=1/4 * (1-atfp^4) / (1-atfp)
-      !!             tb  = tn + atfp*[ tb - 2 tn + ta ]
-      !!             tn  = ta  
-      !!             ta  = ztm       (NB: reset to 0 after eos_bn2 call)
-      !!
-      !! ** Action  : - (tb,sb) and (tn,sn) ready for the next time step
-      !!              - (ta,sa) time averaged (t,s)   (ln_dynhpg_imp = T)
-      !!----------------------------------------------------------------------
-      INTEGER         , INTENT(in   )                               ::   kt       ! ocean time-step index
-      INTEGER         , INTENT(in   )                               ::   kit000   ! first time step index
-      CHARACTER(len=3), INTENT(in   )                               ::   cdtype   ! =TRA or TRC (tracer indicator)
-      INTEGER         , INTENT(in   )                               ::   kjpt     ! number of tracers
-      REAL(wp)        , INTENT(inout), DIMENSION(jpi,jpj,jpk,kjpt)  ::   ptb      ! before tracer fields
-      REAL(wp)        , INTENT(inout), DIMENSION(jpi,jpj,jpk,kjpt)  ::   ptn      ! now tracer fields
-      REAL(wp)        , INTENT(inout), DIMENSION(jpi,jpj,jpk,kjpt)  ::   pta      ! tracer trend
+      !!
+      !! ** Action  : - tsb & tsn ready for the next time step
+      !!----------------------------------------------------------------------
+      INTEGER                              , INTENT(in   ) ::  kt        ! ocean time-step index
+      INTEGER                              , INTENT(in   ) ::  kit000    ! first time step index
+      CHARACTER(len=3)                     , INTENT(in   ) ::  cdtype    ! =TRA or TRC (tracer indicator)
+      INTEGER                              , INTENT(in   ) ::  kjpt      ! number of tracers
+      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::  ptb       ! before tracer fields
+      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::  ptn       ! now tracer fields
+      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::  pta       ! tracer trend
       !
       INTEGER  ::   ji, jj, jk, jn   ! dummy loop indices
-      LOGICAL  ::   ll_tra_hpg       ! local logical
       REAL(wp) ::   ztn, ztd         ! local scalars
       !!----------------------------------------------------------------------
-
+      !
       IF( kt == kit000 )  THEN
          IF(lwp) WRITE(numout,*)
@@ -218,20 +202,14 @@
       ENDIF
       !
-      IF( cdtype == 'TRA' )  THEN   ;   ll_tra_hpg = ln_dynhpg_imp    ! active  tracers case  and  semi-implicit hpg    
-      ELSE                          ;   ll_tra_hpg = .FALSE.          ! passive tracers case or NO semi-implicit hpg
-      ENDIF
-      !
       DO jn = 1, kjpt
          !
          DO jk = 1, jpkm1
-            DO jj = 1, jpj
-               DO ji = 1, jpi
+            DO jj = 2, jpjm1
+               DO ji = fs_2, fs_jpim1
                   ztn = ptn(ji,jj,jk,jn)                                    
-                  ztd = pta(ji,jj,jk,jn) - 2. * ztn + ptb(ji,jj,jk,jn)      !  time laplacian on tracers
-                  !
-                  ptb(ji,jj,jk,jn) = ztn + atfp * ztd                       ! ptb <-- filtered ptn 
-                  ptn(ji,jj,jk,jn) = pta(ji,jj,jk,jn)                       ! ptn <-- pta
-                  !
-                  IF( ll_tra_hpg )   pta(ji,jj,jk,jn) = ztn + rbcp * ztd    ! pta <-- Brown & Campana average
+                  ztd = pta(ji,jj,jk,jn) - 2._wp * ztn + ptb(ji,jj,jk,jn)  ! time laplacian on tracers
+                  !
+                  ptb(ji,jj,jk,jn) = ztn + atfp * ztd                      ! ptb <-- filtered ptn 
+                  ptn(ji,jj,jk,jn) = pta(ji,jj,jk,jn)                      ! ptn <-- pta
                END DO
            END DO
@@ -251,33 +229,23 @@
       !! 
       !! ** Method  : - Apply a thickness weighted Asselin time filter on now fields.
-      !!              - save in (ta,sa) a thickness weighted average over the three 
-      !!             time levels which will be used to compute rdn and thus the semi-
-      !!             implicit hydrostatic pressure gradient (ln_dynhpg_imp = T)
       !!              - swap tracer fields to prepare the next time_step.
-      !!                This can be summurized for tempearture as:
-      !!             ztm = ( e3t_n*tn + rbcp*[ e3t_b*tb - 2 e3t_n*tn + e3t_a*ta ] )   ln_dynhpg_imp = T
-      !!                  /( e3t_n    + rbcp*[ e3t_b    - 2 e3t_n    + e3t_a    ] )   
-      !!             ztm = 0                                                       otherwise
       !!             tb  = ( e3t_n*tn + atfp*[ e3t_b*tb - 2 e3t_n*tn + e3t_a*ta ] )
       !!                  /( e3t_n    + atfp*[ e3t_b    - 2 e3t_n    + e3t_a    ] )
       !!             tn  = ta 
-      !!             ta  = zt        (NB: reset to 0 after eos_bn2 call)
-      !!
-      !! ** Action  : - (tb,sb) and (tn,sn) ready for the next time step
-      !!              - (ta,sa) time averaged (t,s)   (ln_dynhpg_imp = T)
-      !!----------------------------------------------------------------------
-      INTEGER         , INTENT(in   )                               ::  kt       ! ocean time-step index
-      INTEGER         , INTENT(in   )                               ::  kit000   ! first time step index
-      REAL(wp)        , INTENT(in   ), DIMENSION(jpk)               ::  p2dt     ! time-step
-      CHARACTER(len=3), INTENT(in   )                               ::  cdtype   ! =TRA or TRC (tracer indicator)
-      INTEGER         , INTENT(in   )                               ::  kjpt     ! number of tracers
-      REAL(wp)        , INTENT(inout), DIMENSION(jpi,jpj,jpk,kjpt)  ::  ptb      ! before tracer fields
-      REAL(wp)        , INTENT(inout), DIMENSION(jpi,jpj,jpk,kjpt)  ::  ptn      ! now tracer fields
-      REAL(wp)        , INTENT(inout), DIMENSION(jpi,jpj,jpk,kjpt)  ::  pta      ! tracer trend
-      REAL(wp)        , INTENT(in   ), DIMENSION(jpi,jpj,kjpt)      ::  psbc_tc   ! surface tracer content
-      REAL(wp)        , INTENT(in   ), DIMENSION(jpi,jpj,kjpt)      ::  psbc_tc_b ! before surface tracer content
-
-      !!     
-      LOGICAL  ::   ll_tra_hpg, ll_traqsr, ll_rnf   ! local logical
+      !!
+      !! ** Action  : - tsb & tsn ready for the next time step
+      !!----------------------------------------------------------------------
+      INTEGER                              , INTENT(in   ) ::  kt        ! ocean time-step index
+      INTEGER                              , INTENT(in   ) ::  kit000    ! first time step index
+      REAL(wp)                             , INTENT(in   ) ::  p2dt      ! time-step
+      CHARACTER(len=3)                     , INTENT(in   ) ::  cdtype    ! =TRA or TRC (tracer indicator)
+      INTEGER                              , INTENT(in   ) ::  kjpt      ! number of tracers
+      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::  ptb       ! before tracer fields
+      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::  ptn       ! now tracer fields
+      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::  pta       ! tracer trend
+      REAL(wp), DIMENSION(jpi,jpj    ,kjpt), INTENT(in   ) ::  psbc_tc   ! surface tracer content
+      REAL(wp), DIMENSION(jpi,jpj    ,kjpt), INTENT(in   ) ::  psbc_tc_b ! before surface tracer content
+      !
+      LOGICAL  ::   ll_traqsr, ll_rnf, ll_isf   ! local logical
       INTEGER  ::   ji, jj, jk, jn              ! dummy loop indices
       REAL(wp) ::   zfact1, ztc_a , ztc_n , ztc_b , ztc_f , ztc_d    ! local scalar
@@ -292,22 +260,22 @@
       !
       IF( cdtype == 'TRA' )  THEN   
-         ll_tra_hpg = ln_dynhpg_imp    ! active  tracers case  and  semi-implicit hpg
          ll_traqsr  = ln_traqsr        ! active  tracers case  and  solar penetration
          ll_rnf     = ln_rnf           ! active  tracers case  and  river runoffs
-      ELSE                          
-         ll_tra_hpg = .FALSE.          ! passive tracers case or NO semi-implicit hpg
-         ll_traqsr  = .FALSE.          ! active  tracers case and NO solar penetration
-         ll_rnf     = .FALSE.          ! passive tracers or NO river runoffs
+         ll_isf     = ln_isf           ! active  tracers case  and  ice shelf melting
+      ELSE                          ! passive tracers case
+         ll_traqsr  = .FALSE.          ! NO solar penetration
+         ll_rnf     = .FALSE.          ! NO river runoffs ????          !!gm BUG ?  
+         ll_isf     = .FALSE.          ! NO ice shelf melting/freezing  !!gm BUG ?? 
       ENDIF
       !
       DO jn = 1, kjpt      
          DO jk = 1, jpkm1
-            zfact1 = atfp * p2dt(jk)
-            zfact2 = zfact1 / rau0
-            DO jj = 1, jpj
-               DO ji = 1, jpi
-                  ze3t_b = fse3t_b(ji,jj,jk)
-                  ze3t_n = fse3t_n(ji,jj,jk)
-                  ze3t_a = fse3t_a(ji,jj,jk)
+            zfact1 = atfp * p2dt
+            zfact2 = zfact1 * r1_rau0
+            DO jj = 2, jpjm1
+               DO ji = fs_2, fs_jpim1
+                  ze3t_b = e3t_b(ji,jj,jk)
+                  ze3t_n = e3t_n(ji,jj,jk)
+                  ze3t_a = e3t_a(ji,jj,jk)
                   !                                         ! tracer content at Before, now and after
                   ztc_b  = ptb(ji,jj,jk,jn) * ze3t_b
@@ -321,24 +289,36 @@
                   ztc_f  = ztc_n  + atfp * ztc_d
                   !
-                  IF( jk == 1 ) THEN           ! first level 
-                     ze3t_f = ze3t_f - zfact2 * ( emp_b(ji,jj) - emp(ji,jj) + rnf(ji,jj) - rnf_b(ji,jj) )
+                  IF( jk == mikt(ji,jj) ) THEN           ! first level 
+                     ze3t_f = ze3t_f - zfact2 * ( (emp_b(ji,jj)    - emp(ji,jj)   )  &
+                            &                   - (rnf_b(ji,jj)    - rnf(ji,jj)   )  &
+                            &                   + (fwfisf_b(ji,jj) - fwfisf(ji,jj))  )
                      ztc_f  = ztc_f  - zfact1 * ( psbc_tc(ji,jj,jn) - psbc_tc_b(ji,jj,jn) )
                   ENDIF
-
-                  IF( ll_traqsr .AND. jn == jp_tem .AND. jk <= nksr )   &     ! solar penetration (temperature only)
+                  !
+                  ! solar penetration (temperature only)
+                  IF( ll_traqsr .AND. jn == jp_tem .AND. jk <= nksr )                            & 
                      &     ztc_f  = ztc_f  - zfact1 * ( qsr_hc(ji,jj,jk) - qsr_hc_b(ji,jj,jk) ) 
-
-                  IF( ll_rnf .AND. jk <= nk_rnf(ji,jj) )   & 		      ! river runoffs
+                     !
+                  ! river runoff
+                  IF( ll_rnf .AND. jk <= nk_rnf(ji,jj) )                                          &
                      &     ztc_f  = ztc_f  - zfact1 * ( rnf_tsc(ji,jj,jn) - rnf_tsc_b(ji,jj,jn) ) & 
-                     &                              * fse3t_n(ji,jj,jk) / h_rnf(ji,jj)
-
+                     &                              * e3t_n(ji,jj,jk) / h_rnf(ji,jj)
+                     !
+                  ! ice shelf
+                  IF( ll_isf ) THEN
+                     ! level fully include in the Losch_2008 ice shelf boundary layer
+                     IF ( jk >= misfkt(ji,jj) .AND. jk < misfkb(ji,jj) )                          &
+                        ztc_f  = ztc_f  - zfact1 * ( risf_tsc(ji,jj,jn) - risf_tsc_b(ji,jj,jn) )  &
+                               &                 * e3t_n(ji,jj,jk) * r1_hisf_tbl (ji,jj)
+                     ! level partially include in Losch_2008 ice shelf boundary layer 
+                     IF ( jk == misfkb(ji,jj) )                                                   &
+                        ztc_f  = ztc_f  - zfact1 * ( risf_tsc(ji,jj,jn) - risf_tsc_b(ji,jj,jn) )  &
+                               &                 * e3t_n(ji,jj,jk) * r1_hisf_tbl (ji,jj) * ralpha(ji,jj)
+                  END IF
+                  !
                   ze3t_f = 1.e0 / ze3t_f
                   ptb(ji,jj,jk,jn) = ztc_f * ze3t_f       ! ptb <-- ptn filtered
                   ptn(ji,jj,jk,jn) = pta(ji,jj,jk,jn)     ! ptn <-- pta
                   !
-                  IF( ll_tra_hpg ) THEN        ! semi-implicit hpg (T & S only)
-                     ze3t_d           = 1.e0   / ( ze3t_n + rbcp * ze3t_d )
-                     pta(ji,jj,jk,jn) = ze3t_d * ( ztc_n  + rbcp * ztc_d  )   ! ta <-- Brown & Campana average
-                  ENDIF
                END DO
             END DO
