Index: branches/UKMO/dev_r5518_obs_oper_update_icethick/NEMOGCM/NEMO/OPA_SRC/TRA/tranxt.F90
===================================================================
--- branches/UKMO/dev_r5518_obs_oper_update_icethick/NEMOGCM/NEMO/OPA_SRC/TRA/tranxt.F90	(revision 9867)
+++ branches/UKMO/dev_r5518_obs_oper_update_icethick/NEMOGCM/NEMO/OPA_SRC/TRA/tranxt.F90	(revision 9987)
@@ -28,4 +28,5 @@
    USE sbc_oce         ! surface boundary condition: ocean
    USE sbcrnf          ! river runoffs
+   USE sbcisf          ! ice shelf melting/freezing
    USE zdf_oce         ! ocean vertical mixing
    USE domvvl          ! variable volume
@@ -46,5 +47,4 @@
    USE timing          ! Timing
 #if defined key_agrif
-   USE agrif_opa_update
    USE agrif_opa_interp
 #endif
@@ -110,4 +110,8 @@
       ! 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 )
@@ -115,7 +119,4 @@
 #if defined key_bdy 
       IF( lk_bdy )   CALL bdy_tra( kt )  ! BDY open boundaries
-#endif
-#if defined key_agrif
-      CALL Agrif_tra                     ! AGRIF zoom boundaries
 #endif
  
@@ -126,12 +127,35 @@
 
       ! trends computation initialisation
-      IF( l_trdtra )   THEN                    ! store now fields before applying the Asselin filter
+      IF( l_trdtra )   THEN                    
          CALL wrk_alloc( jpi, jpj, jpk, ztrdt, ztrds )
-         ztrdt(:,:,:) = tsn(:,:,:,jp_tem) 
-         ztrds(:,:,:) = tsn(:,:,:,jp_sal)
+         ztrdt(:,:,jpk) = 0._wp
+         ztrds(:,:,jpk) = 0._wp
          IF( ln_traldf_iso ) THEN              ! diagnose the "pure" Kz diffusive trend 
             CALL trd_tra( kt, 'TRA', jp_tem, jptra_zdfp, ztrdt )
             CALL trd_tra( kt, 'TRA', jp_sal, jptra_zdfp, ztrds )
          ENDIF
+         ! total trend for the non-time-filtered variables.
+         ! G Nurser 23 Mar 2017. Recalculate trend as Delta(e3t*T)/e3tn; e3tn cancel from tsn terms
+         IF( lk_vvl ) THEN
+            DO jk = 1, jpkm1
+               zfact = 1.0 / rdttra(jk)
+               ztrdt(:,:,jk) = ( tsa(:,:,jk,jp_tem)*fse3t_a(:,:,jk) / fse3t_n(:,:,jk) - tsn(:,:,jk,jp_tem)) * zfact
+               ztrds(:,:,jk) = ( tsa(:,:,jk,jp_sal)*fse3t_a(:,:,jk) / fse3t_n(:,:,jk) - tsn(:,:,jk,jp_sal)) * zfact
+            END DO
+         ELSE
+            DO jk = 1, jpkm1
+               zfact = 1.0 / rdttra(jk)
+               ztrdt(:,:,jk) = ( tsa(:,:,jk,jp_tem) - tsn(:,:,jk,jp_tem) ) * zfact 
+               ztrds(:,:,jk) = ( tsa(:,:,jk,jp_sal) - tsn(:,:,jk,jp_sal) ) * zfact 
+            END DO
+         END IF
+         CALL trd_tra( kt, 'TRA', jp_tem, jptra_tot, ztrdt )
+         CALL trd_tra( kt, 'TRA', jp_sal, jptra_tot, ztrds )
+         IF( .NOT.lk_vvl )  THEN
+            ! Store now fields before applying the Asselin filter 
+            ! in order to calculate Asselin filter trend later.
+            ztrdt(:,:,:) = tsn(:,:,:,jp_tem) 
+            ztrds(:,:,:) = tsn(:,:,:,jp_sal)
+         END IF
       ENDIF
 
@@ -142,4 +166,11 @@
             END DO
          END DO
+         IF (l_trdtra.AND.lk_vvl) THEN      ! Zero Asselin filter contribution must be explicitly written out since for vvl
+                                            ! Asselin filter is output by tra_nxt_vvl that is not called on this time step
+            ztrdt(:,:,:) = 0._wp
+            ztrds(:,:,:) = 0._wp
+            CALL trd_tra( kt, 'TRA', jp_tem, jptra_atf, ztrdt )
+            CALL trd_tra( kt, 'TRA', jp_sal, jptra_atf, ztrds )
+         END IF
       ELSE                                            ! Leap-Frog + Asselin filter time stepping
          !
