Index: trunk/NEMOGCM/NEMO/OPA_SRC/DIA/diahsb.F90
===================================================================
--- trunk/NEMOGCM/NEMO/OPA_SRC/DIA/diahsb.F90	(revision 4624)
+++ trunk/NEMOGCM/NEMO/OPA_SRC/DIA/diahsb.F90	(revision 4990)
@@ -9,21 +9,27 @@
 
    !!----------------------------------------------------------------------
+   !!   dia_hsb       : Diagnose the conservation of ocean heat and salt contents, and volume
+   !!   dia_hsb_rst   : Read or write DIA file in restart file
+   !!   dia_hsb_init  : Initialization of the conservation diagnostic
+   !!----------------------------------------------------------------------
    USE oce             ! ocean dynamics and tracers
    USE dom_oce         ! ocean space and time domain
    USE phycst          ! physical constants
    USE sbc_oce         ! surface thermohaline fluxes
-   USE in_out_manager  ! I/O manager
+   USE sbcrnf          ! river runoff
+   USE sbcisf          ! ice shelves
    USE domvvl          ! vertical scale factors
    USE traqsr          ! penetrative solar radiation
    USE trabbc          ! bottom boundary condition 
-   USE lib_mpp         ! distributed memory computing library
    USE trabbc          ! bottom boundary condition
    USE bdy_par         ! (for lk_bdy)
+   USE restart         ! ocean restart
+   !
+   USE iom             ! I/O manager
+   USE in_out_manager  ! I/O manager
+   USE lib_fortran     ! glob_sum
+   USE lib_mpp         ! distributed memory computing library
    USE timing          ! preformance summary
-   USE iom             ! I/O manager
-   USE lib_fortran     ! glob_sum
-   USE restart         ! ocean restart
-   USE wrk_nemo         ! work arrays
-   USE sbcrnf         ! river runoffd
+   USE wrk_nemo        ! work arrays
 
    IMPLICIT NONE
@@ -36,15 +42,15 @@
    LOGICAL, PUBLIC ::   ln_diahsb   !: check the heat and salt budgets
 
-   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
+   REAL(wp) ::   surf_tot              ! ocean surface
+   REAL(wp) ::   frc_t, frc_s, frc_v   ! global forcing trends
+   REAL(wp) ::   frc_wn_t, frc_wn_s    ! global forcing trends
+   !
+   REAL(wp), DIMENSION(:,:)  , ALLOCATABLE ::   surf          , ssh_ini          !
+   REAL(wp), DIMENSION(:,:)  , ALLOCATABLE ::   ssh_hc_loc_ini, ssh_sc_loc_ini   !
+   REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::   hc_loc_ini, sc_loc_ini, e3t_ini  !
 
    !! * Substitutions
 #  include "domzgr_substitute.h90"
 #  include "vectopt_loop_substitute.h90"
-
    !!----------------------------------------------------------------------
    !! NEMO/OPA 3.3 , NEMO Consortium (2010)
@@ -52,5 +58,4 @@
    !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
    !!----------------------------------------------------------------------
-
 CONTAINS
 
@@ -67,27 +72,37 @@
       !!---------------------------------------------------------------------------
       INTEGER, INTENT(in) ::   kt   ! ocean time-step index
-      !!
-      INTEGER    ::   jk                          ! dummy loop indice
-      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   !    -     -
+      !
+      INTEGER    ::   ji, jj, jk                  ! dummy loop indice
+      REAL(wp)   ::   zdiff_hc    , zdiff_sc      ! heat and salt content variations
+      REAL(wp)   ::   zdiff_hc1   , zdiff_sc1     !  -         -     -        - 
+      REAL(wp)   ::   zdiff_v1    , zdiff_v2      ! volume variation
+      REAL(wp)   ::   zerr_hc1    , zerr_sc1      ! heat and salt content misfit
+      REAL(wp)   ::   zvol_tot                    ! volume
+      REAL(wp)   ::   z_frc_trd_t , z_frc_trd_s   !    -     -
+      REAL(wp)   ::   z_frc_trd_v                 !    -     -
+      REAL(wp)   ::   z_wn_trd_t , z_wn_trd_s     !    -     -
+      REAL(wp)   ::   z_ssh_hc , z_ssh_sc         !    -     -
+      REAL(wp), DIMENSION(:,:), POINTER ::   z2d0, z2d1
       !!---------------------------------------------------------------------------
       IF( nn_timing == 1 )   CALL timing_start('dia_hsb')      
-
+      CALL wrk_alloc( jpi,jpj,   z2d0, z2d1 )
+      !
+      tsn(:,:,:,1) = tsn(:,:,:,1) * tmask(:,:,:) ; tsb(:,:,:,1) = tsb(:,:,:,1) * tmask(:,:,:) ;
+      tsn(:,:,:,2) = tsn(:,:,:,2) * tmask(:,:,:) ; tsb(:,:,:,2) = tsb(:,:,:,2) * tmask(:,:,:) ;
       ! ------------------------- !
       ! 1 - Trends due to forcing !
       ! ------------------------- !
