Index: /trunk/NEMOGCM/NEMO/OPA_SRC/DIA/diahsb.F90
===================================================================
--- /trunk/NEMOGCM/NEMO/OPA_SRC/DIA/diahsb.F90	(revision 4557)
+++ /trunk/NEMOGCM/NEMO/OPA_SRC/DIA/diahsb.F90	(revision 4558)
@@ -36,8 +36,10 @@
    LOGICAL, PUBLIC ::   ln_diahsb   !: check the heat and salt budgets
 
-   REAL(wp), SAVE                                ::   frc_t      , frc_s     , frc_v   ! global forcing trends
-   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   ssh_ini              !
-   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   hc_loc_ini, sc_loc_ini, e3t_ini  !
-   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   hcssh_loc_ini, scssh_loc_ini     !
+   REAL(dp)                                ::   surf_tot                !
+   REAL(dp)                                ::   frc_t      , frc_s     , frc_v   ! global forcing trends
+   REAL(dp)                                ::   frc_wn_t      , frc_wn_s ! global forcing trends
+   REAL(dp), DIMENSION(:,:)  , ALLOCATABLE ::   surf      , ssh_ini              !
+   REAL(dp), DIMENSION(:,:,:), ALLOCATABLE ::   hc_loc_ini, sc_loc_ini, e3t_ini  !
+   REAL(dp), DIMENSION(:,:)  , ALLOCATABLE ::   ssh_hc_loc_ini, ssh_sc_loc_ini
 
    !! * Substitutions
@@ -67,92 +69,83 @@
       !!
       INTEGER    ::   jk                          ! dummy loop indice
-      REAL(wp)   ::   zdiff_hc    , zdiff_sc      ! heat and salt content variations
-      REAL(wp)   ::   zdiff_v1    , zdiff_v2      ! volume variation
-      REAL(wp)   ::   z_hc        , z_sc          ! heat and salt content
-      REAL(wp)   ::   z_v1        , z_v2          ! volume
-      REAL(wp)   ::   zdeltat                     !    -     -
-      REAL(wp)   ::   z_frc_trd_t , z_frc_trd_s   !    -     -
-      REAL(wp)   ::   z_frc_trd_v                 !    -     -
-      REAL(wp), POINTER, DIMENSION(:,:)   ::   zsurf              !
+      REAL(dp)   ::   zdiff_hc    , zdiff_sc      ! heat and salt content variations
+      REAL(dp)   ::   zdiff_hc1   , zdiff_sc1     ! -   -   -   -   -   -   -   - 
+      REAL(dp)   ::   zdiff_v1    , zdiff_v2      ! volume variation
+      REAL(dp)   ::   zerr_hc1    , zerr_sc1       ! heat and salt content misfit
+      REAL(dp)   ::   zvol_tot                    ! volume
+      REAL(dp)   ::   z_frc_trd_t , z_frc_trd_s   !    -     -
+      REAL(dp)   ::   z_frc_trd_v                 !    -     -
+      REAL(dp)   ::   z_wn_trd_t , z_wn_trd_s   !    -     -
+      REAL(dp)   ::   z_ssh_hc , z_ssh_sc   !    -     -
       !!---------------------------------------------------------------------------
       IF( nn_timing == 1 )   CALL timing_start('dia_hsb')      
 
-      CALL wrk_alloc( jpi, jpj, zsurf )
-  
-      zsurf(:,:) = e1t(:,:) * e2t(:,:) * tmask(:,:,1) * tmask_i(:,:)      ! masked surface grid cell area
-      
       ! ------------------------- !
       ! 1 - Trends due to forcing !
       ! ------------------------- !