@@ -148,13 +179,8 @@
          ELSE                 ;   CALL tra_nxt_fix( kt, nit000,         'TRA', tsb, tsn, tsa, jpts )  ! fixed    volume level 
          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
-      IF( l_trdtra ) THEN      ! trend of the Asselin filter (tb filtered - tb)/dt     
+      ENDIF     
+      !
+     ! trends computation
+      IF( l_trdtra.AND..NOT.lk_vvl) THEN      ! trend of the Asselin filter (tb filtered - tb)/dt     
          DO jk = 1, jpkm1
             zfact = 1._wp / r2dtra(jk)             
@@ -164,6 +190,6 @@
          CALL trd_tra( kt, 'TRA', jp_tem, jptra_atf, ztrdt )
          CALL trd_tra( kt, 'TRA', jp_sal, jptra_atf, ztrds )
-         CALL wrk_dealloc( jpi, jpj, jpk, ztrdt, ztrds )
       END IF
+      IF( l_trdtra) CALL wrk_dealloc( jpi, jpj, jpk, ztrdt, ztrds )
       !
       !                        ! control print
@@ -279,8 +305,9 @@
 
       !!     
-      LOGICAL  ::   ll_tra_hpg, ll_traqsr, ll_rnf   ! local logical
+      LOGICAL  ::   ll_tra_hpg, 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
+      REAL(wp) ::   zfact, zfact1, ztc_a , ztc_n , ztc_b , ztc_f , ztc_d    ! local scalar
       REAL(wp) ::   zfact2, ze3t_b, ze3t_n, ze3t_a, ze3t_f, ze3t_d   !   -      -
+      REAL(wp), POINTER, DIMENSION(:,:,:,:) :: ztrd_atf
       !!----------------------------------------------------------------------
       !
@@ -295,12 +322,23 @@
          ll_traqsr  = ln_traqsr        ! active  tracers case  and  solar penetration
          ll_rnf     = ln_rnf           ! active  tracers case  and  river runoffs
+         IF (nn_isf .GE. 1) THEN 
+            ll_isf = .TRUE.            ! active  tracers case  and  ice shelf melting/freezing
+         ELSE
+            ll_isf = .FALSE.
+         END IF
       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
-      ENDIF
-      !
+         ll_isf     = .FALSE.          ! passive tracers or NO ice shelf melting/freezing
+      ENDIF
+      !
+      IF( ( l_trdtra .and. cdtype == 'TRA' ) .OR. ( l_trdtrc .and. cdtype == 'TRC' ) )   THEN
+         CALL wrk_alloc( jpi, jpj, jpk, kjpt, ztrd_atf )
+         ztrd_atf(:,:,:,:) = 0.0_wp
+      ENDIF
       DO jn = 1, kjpt      
          DO jk = 1, jpkm1
+            zfact = 1._wp / r2dtra(jk)
             zfact1 = atfp * p2dt(jk)
             zfact2 = zfact1 / rau0
@@ -321,15 +359,31 @@
                   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)
+
+                  ! 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) )  &
+                               &                 * fse3t_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) )  &
+                               &                 * fse3t_n(ji,jj,jk) * r1_hisf_tbl (ji,jj) * ralpha(ji,jj)
+                  END IF
 
                   ze3t_f = 1.e0 / ze3t_f
@@ -340,4 +394,7 @@
                      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
+                  IF( ( l_trdtra .and. cdtype == 'TRA' ) .OR. ( l_trdtrc .and. cdtype == 'TRC' ) ) THEN
+                     ztrd_atf(ji,jj,jk,jn) = (ztc_f - ztc_n) * zfact/ze3t_n
                   ENDIF
                END DO
@@ -347,4 +404,16 @@
       END DO
       !
+      IF( l_trdtra .and. cdtype == 'TRA' ) THEN 
+         CALL trd_tra( kt, cdtype, jp_tem, jptra_atf, ztrd_atf(:,:,:,jp_tem) )
+         CALL trd_tra( kt, cdtype, jp_sal, jptra_atf, ztrd_atf(:,:,:,jp_sal) )
+         CALL wrk_dealloc( jpi, jpj, jpk, kjpt, ztrd_atf )
+      ENDIF
+      IF( l_trdtrc .and. cdtype == 'TRC' ) THEN
+         DO jn = 1, kjpt
+            CALL trd_tra( kt, cdtype, jn, jptra_atf, ztrd_atf(:,:,:,jn) )
+         END DO
+         CALL wrk_dealloc( jpi, jpj, jpk, kjpt, ztrd_atf )
+      ENDIF
+
    END SUBROUTINE tra_nxt_vvl
 