-      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
+      z_frc_trd_v = r1_rau0 * glob_sum( - ( emp(:,:) - rnf(:,:) + rdivisf * fwfisf(:,:) ) * 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 geothermal ice shelf
+      IF( nn_isf .GE. 1 )  THEN
+          z_frc_trd_t = z_frc_trd_t &
+              &   + glob_sum( ( risf_tsc(:,:,jp_tem) - rdivisf * fwfisf(:,:) * (-1.9) * r1_rau0 ) * surf(:,:) )
+          z_frc_trd_s = z_frc_trd_s + (1.0_wp - rdivisf) * glob_sum( risf_tsc(:,:,jp_sal) * surf(:,:) )
+      ENDIF
 
       ! Add penetrative solar radiation
@@ -97,6 +112,13 @@
       !
       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) )
+         z2d0=0.0_wp ; z2d1=0.0_wp
+         DO ji=1,jpi
+            DO jj=1,jpj
+              z2d0(ji,jj) = surf(ji,jj) * wn(ji,jj,mikt(ji,jj)) * tsb(ji,jj,mikt(ji,jj),jp_tem)
+              z2d1(ji,jj) = surf(ji,jj) * wn(ji,jj,mikt(ji,jj)) * tsb(ji,jj,mikt(ji,jj),jp_sal)
+            ENDDO
+         ENDDO
+         z_wn_trd_t = - glob_sum( z2d0 ) 
+         z_wn_trd_s = - glob_sum( z2d1 )
       ENDIF
 
@@ -113,7 +135,7 @@
       ! 2 -  Content variations !
       ! ------------------------ !
-      zdiff_v2 = 0.d0
-      zdiff_hc = 0.d0
-      zdiff_sc = 0.d0
+      zdiff_v2 = 0._wp
+      zdiff_hc = 0._wp
+      zdiff_sc = 0._wp
 
       ! volume variation (calculated with ssh)
@@ -122,6 +144,13 @@
       ! 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(:,:) ) )
+         z2d0 = 0._wp   ;   z2d1 = 0._wp
+         DO ji = 1, jpi
+            DO jj = 1, jpj
+              z2d0(ji,jj) = surf(ji,jj) * ( tsn(ji,jj,mikt(ji,jj),jp_tem) * sshn(ji,jj) - ssh_hc_loc_ini(ji,jj) ) 
+              z2d1(ji,jj) = surf(ji,jj) * ( tsn(ji,jj,mikt(ji,jj),jp_sal) * sshn(ji,jj) - ssh_sc_loc_ini(ji,jj) ) 
+            END DO
+         END DO
+         z_ssh_hc = glob_sum( z2d0 ) 
+         z_ssh_sc = glob_sum( z2d1 ) 
       ENDIF
 
@@ -153,8 +182,15 @@
       ! 3 - Diagnostics writing !
       ! ----------------------- !
-      zvol_tot   = 0.d0                                                   ! total ocean volume
+      zvol_tot = 0._wp                    ! total ocean volume (calculated with scale factors)
       DO jk = 1, jpkm1
          zvol_tot  = zvol_tot + glob_sum( surf(:,:) * tmask(:,:,jk) * fse3t_n(:,:,jk) )
       END DO
+
+!!gm to be added ?
+!      IF( .NOT. lk_vvl ) THEN            ! fixed volume, add the ssh contribution
+!        zvol_tot = zvol_tot + glob_sum( surf(:,:) * sshn(:,:) )
+!      ENDIF
+!!gm end
+
 
       IF( lk_vvl ) THEN
@@ -183,90 +219,10 @@
       IF( lrst_oce )   CALL dia_hsb_rst( kt, 'WRITE' )
 
+      CALL wrk_dealloc( jpi,jpj,   z2d0, z2d1 )
+
       IF( nn_timing == 1 )   CALL timing_stop('dia_hsb')
-!
+      !
    END SUBROUTINE dia_hsb
 