-      z_frc_trd_v = r1_rau0 * glob_sum( - ( emp(:,:) - rnf(:,:) ) * zsurf(:,:) ) ! volume fluxes
-      z_frc_trd_t =           glob_sum( sbc_tsc(:,:,jp_tem) * zsurf(:,:) )       ! heat fluxes
-      z_frc_trd_s =           glob_sum( sbc_tsc(:,:,jp_sal) * zsurf(:,:) )       ! salt fluxes
-      !
-      IF( ln_rnf    )   z_frc_trd_t = z_frc_trd_t + glob_sum( rnf_tsc(:,:,jp_tem) * zsurf(:,:) )
-      IF( ln_rnf_sal)   z_frc_trd_s = z_frc_trd_s + glob_sum( rnf_tsc(:,:,jp_sal) * zsurf(:,:) )
+      z_frc_trd_v = r1_rau0 * glob_sum( - ( emp(:,:) - rnf(:,:) ) * surf(:,:) ) ! volume fluxes
+      z_frc_trd_t =           glob_sum( sbc_tsc(:,:,jp_tem) * surf(:,:) )       ! heat fluxes
+      z_frc_trd_s =           glob_sum( sbc_tsc(:,:,jp_sal) * surf(:,:) )       ! salt fluxes
+      ! Add runoff heat & salt input
+      IF( ln_rnf    )   z_frc_trd_t = z_frc_trd_t + glob_sum( rnf_tsc(:,:,jp_tem) * surf(:,:) )
+      IF( ln_rnf_sal)   z_frc_trd_s = z_frc_trd_s + glob_sum( rnf_tsc(:,:,jp_sal) * surf(:,:) )
 
       ! Add penetrative solar radiation
-      IF( ln_traqsr )   z_frc_trd_t = z_frc_trd_t + r1_rau0_rcp * glob_sum( qsr     (:,:) * zsurf(:,:) )
+      IF( ln_traqsr )   z_frc_trd_t = z_frc_trd_t + r1_rau0_rcp * glob_sum( qsr     (:,:) * surf(:,:) )
       ! Add geothermal heat flux
-      IF( ln_trabbc )   z_frc_trd_t = z_frc_trd_t +               glob_sum( qgh_trd0(:,:) * zsurf(:,:) )
-      !
+      IF( ln_trabbc )   z_frc_trd_t = z_frc_trd_t +               glob_sum( qgh_trd0(:,:) * surf(:,:) )
+      !
+      IF( .NOT. lk_vvl ) THEN
+         z_wn_trd_t = - glob_sum( surf(:,:) * wn(:,:,1) * tsb(:,:,1,jp_tem) )
+         z_wn_trd_s = - glob_sum( surf(:,:) * wn(:,:,1) * tsb(:,:,1,jp_sal) )
+      ENDIF
+
       frc_v = frc_v + z_frc_trd_v * rdt
       frc_t = frc_t + z_frc_trd_t * rdt
       frc_s = frc_s + z_frc_trd_s * rdt
+      !                                          ! Advection flux through fixed surface (z=0)
+      IF( .NOT. lk_vvl ) THEN
+         frc_wn_t = frc_wn_t + z_wn_trd_t * rdt
+         frc_wn_s = frc_wn_s + z_wn_trd_s * rdt
+      ENDIF
 
       ! ------------------------ !
-      ! 2a -  Content variations !
+      ! 2 -  Content variations !
       ! ------------------------ !
-      zdiff_v2 = 0._wp
-      zdiff_hc = 0._wp
-      zdiff_sc = 0._wp
+      zdiff_v2 = 0.d0
+      zdiff_hc = 0.d0
+      zdiff_sc = 0.d0
+
       ! volume variation (calculated with ssh)
-      zdiff_v1 = glob_sum( zsurf(:,:) * ( sshn(:,:) - ssh_ini(:,:) ) )
+      zdiff_v1 = glob_sum( surf(:,:) * ( sshn(:,:) - ssh_ini(:,:) ) )
+
+      ! heat & salt content variation (associated with ssh)
+      IF( .NOT. lk_vvl ) THEN
+         z_ssh_hc = glob_sum( surf(:,:) * ( tsn(:,:,1,jp_tem) * sshn(:,:) - ssh_hc_loc_ini(:,:) ) )
+         z_ssh_sc = glob_sum( surf(:,:) * ( tsn(:,:,1,jp_sal) * sshn(:,:) - ssh_sc_loc_ini(:,:) ) )
+      ENDIF
+
       DO jk = 1, jpkm1
          ! volume variation (calculated with scale factors)
-         zdiff_v2 = zdiff_v2 + glob_sum( zsurf(:,:) * tmask(:,:,jk) * ( fse3t_n(:,:,jk) - e3t_ini(:,:,jk) ) )
+         zdiff_v2 = zdiff_v2 + glob_sum( surf(:,:) * tmask(:,:,jk) &
+            &                           * ( fse3t_n(:,:,jk) - e3t_ini(:,:,jk) ) )
          ! heat content variation
-         zdiff_hc = zdiff_hc + glob_sum( zsurf(:,:) * tmask(:,:,jk) * ( fse3t_n(:,:,jk) * tsn(:,:,jk,jp_tem)   &
-            &                           - hc_loc_ini(:,:,jk) ) )
+         zdiff_hc = zdiff_hc + glob_sum(  surf(:,:) * tmask(:,:,jk) & 
+            &                           * ( fse3t_n(:,:,jk) * tsn(:,:,jk,jp_tem) - hc_loc_ini(:,:,jk) ) )
          ! salt content variation
-         zdiff_sc = zdiff_sc + glob_sum( zsurf(:,:) * tmask(:,:,jk) * ( fse3t_n(:,:,jk) * tsn(:,:,jk,jp_sal)   &
-            &                           - sc_loc_ini(:,:,jk) ) )
+         zdiff_sc = zdiff_sc + glob_sum(  surf(:,:) * tmask(:,:,jk)   &
+            &                           * ( fse3t_n(:,:,jk) * tsn(:,:,jk,jp_sal) - sc_loc_ini(:,:,jk) ) )
       ENDDO
 
       ! Substract forcing from heat content, salt content and volume variations
-      !frc_v = zdiff_v2 - frc_v
-      !frc_t = zdiff_hc - frc_t
-      !frc_s = zdiff_sc - frc_s
-      
-      ! add ssh if not vvl
+      zdiff_v1 = zdiff_v1 - frc_v
+      IF( lk_vvl )   zdiff_v2 = zdiff_v2 - frc_v
+      zdiff_hc = zdiff_hc - frc_t
+      zdiff_sc = zdiff_sc - frc_s
       IF( .NOT. lk_vvl ) THEN
-        zdiff_v2 = zdiff_v2 + zdiff_v1
-        zdiff_hc = zdiff_hc + glob_sum( zsurf(:,:) * ( sshn(:,:) * tsn(:,:,1,jp_tem)   &
-               &                           - hcssh_loc_ini(:,:) ) )
-        zdiff_sc = zdiff_sc + glob_sum( zsurf(:,:) * ( sshn(:,:) * tsn(:,:,1,jp_sal)   &
-               &                           - scssh_loc_ini(:,:) ) )
-      ENDIF
-      !
-      ! ----------------------- !
-      ! 2b -  Content           !
-      ! ----------------------- !
-      z_v2 = 0._wp
-      z_hc = 0._wp
-      z_sc = 0._wp
-      ! volume (calculated with ssh)
-      z_v1 = glob_sum( zsurf(:,:) * sshn(:,:) )
-      DO jk = 1, jpkm1
-         ! volume (calculated with scale factors)
-         z_v2 = z_v2 + glob_sum( zsurf(:,:) * tmask(:,:,jk) * fse3t_n(:,:,jk) )
-         ! heat content
-         z_hc = z_hc + glob_sum( zsurf(:,:) * tmask(:,:,jk) * fse3t_n(:,:,jk) * tsn(:,:,jk,jp_tem) )
-         ! salt content
-         z_sc = z_sc + glob_sum( zsurf(:,:) * tmask(:,:,jk) * fse3t_n(:,:,jk) * tsn(:,:,jk,jp_sal) )
-      ENDDO
-      ! add ssh if not vvl
-      IF( .NOT. lk_vvl ) THEN
-        z_v2 = z_v2 + z_v1
-        z_hc = z_hc + glob_sum( zsurf(:,:) * sshn(:,:) * tsn(:,:,1,jp_tem) )
-        z_sc = z_sc + glob_sum( zsurf(:,:) * sshn(:,:) * tsn(:,:,1,jp_sal) )
+         zdiff_hc1 = zdiff_hc + z_ssh_hc 
+         zdiff_sc1 = zdiff_sc + z_ssh_sc
+         zerr_hc1  = z_ssh_hc - frc_wn_t
+         zerr_sc1  = z_ssh_sc - frc_wn_s
       ENDIF
 