-
-   SUBROUTINE dia_hsb_init
-      !!---------------------------------------------------------------------------
-      !!                  ***  ROUTINE dia_hsb  ***
-      !!     
-      !! ** Purpose: Initialization for the heat salt volume budgets
-      !!	
-      !! ** Method : Compute initial heat content, salt content and volume
-      !!
-      !! ** Action : - Compute initial heat content, salt content and volume
-      !!             - Initialize forcing trends
-      !!             - Compute coefficients for conversion
-      !!---------------------------------------------------------------------------
-      INTEGER            ::   jk       ! dummy loop indice
-      INTEGER            ::   ierror   ! local integer
-      !!
-      NAMELIST/namhsb/ ln_diahsb
-      !
-      INTEGER  ::   ios
-      !!----------------------------------------------------------------------
-
-      IF(lwp) THEN
-         WRITE(numout,*)
-         WRITE(numout,*) 'dia_hsb_init : check the heat and salt budgets'
-         WRITE(numout,*) '~~~~~~~~ '
-      ENDIF
-
-      REWIND( numnam_ref )              ! Namelist namhsb in reference namelist
-      READ  ( numnam_ref, namhsb, IOSTAT = ios, ERR = 901)
-901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namhsb in reference namelist', lwp )
-
-      REWIND( numnam_cfg )              ! Namelist namhsb in configuration namelist
-      READ  ( numnam_cfg, namhsb, IOSTAT = ios, ERR = 902 )
-902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namhsb in configuration namelist', lwp )
-      IF(lwm) WRITE ( numond, namhsb )
-
-      !
-      IF(lwp) THEN                   ! Control print
-         WRITE(numout,*)
-         WRITE(numout,*) 'dia_hsb_init : check the heat and salt budgets'
-         WRITE(numout,*) '~~~~~~~~~~~~'
-         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
-!      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
 
    SUBROUTINE dia_hsb_rst( kt, cdrw )
@@ -281,6 +237,6 @@
      CHARACTER(len=*), INTENT(in) ::   cdrw   ! "READ"/"WRITE" flag
      !
-     INTEGER ::   jk   ! 
-     INTEGER ::   id1   ! local integers
+     INTEGER ::   ji, jj, jk   ! dummy loop indices
+     INTEGER ::   id1          ! local integers
      !!----------------------------------------------------------------------
      !
@@ -317,12 +273,16 @@
              sc_loc_ini(:,:,jk) = tsn(:,:,jk,jp_sal) * fse3t_n(:,:,jk)   ! initial salt content
           END DO
-          frc_v = 0.d0                                           ! volume       trend due to forcing
-          frc_t = 0.d0                                           ! heat content   -    -   -    -   
-          frc_s = 0.d0                                           ! salt content   -    -   -    -        
+          frc_v = 0._wp                                           ! volume       trend due to forcing
+          frc_t = 0._wp                                           ! heat content   -    -   -    -   
+          frc_s = 0._wp                                           ! 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
+             DO ji=1,jpi
+                DO jj=1,jpj
+                   ssh_hc_loc_ini(ji,jj) = tsn(ji,jj,mikt(ji,jj),jp_tem) * sshn(ji,jj)   ! initial heat content in ssh
+                   ssh_sc_loc_ini(ji,jj) = tsn(ji,jj,mikt(ji,jj),jp_sal) * sshn(ji,jj)   ! initial salt content in ssh
+                ENDDO
+             ENDDO
+             frc_wn_t = 0._wp                                       ! initial heat content misfit due to free surface
+             frc_wn_s = 0._wp                                       ! initial salt content misfit due to free surface
           ENDIF
        ENDIF
@@ -354,4 +314,86 @@
    END SUBROUTINE dia_hsb_rst
 
+
+   SUBROUTINE dia_hsb_init
+      !!---------------------------------------------------------------------------
+      !!                  ***  ROUTINE dia_hsb  ***
+      !!     
+      !! ** Purpose: Initialization for the heat salt volume budgets
+      !!	
+      !! ** Method : Compute initial heat content, salt content and volume
+      !!
+      !! ** Action : - Compute initial heat content, salt content and volume
+      !!             - Initialize forcing trends
+      !!             - Compute coefficients for conversion
+      !!---------------------------------------------------------------------------
+      INTEGER ::   jk       ! dummy loop indice
+      INTEGER ::   ierror   ! local integer
+      INTEGER ::   ios
+      !
+      NAMELIST/namhsb/ ln_diahsb
+      !!----------------------------------------------------------------------
+
+      IF(lwp) THEN
+         WRITE(numout,*)
+         WRITE(numout,*) 'dia_hsb_init : check the heat and salt budgets'
+         WRITE(numout,*) '~~~~~~~~ '
+      ENDIF
+
+      REWIND( numnam_ref )              ! Namelist namhsb in reference namelist
+      READ  ( numnam_ref, namhsb, IOSTAT = ios, ERR = 901)
+901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namhsb in reference namelist', lwp )
+
+      REWIND( numnam_cfg )              ! Namelist namhsb in configuration namelist
+      READ  ( numnam_cfg, namhsb, IOSTAT = ios, ERR = 902 )
+902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namhsb in configuration namelist', lwp )
+      IF(lwm) WRITE ( numond, namhsb )
+
+      !
+      IF(lwp) THEN                   ! Control print
+         WRITE(numout,*)
+         WRITE(numout,*) 'dia_hsb_init : check the heat and salt budgets'
+         WRITE(numout,*) '~~~~~~~~~~~~'
+         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
+         !      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_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
+
    !!======================================================================
 END MODULE diahsb