@@ -160,19 +153,34 @@
       ! 3 - Diagnostics writing !
       ! ----------------------- !
-      zdeltat  = 1.e0 / ( ( kt - nit000 + 1 ) * rdt )
-!
-      CALL iom_put( 'bgtemper' , z_hc / z_v2 )                      ! Temperature (C) 
-      CALL iom_put( 'bgsaline' , z_sc / z_v2 )                      ! Salinity (psu)
-      CALL iom_put( 'bgheatco' , zdiff_hc * rau0 * rcp * 1.e-9_wp ) ! Heat content variation (10^9 J)
-      CALL iom_put( 'bgsaltco' , zdiff_sc * 1.e-9 )                 ! Salt content variation (psu*km3) 
-      CALL iom_put( 'bgvolssh' , zdiff_v1 * 1.e-9 )                    ! volume ssh (km3)  
-      CALL iom_put( 'bgsshtot' , zdiff_v1 / glob_sum(zsurf) )          ! ssh (m)  
-      CALL iom_put( 'bgvoltot' , zdiff_v2 * 1.e-9 )                 ! volume total (km3) 
-      CALL iom_put( 'bgfrcvol' , frc_v * 1.e-9 )                     ! vol - surface forcing (volume) 
-      CALL iom_put( 'bgfrctem' , frc_t * rau0 * rcp * 1.e-9_wp ) ! hc  - surface forcing (heat content) 
-      CALL iom_put( 'bgfrcsal' , frc_s * 1.e-9 )                     ! sc  - surface forcing (salt content) 
-      !
-      CALL wrk_dealloc( jpi, jpj, zsurf )
-      !
+      zvol_tot   = 0.d0                                                   ! total ocean volume
+      DO jk = 1, jpkm1
+         zvol_tot  = zvol_tot + glob_sum( surf(:,:) * tmask(:,:,jk) * fse3t_n(:,:,jk) )
+      END DO
+
+      IF( lk_vvl ) THEN
+        CALL iom_put( 'bgtemper' , zdiff_hc / zvol_tot )              ! Temperature variation (C) 
+        CALL iom_put( 'bgsaline' , zdiff_sc / zvol_tot )              ! Salinity    variation (psu)
+        CALL iom_put( 'bgheatco' , zdiff_hc * 1.e-20 * rau0 * rcp )   ! Heat content variation (1.e20 J) 
+        CALL iom_put( 'bgsaltco' , zdiff_sc * 1.e-9    )              ! Salt content variation (psu*km3)
+        CALL iom_put( 'bgvolssh' , zdiff_v1 * 1.e-9    )              ! volume ssh variation (km3)  
+        CALL iom_put( 'bgvole3t' , zdiff_v2 * 1.e-9    )              ! volume e3t variation (km3)  
+        CALL iom_put( 'bgfrcvol' , frc_v    * 1.e-9    )              ! vol - surface forcing (km3) 
+        CALL iom_put( 'bgfrctem' , frc_t / zvol_tot    )              ! hc  - surface forcing (C) 
+        CALL iom_put( 'bgfrcsal' , frc_s / zvol_tot    )              ! sc  - surface forcing (psu) 
+      ELSE
+        CALL iom_put( 'bgtemper' , zdiff_hc1 / zvol_tot)              ! Heat content variation (C) 
+        CALL iom_put( 'bgsaline' , zdiff_sc1 / zvol_tot)              ! Salt content variation (psu)
+        CALL iom_put( 'bgheatco' , zdiff_hc1 * 1.e-20 * rau0 * rcp )  ! Heat content variation (1.e20 J) 
+        CALL iom_put( 'bgsaltco' , zdiff_sc1 * 1.e-9    )             ! Salt content variation (psu*km3)
+        CALL iom_put( 'bgvolssh' , zdiff_v1 * 1.e-9    )              ! volume ssh variation (km3)  
+        CALL iom_put( 'bgfrcvol' , frc_v    * 1.e-9    )              ! vol - surface forcing (km3) 
+        CALL iom_put( 'bgfrctem' , frc_t / zvol_tot    )              ! hc  - surface forcing (C) 
+        CALL iom_put( 'bgfrcsal' , frc_s / zvol_tot    )              ! sc  - surface forcing (psu) 
+        CALL iom_put( 'bgmistem' , zerr_hc1 / zvol_tot )              ! hc  - error due to free surface (C)
+        CALL iom_put( 'bgmissal' , zerr_sc1 / zvol_tot )              ! sc  - error due to free surface (psu)
+      ENDIF
+      !
+      IF( lrst_oce )   CALL dia_hsb_rst( kt, 'WRITE' )
+
       IF( nn_timing == 1 )   CALL timing_stop('dia_hsb')
 !
@@ -222,45 +230,41 @@
          WRITE(numout,*) '   Namelist namhsb : set hsb parameters'
          WRITE(numout,*) '      Switch for hsb diagnostic (T) or not (F)  ln_diahsb  = ', ln_diahsb
+         WRITE(numout,*)
       ENDIF
 
       IF( .NOT. ln_diahsb )   RETURN
-
-         ! ------------------- !
-         ! 1 - Allocate memory !
-         ! ------------------- !
-         ALLOCATE( hc_loc_ini(jpi,jpj,jpk), STAT=ierror )
-         IF( ierror > 0 ) THEN
-            CALL ctl_stop( 'dia_hsb: unable to allocate hc_loc_ini' )   ;   RETURN
-         ENDIF
-         ALLOCATE( sc_loc_ini(jpi,jpj,jpk), STAT=ierror )
-         IF( ierror > 0 ) THEN
-            CALL ctl_stop( 'dia_hsb: unable to allocate sc_loc_ini' )   ;   RETURN
-         ENDIF
-         ALLOCATE( hcssh_loc_ini(jpi,jpj), STAT=ierror )
-         IF( ierror > 0 ) THEN
-            CALL ctl_stop( 'dia_hsb: unable to allocate hcssh_loc_ini' )   ;   RETURN
-         ENDIF
-         ALLOCATE( scssh_loc_ini(jpi,jpj), STAT=ierror )
-         IF( ierror > 0 ) THEN
-            CALL ctl_stop( 'dia_hsb: unable to allocate scssh_loc_ini' )   ;   RETURN
-         ENDIF
-         ALLOCATE( e3t_ini(jpi,jpj,jpk)   , STAT=ierror )
-         IF( ierror > 0 ) THEN
-            CALL ctl_stop( 'dia_hsb: unable to allocate e3t_ini' )      ;   RETURN
-         ENDIF
-         ALLOCATE( ssh_ini(jpi,jpj)       , STAT=ierror )
-         IF( ierror > 0 ) THEN
-            CALL ctl_stop( 'dia_hsb: unable to allocate ssh_ini' )      ;   RETURN
-         ENDIF
-         
-         ! ----------------------------------------------- !
-         ! 2 - Time independant variables and file opening !
-         ! ----------------------------------------------- !
-         IF(lwp) WRITE(numout,*) "dia_hsb: heat salt volume budgets activated"
-         IF( lk_bdy ) THEN
-            CALL ctl_warn( 'dia_hsb does not take open boundary fluxes into account' )         
-         ENDIF
-         !
-         CALL dia_hsb_rst( nit000, 'READ' )  !* read or initialize all required files
+!      IF( .NOT. lk_mpp_rep ) &
+!        CALL ctl_stop (' Your global mpp_sum if performed in single precision - 64 bits -', &
+!             &         ' whereas the global sum to be precise must be done in double precision ',&
+!             &         ' please add key_mpp_rep')
+
+      ! ------------------- !
+      ! 1 - Allocate memory !
+      ! ------------------- !
+      ALLOCATE( hc_loc_ini(jpi,jpj,jpk), sc_loc_ini(jpi,jpj,jpk), &
+         &      e3t_ini(jpi,jpj,jpk), surf(jpi,jpj),  ssh_ini(jpi,jpj), STAT=ierror )
+      IF( ierror > 0 ) THEN
+         CALL ctl_stop( 'dia_hsb: unable to allocate hc_loc_ini' )   ;   RETURN
+      ENDIF
+
+      IF(.NOT. lk_vvl ) ALLOCATE( ssh_hc_loc_ini(jpi,jpj), ssh_sc_loc_ini(jpi,jpj),STAT=ierror )
+      IF( ierror > 0 ) THEN
+         CALL ctl_stop( 'dia_hsb: unable to allocate hc_loc_ini' )   ;   RETURN
+      ENDIF
+
+      ! ----------------------------------------------- !
+      ! 2 - Time independant variables and file opening !
+      ! ----------------------------------------------- !
+      IF(lwp) WRITE(numout,*) "dia_hsb: heat salt volume budgets activated"
+      IF(lwp) WRITE(numout,*) '~~~~~~~'
+      surf(:,:) = e1t(:,:) * e2t(:,:) * tmask(:,:,1) * tmask_i(:,:)      ! masked surface grid cell area
+      surf_tot  = glob_sum( surf(:,:) )                                       ! total ocean surface area
+
+      IF( lk_bdy ) CALL ctl_warn( 'dia_hsb does not take open boundary fluxes into account' )         
+      !
+      ! ---------------------------------- !
+      ! 4 - initial conservation variables !
+      ! ---------------------------------- !
+      CALL dia_hsb_rst( nit000, 'READ' )  !* read or initialize all required files
       !
    END SUBROUTINE dia_hsb_init
@@ -285,15 +289,26 @@
            !id1 = iom_varid( numror, 'frc_vol'  , ldstop = .FALSE. )
            !
+           IF(lwp) WRITE(numout,*) '~~~~~~~'
+           IF(lwp) WRITE(numout,*) ' dia_hsb_rst at it= ', kt,' date= ', ndastp
+           IF(lwp) WRITE(numout,*) '~~~~~~~'
            CALL iom_get( numror, 'frc_v', frc_v )
            CALL iom_get( numror, 'frc_t', frc_t )
            CALL iom_get( numror, 'frc_s', frc_s )
-
+           IF( .NOT. lk_vvl ) THEN
+              CALL iom_get( numror, 'frc_wn_t', frc_wn_t )
+              CALL iom_get( numror, 'frc_wn_s', frc_wn_s )
+           ENDIF
            CALL iom_get( numror, jpdom_autoglo, 'ssh_ini', ssh_ini )
            CALL iom_get( numror, jpdom_autoglo, 'e3t_ini', e3t_ini )
            CALL iom_get( numror, jpdom_autoglo, 'hc_loc_ini', hc_loc_ini )
            CALL iom_get( numror, jpdom_autoglo, 'sc_loc_ini', sc_loc_ini )
-           CALL iom_get( numror, jpdom_autoglo, 'hcssh_loc_ini', hcssh_loc_ini )
-           CALL iom_get( numror, jpdom_autoglo, 'scssh_loc_ini', scssh_loc_ini )
+           IF( .NOT. lk_vvl ) THEN
+              CALL iom_get( numror, jpdom_autoglo, 'ssh_hc_loc_ini', ssh_hc_loc_ini )
+              CALL iom_get( numror, jpdom_autoglo, 'ssh_sc_loc_ini', ssh_sc_loc_ini )
+           ENDIF
        ELSE
+          IF(lwp) WRITE(numout,*) '~~~~~~~'
+          IF(lwp) WRITE(numout,*) ' dia_hsb at initial state '
+          IF(lwp) WRITE(numout,*) '~~~~~~~'
           ssh_ini(:,:) = sshn(:,:)                                       ! initial ssh
           DO jk = 1, jpk
@@ -302,24 +317,36 @@
              sc_loc_ini(:,:,jk) = tsn(:,:,jk,jp_sal) * fse3t_n(:,:,jk)   ! initial salt content
           END DO
-          hcssh_loc_ini(:,:) = tsn(:,:,1,jp_tem) * sshn(:,:)   ! initial heat content in ssh
-          scssh_loc_ini(:,:) = tsn(:,:,1,jp_sal) * sshn(:,:)   ! initial salt content in ssh
-          frc_v = 0._wp                                           
-          frc_t = 0._wp                                           
-          frc_s = 0._wp                                                  
+          frc_v = 0.d0                                           ! volume       trend due to forcing
+          frc_t = 0.d0                                           ! heat content   -    -   -    -   
+          frc_s = 0.d0                                           ! salt content   -    -   -    -        
+          IF( .NOT. lk_vvl ) THEN
+             ssh_hc_loc_ini(:,:) = tsn(:,:,1,jp_tem) * sshn(:,:)   ! initial heat content in ssh
+             ssh_sc_loc_ini(:,:) = tsn(:,:,1,jp_sal) * sshn(:,:)   ! initial salt content in ssh
+             frc_wn_t = 0.d0                                       ! initial heat content misfit due to free surface
+             frc_wn_s = 0.d0                                       ! initial salt content misfit due to free surface
+          ENDIF
        ENDIF
 
      ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN   ! Create restart file
         !                                   ! -------------------
-        IF(lwp) WRITE(numout,*) '---- dia-rst ----'
+        IF(lwp) WRITE(numout,*) '~~~~~~~'
+        IF(lwp) WRITE(numout,*) ' dia_hsb_rst at it= ', kt,' date= ', ndastp
+        IF(lwp) WRITE(numout,*) '~~~~~~~'
+
         CALL iom_rstput( kt, nitrst, numrow, 'frc_v'   , frc_v     )
         CALL iom_rstput( kt, nitrst, numrow, 'frc_t'   , frc_t     )
         CALL iom_rstput( kt, nitrst, numrow, 'frc_s'   , frc_s     )
-        
+        IF( .NOT. lk_vvl ) THEN
+           CALL iom_rstput( kt, nitrst, numrow, 'frc_wn_t', frc_wn_t )
+           CALL iom_rstput( kt, nitrst, numrow, 'frc_wn_s', frc_wn_s )
+        ENDIF
         CALL iom_rstput( kt, nitrst, numrow, 'ssh_ini', ssh_ini )
         CALL iom_rstput( kt, nitrst, numrow, 'e3t_ini', e3t_ini )
         CALL iom_rstput( kt, nitrst, numrow, 'hc_loc_ini', hc_loc_ini )
         CALL iom_rstput( kt, nitrst, numrow, 'sc_loc_ini', sc_loc_ini )
-        CALL iom_rstput( kt, nitrst, numrow, 'hcssh_loc_ini', hcssh_loc_ini )
-        CALL iom_rstput( kt, nitrst, numrow, 'scssh_loc_ini', scssh_loc_ini )
+        IF( .NOT. lk_vvl ) THEN
+           CALL iom_rstput( kt, nitrst, numrow, 'ssh_hc_loc_ini', ssh_hc_loc_ini )
+           CALL iom_rstput( kt, nitrst, numrow, 'ssh_sc_loc_ini', ssh_sc_loc_ini )
+        ENDIF
         !
      ENDIF
