Index: branches/nemo_v3_3_beta/NEMOGCM/NEMO/LIM_SRC_3/limdia.F90
===================================================================
--- branches/nemo_v3_3_beta/NEMOGCM/NEMO/LIM_SRC_3/limdia.F90	(revision 2313)
+++ branches/nemo_v3_3_beta/NEMOGCM/NEMO/LIM_SRC_3/limdia.F90	(revision 2439)
@@ -2,73 +2,58 @@
    !!======================================================================
    !!                       ***  MODULE limdia   ***
-   !!                      diagnostics of ice model 
+   !!  LIM-3 sea ice model :   diagnostics of ice model 
    !!======================================================================
-   !! To add a new field
-   !! 1) in lim_dia : add its definition for both hemispheres if wished
-   !! 2) add the new titles in lim_dia_init
+   !! History :  3.2  ! 2007-01  (M. Vancoppenolle)  Code adapted from LIM-2
+   !!             -   ! 2008-03  (M. Vancoppenolle)  add lim_dia_init
    !!----------------------------------------------------------------------
 #if defined key_lim3
    !!----------------------------------------------------------------------
-   !!   'key_lim3' :                                   LIM3 sea-ice model
-   !!----------------------------------------------------------------------
-   !!   lim_dia      : computation of the time evolution of keys var.
-   !!   lim_dia_init : initialization and namelist read
-   !!----------------------------------------------------------------------
-   !! * Modules used
-   USE phycst
-   USE in_out_manager
-   USE par_ice         ! ice parameters
-   USE sbc_ice         ! ice variables
-   USE daymod
-   USE dom_ice
-   USE ice
-   USE iceini
-   USE limistate
-   USE dom_oce
-   USE sbc_oce         ! Surface boundary condition: ocean fields
-   USE dom_oce
-   USE lib_mpp
-
+   !!   'key_lim3'                                       LIM3 sea-ice model
+   !!----------------------------------------------------------------------
+   !!   lim_dia        : computation and output of the time evolution of keys variables
+   !!   lim_dia_init   : initialization and namelist read
+   !!----------------------------------------------------------------------
+   USE ice             ! LIM-3: sea-ice variable
+   USE par_ice         ! LIM-3: ice parameters
+!   USE ice_oce         ! LIM-3: ice-ocean share variables
+   USE iceini          ! LIM-3: sea-ice initialization
+   USE limistate       ! LIM-3: sea-ice initial state
+   USE dom_ice         ! LIM-3: sea-ice domain
+   USE dom_oce         ! ocean domain
+   USE sbc_oce         ! surface boundary condition: ocean fields
+   USE daymod          ! model calendar
+   USE phycst          ! physical constant
+   USE in_out_manager  ! I/O manager
+   USE lib_mpp         ! MPP library
+   
    IMPLICIT NONE
    PRIVATE
 
-   !! * Routine accessibility
    PUBLIC lim_dia       ! called by ice_step
 
-   !! * Shared module variables
-   INTEGER, PUBLIC  ::  &
-      ntmoy   = 1 ,     &  !: instantaneous values of ice evolution or averaging ntmoy
-      ninfo   = 1          !: frequency of ouputs on file ice_evolu in case of averaging
-
-   !! * Module variables
-   INTEGER, PARAMETER ::   &  ! Parameters for outputs to files "evolu"
-      jpinfmx = 100         ,    &  ! maximum number of key variables
-      jpchinf = 5           ,    &  ! ???
-      jpchsep = jpchinf + 2         ! ???
-
-   INTEGER          ::  &
-      nfrinf  = 4    ,  &  ! number of variables written in one line 
-      nferme ,          &  ! last time step at which the var. are written on file
-      nvinfo ,          &  ! number of total variables 
-      nbvt   ,          &  ! number of time variables
-      naveg                ! number of step for accumulation before averaging
-
-   CHARACTER(len=8) ::  &
-      fmtinf  = '1PE13.5 ' ! format of the output values  
-
-   CHARACTER(len=30) :: &
-      fmtw  ,           &  ! formats
-      fmtr  ,           &  ! ???
-      fmtitr               ! ???
-
-   CHARACTER(len=jpchsep), DIMENSION(jpinfmx) ::   &
-      titvar               ! title of key variables
-
-   REAL(wp), DIMENSION(jpinfmx) ::  &
-      vinfom               ! temporary working space
-   REAL(wp) :: &
-      epsi06 = 1.e-06
-   REAL(wp), DIMENSION(jpi,jpj) ::   &
-      aire                 ! masked grid cell area
+   INTEGER, PUBLIC  ::   ntmoy   = 1   !: instantaneous values of ice evolution or averaging ntmoy
+   INTEGER, PUBLIC  ::   ninfo   = 1   !: frequency of ouputs on file ice_evolu in case of averaging
+
+   !                                              !!! Parameters for outputs to files "evolu"
+   INTEGER, PARAMETER ::   jpinfmx = 100           ! maximum number of key variables
+   INTEGER, PARAMETER ::   jpchinf = 5             ! ???
+   INTEGER, PARAMETER ::   jpchsep = jpchinf + 2   ! ???
+
+   INTEGER  ::   nfrinf  = 4         ! number of variables written in one line 
+   INTEGER  ::   nferme              ! last time step at which the var. are written on file
+   INTEGER  ::   nvinfo              ! number of total variables 
+   INTEGER  ::   nbvt                ! number of time variables
+   INTEGER  ::   naveg               ! number of step for accumulation before averaging
+   REAL(wp) ::   epsi06 = 1.e-6_wp   ! small number
+
+   CHARACTER(len= 8) ::   fmtinf = '1PE13.5 '   ! format of the output values  
+   CHARACTER(len=30) ::   fmtw                  ! formats
+   CHARACTER(len=30) ::   fmtr                  ! ???
+   CHARACTER(len=30) ::   fmtitr                ! ???
+
+   CHARACTER(len=jpchsep), DIMENSION(jpinfmx) ::   titvar   ! title of key variables
+
+   REAL(wp), DIMENSION(jpinfmx) ::   vinfom     ! temporary working space
+   REAL(wp), DIMENSION(jpi,jpj) ::   aire       ! masked grid cell area
 
    !! * Substitutions
@@ -77,8 +62,6 @@
    !! NEMO/LIM3 3.3 , UCL - NEMO Consortium (2010)
    !! $Id$
-   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
-   !!----------------------------------------------------------------------
-
-
+   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
+   !!----------------------------------------------------------------------
 CONTAINS
 
@@ -89,42 +72,29 @@
       !! ** Purpose : Computation and outputs on file ice.evolu 
       !!      the temporal evolution of some key variables
-      !!
-      !! History :
-      !!   8.0  !  97-06  (Louvain-La-Neuve)  Original code
-      !!   8.5  !  02-09  (C. Ethe , G. Madec )  F90: Free form and module
-      !!   9.0  !  05-06  (M. Vancoppenolle LLN) Rehabilitation and fixing a few bugs
-      !!           01-07  (M. Vancoppenolle LLN) added new fields + fram strait
-      !!                                         export
-      !!-------------------------------------------------------------------
-      !! * Local variables
-      INTEGER  ::   jv,ji,jj,jl ! dummy loop indices
-      REAL(wp), DIMENSION(jpinfmx) ::  & 
-         vinfor           ! temporary working space 
-      REAL(wp) ::    &
-         zshift_date   , & ! date from the minimum ice extent
-         zday, zday_min, & ! current day, day of minimum extent
-         zafy, zamy,     & ! temporary area of fy and my ice
-         zindb
+      !!-------------------------------------------------------------------
+      INTEGER  ::   jv, ji, jj, jl   ! dummy loop indices
+      REAL(wp) ::   zshift_date      ! date from the minimum ice extent
+      REAL(wp) ::   zday, zday_min   ! current day, day of minimum extent
+      REAL(wp) ::   zafy, zamy       ! temporary area of fy and my ice
+      REAL(wp) ::   zindb
+      REAL(wp), DIMENSION(jpinfmx) ::   vinfor           ! temporary working space 
       !!-------------------------------------------------------------------
 
       ! 0) date from the minimum of ice extent
       !---------------------------------------
-      zday_min = 273.0        ! zday_min = date of minimum extent, here September 30th
-      zday = FLOAT(numit-nit000) * rdt_ice / ( 86400.0 * FLOAT(nn_fsbc) )
-      IF (zday.GT.zday_min) THEN 
-         zshift_date  =  zday - zday_min
-      ELSE
-         zshift_date  =  zday - (365.0 - zday_min)
+      zday_min = 273._wp        ! zday_min = date of minimum extent, here September 30th
+      zday = REAL(numit-nit000,wp) * rdt_ice / ( 86400._wp * REAL(nn_fsbc,wp) )
+      !
+      IF( zday > zday_min ) THEN   ;   zshift_date  =  zday - zday_min
+      ELSE                         ;   zshift_date  =  zday - (365.0 - zday_min)
       ENDIF
 
       IF( numit == nstart )   CALL lim_dia_init   ! initialisation of ice_evolu file      
 
-      ! temporal diagnostics 
-      vinfor(1) = REAL(numit)
+      vinfor(1) = REAL(numit)       ! time diagnostics 
       vinfor(2) = nyear
 
-      ! put everything to zero
-      DO jv = nbvt + 1, nvinfo
-         vinfor(jv) = 0.0
+      DO jv = nbvt + 1, nvinfo      ! put everything to zero
+         vinfor(jv) = 0._wp
       END DO
 
@@ -137,27 +107,27 @@
          DO ji = fs_2, fs_jpim1   ! vector opt.
             IF( tms(ji,jj) == 1 ) THEN
-               vinfor(3)  = vinfor(3)  + at_i(ji,jj)*aire(ji,jj) / 1.0e12 !ice area
-               IF (at_i(ji,jj).GT.0.15) vinfor(5) = vinfor(5) + aire(ji,jj) / 1.0e12 !ice extent
-               vinfor(7)  = vinfor(7)  + vt_i(ji,jj)*aire(ji,jj) / 1.0e12 !ice volume
-               vinfor(9)  = vinfor(9)  + vt_s(ji,jj)*aire(ji,jj) / 1.0e12 !snow volume
-               vinfor(15) = vinfor(15) + ot_i(ji,jj) *vt_i(ji,jj)*aire(ji,jj) / 1.0e12 !mean age
-               vinfor(29) = vinfor(29) + smt_i(ji,jj)*vt_i(ji,jj)*aire(ji,jj) / 1.0e12 !mean salinity
+               vinfor(3)  = vinfor(3)  + at_i(ji,jj)*aire(ji,jj) * 1.e-12_wp !ice area
+               IF (at_i(ji,jj).GT.0.15) vinfor(5) = vinfor(5) + aire(ji,jj) * 1.e-12_wp !ice extent
+               vinfor(7)  = vinfor(7)  + vt_i(ji,jj)*aire(ji,jj) * 1.e-12_wp !ice volume
+               vinfor(9)  = vinfor(9)  + vt_s(ji,jj)*aire(ji,jj) * 1.e-12_wp !snow volume
+               vinfor(15) = vinfor(15) + ot_i(ji,jj) *vt_i(ji,jj)*aire(ji,jj) * 1.e-12_wp !mean age
+               vinfor(29) = vinfor(29) + smt_i(ji,jj)*vt_i(ji,jj)*aire(ji,jj) * 1.e-12_wp !mean salinity
                ! the computation of this diagnostic is not reliable
                vinfor(31) = vinfor(31) + vt_i(ji,jj)*( u_ice(ji,jj)*u_ice(ji,jj) + & 
                   v_ice(ji,jj)*v_ice(ji,jj) )*aire(ji,jj)/1.0e12 
-               vinfor(53) = vinfor(53) + emps(ji,jj)*aire(ji,jj) / 1.0e12 !salt flux
-               vinfor(55) = vinfor(55) + fsbri(ji,jj)*aire(ji,jj) / 1.0e12 !brine drainage flux
-               vinfor(57) = vinfor(57) + fseqv(ji,jj)*aire(ji,jj) / 1.0e12 !equivalent salt flux
-               vinfor(59) = vinfor(59) +(sst_m(ji,jj)+rt0)*at_i(ji,jj)*aire(ji,jj) / 1.0e12  !SST
-               vinfor(61) = vinfor(61) + sss_m(ji,jj)*at_i(ji,jj)*aire(ji,jj) / 1.0e12  !SSS
-               vinfor(65) = vinfor(65) + et_s(ji,jj)/1.0e9*aire(ji,jj) / 1.0e12  ! snow temperature
-               vinfor(67) = vinfor(67) + et_i(ji,jj)/1.0e9*aire(ji,jj) / 1.0e12       ! ice heat content
-               vinfor(69) = vinfor(69) + v_i(ji,jj,1)*aire(ji,jj) / 1.0e12 !ice volume
-               vinfor(71) = vinfor(71) + v_i(ji,jj,2)*aire(ji,jj) / 1.0e12 !ice volume
-               vinfor(73) = vinfor(73) + v_i(ji,jj,3)*aire(ji,jj) / 1.0e12 !ice volume
-               vinfor(75) = vinfor(75) + v_i(ji,jj,4)*aire(ji,jj) / 1.0e12 !ice volume
-               vinfor(77) = vinfor(77) + v_i(ji,jj,5)*aire(ji,jj) / 1.0e12 !ice volume
+               vinfor(53) = vinfor(53) + emps(ji,jj)*aire(ji,jj) * 1.e-12_wp !salt flux
+               vinfor(55) = vinfor(55) + fsbri(ji,jj)*aire(ji,jj) * 1.e-12_wp !brine drainage flux
+               vinfor(57) = vinfor(57) + fseqv(ji,jj)*aire(ji,jj) * 1.e-12_wp !equivalent salt flux
+               vinfor(59) = vinfor(59) +(sst_m(ji,jj)+rt0)*at_i(ji,jj)*aire(ji,jj) * 1.e-12_wp  !SST
+               vinfor(61) = vinfor(61) + sss_m(ji,jj)*at_i(ji,jj)*aire(ji,jj) * 1.e-12_wp  !SSS
+               vinfor(65) = vinfor(65) + et_s(ji,jj)/1.0e9*aire(ji,jj) * 1.e-12_wp  ! snow temperature
+               vinfor(67) = vinfor(67) + et_i(ji,jj)/1.0e9*aire(ji,jj) * 1.e-12_wp       ! ice heat content
+               vinfor(69) = vinfor(69) + v_i(ji,jj,1)*aire(ji,jj) * 1.e-12_wp !ice volume
+               vinfor(71) = vinfor(71) + v_i(ji,jj,2)*aire(ji,jj) * 1.e-12_wp !ice volume
+               vinfor(73) = vinfor(73) + v_i(ji,jj,3)*aire(ji,jj) * 1.e-12_wp !ice volume
+               vinfor(75) = vinfor(75) + v_i(ji,jj,4)*aire(ji,jj) * 1.e-12_wp !ice volume
+               vinfor(77) = vinfor(77) + v_i(ji,jj,5)*aire(ji,jj) * 1.e-12_wp !ice volume
                vinfor(79) = 0.0
-               vinfor(81) = vinfor(81) + emp(ji,jj)*aire(ji,jj) / 1.0e12 ! mass flux
+               vinfor(81) = vinfor(81) + emp(ji,jj)*aire(ji,jj) * 1.e-12_wp ! mass flux
             ENDIF
          END DO
@@ -168,5 +138,5 @@
             DO ji = fs_2, fs_jpim1   ! vector opt.
                IF( tms(ji,jj) == 1 ) THEN
-                  vinfor(11) = vinfor(11) + v_i(ji,jj,jl)*aire(ji,jj) / 1.0e12 !undef def ice volume
+                  vinfor(11) = vinfor(11) + v_i(ji,jj,jl)*aire(ji,jj) * 1.e-12_wp !undef def ice volume
                ENDIF
             END DO
@@ -174,5 +144,5 @@
       END DO
 
-      vinfor(13) = 0.0
+      vinfor(13) = 0._wp
 
       vinfor(15) = vinfor(15) / MAX(vinfor(7),epsi06) ! these have to be divided by total ice volume to have the
@@ -198,6 +168,6 @@
             DO ji = fs_2, fs_jpim1   ! vector opt.
                IF( tms(ji,jj) == 1 ) THEN
-                  vinfor(33) = vinfor(33) + d_v_i_trp(ji,jj,jl)*aire(ji,jj) / 1.0e12 !ice volume
-                  vinfor(35) = vinfor(35) + d_v_i_thd(ji,jj,jl)*aire(ji,jj) / 1.0e12 !ice volume
+                  vinfor(33) = vinfor(33) + d_v_i_trp(ji,jj,jl)*aire(ji,jj) * 1.e-12_wp !ice volume
+                  vinfor(35) = vinfor(35) + d_v_i_thd(ji,jj,jl)*aire(ji,jj) * 1.e-12_wp !ice volume
                ENDIF
             END DO
@@ -208,10 +178,10 @@
          DO ji = fs_2, fs_jpim1   ! vector opt.
             IF( tms(ji,jj) == 1 ) THEN
-               vinfor(37) = vinfor(37) + diag_sni_gr(ji,jj)*aire(ji,jj) / 1.0e12 !th growth rates
-               vinfor(39) = vinfor(39) + diag_lat_gr(ji,jj)*aire(ji,jj) / 1.0e12 
-               vinfor(41) = vinfor(41) + diag_bot_gr(ji,jj)*aire(ji,jj) / 1.0e12
-               vinfor(43) = vinfor(43) + diag_dyn_gr(ji,jj)*aire(ji,jj) / 1.0e12 
-               vinfor(45) = vinfor(45) + dv_dt_thd(ji,jj,5)*aire(ji,jj) / 1.0e12
-               vinfor(47) = vinfor(47) + v_newice(ji,jj) *aire(ji,jj) / 1.0e12 / rdt_ice ! volume acc in OW
+               vinfor(37) = vinfor(37) + diag_sni_gr(ji,jj)*aire(ji,jj) * 1.e-12_wp !th growth rates
+               vinfor(39) = vinfor(39) + diag_lat_gr(ji,jj)*aire(ji,jj) * 1.e-12_wp 
+               vinfor(41) = vinfor(41) + diag_bot_gr(ji,jj)*aire(ji,jj) * 1.e-12_wp
+               vinfor(43) = vinfor(43) + diag_dyn_gr(ji,jj)*aire(ji,jj) * 1.e-12_wp 
+               vinfor(45) = vinfor(45) + dv_dt_thd(ji,jj,5)*aire(ji,jj) * 1.e-12_wp
+               vinfor(47) = vinfor(47) + v_newice(ji,jj) *aire(ji,jj) * 1.e-12_wp / rdt_ice ! volume acc in OW
             ENDIF
          END DO
@@ -222,5 +192,5 @@
             DO ji = fs_2, fs_jpim1   ! vector opt.
                IF( tms(ji,jj) == 1 ) THEN
-                  vinfor(63) = vinfor(63) + t_su(ji,jj,jl)*a_i(ji,jj,jl)*aire(ji,jj) / 1.0e12
+                  vinfor(63) = vinfor(63) + t_su(ji,jj,jl)*a_i(ji,jj,jl)*aire(ji,jj) * 1.e-12_wp
                ENDIF
             END DO
@@ -238,21 +208,21 @@
                DO jl = 1, jpl
                   IF ((o_i(ji,jj,jl) - zshift_date).LT.0.0) THEN
-                     vinfor(17) = vinfor(17) + a_i(ji,jj,jl)*aire(ji,jj) / 1.0e12 ! FY ice area
-                     vinfor(25) = vinfor(25) + v_i(ji,jj,jl)*aire(ji,jj) / 1.0e12 ! FY ice volume
-                     vinfor(49) = vinfor(49) + sm_i(ji,jj,jl)*v_i(ji,jj,jl)*aire(ji,jj) / 1.0e12 !FY ice salinity
+                     vinfor(17) = vinfor(17) + a_i(ji,jj,jl)*aire(ji,jj) * 1.e-12_wp ! FY ice area
+                     vinfor(25) = vinfor(25) + v_i(ji,jj,jl)*aire(ji,jj) * 1.e-12_wp ! FY ice volume
+                     vinfor(49) = vinfor(49) + sm_i(ji,jj,jl)*v_i(ji,jj,jl)*aire(ji,jj) * 1.e-12_wp !FY ice salinity
                      zafy = zafy + a_i(ji,jj,jl)
                   ENDIF
                   IF ((o_i(ji,jj,jl) - zshift_date).GT.0.0) THEN
-                     vinfor(19) = vinfor(19) + a_i(ji,jj,jl)*aire(ji,jj) / 1.0e12    ! MY ice area
-                     vinfor(27) = vinfor(27) + v_i(ji,jj,jl)*aire(ji,jj) / 1.0e12 ! MY ice volume
-                     vinfor(51) = vinfor(51) + sm_i(ji,jj,jl)*v_i(ji,jj,jl)*aire(ji,jj) / 1.0e12 !MY ice salinity
+                     vinfor(19) = vinfor(19) + a_i(ji,jj,jl)*aire(ji,jj) * 1.e-12_wp    ! MY ice area
+                     vinfor(27) = vinfor(27) + v_i(ji,jj,jl)*aire(ji,jj) * 1.e-12_wp ! MY ice volume
+                     vinfor(51) = vinfor(51) + sm_i(ji,jj,jl)*v_i(ji,jj,jl)*aire(ji,jj) * 1.e-12_wp !MY ice salinity
                      zamy = zamy + a_i(ji,jj,jl)
                   ENDIF
                END DO
                IF ((at_i(ji,jj).GT.0.15).AND.(zafy.GT.zamy)) THEN
-                  vinfor(21) = vinfor(21) + aire(ji,jj) / 1.0e12 ! Seasonal ice extent
+                  vinfor(21) = vinfor(21) + aire(ji,jj) * 1.e-12_wp ! Seasonal ice extent
                ENDIF
                IF ((at_i(ji,jj).GT.0.15).AND.(zafy.LE.zamy)) THEN
-                  vinfor(23) = vinfor(23) + aire(ji,jj) / 1.0e12 ! Perennial ice extent
+                  vinfor(23) = vinfor(23) + aire(ji,jj) * 1.e-12_wp ! Perennial ice extent
                ENDIF
             ENDIF
@@ -274,7 +244,7 @@
       DO ji = 134, 138
          vinfor(83) = vinfor(83) - v_ice(ji,jj) * & 
-            e1t(ji,jj)*at_i(ji,jj)*rdt_ice / 1.0e12
+            e1t(ji,jj)*at_i(ji,jj)*rdt_ice * 1.e-12_wp
          vinfor(84) = vinfor(84) - v_ice(ji,jj) * & 
-            e1t(ji,jj)*vt_i(ji,jj)*rdt_ice / 1.0e12
+            e1t(ji,jj)*vt_i(ji,jj)*rdt_ice * 1.e-12_wp
       END DO
 
@@ -287,28 +257,28 @@
          DO ji = fs_2, fs_jpim1   ! vector opt.
             IF( tms(ji,jj) == 1 ) THEN
-               vinfor(4)  = vinfor(4)  + at_i(ji,jj)*aire(ji,jj) / 1.0e12 !ice area
-               IF (at_i(ji,jj).GT.0.15) vinfor(6) = vinfor(6) + aire(ji,jj) / 1.0e12 !ice extent
-               vinfor(8)  = vinfor(8)  + vt_i(ji,jj)*aire(ji,jj) / 1.0e12 !ice volume
-               vinfor(10) = vinfor(10) + vt_s(ji,jj)*aire(ji,jj) / 1.0e12 !snow volume
-               vinfor(16) = vinfor(16) + ot_i(ji,jj)*vt_i(ji,jj)*aire(ji,jj) / 1.0e12 !mean age
-               vinfor(30) = vinfor(30) + smt_i(ji,jj)*vt_i(ji,jj)*aire(ji,jj) / 1.0e12 !mean salinity
+               vinfor(4)  = vinfor(4)  + at_i(ji,jj)*aire(ji,jj) * 1.e-12_wp !ice area
+               IF (at_i(ji,jj).GT.0.15) vinfor(6) = vinfor(6) + aire(ji,jj) * 1.e-12_wp !ice extent
+               vinfor(8)  = vinfor(8)  + vt_i(ji,jj)*aire(ji,jj) * 1.e-12_wp !ice volume
+               vinfor(10) = vinfor(10) + vt_s(ji,jj)*aire(ji,jj) * 1.e-12_wp !snow volume
+               vinfor(16) = vinfor(16) + ot_i(ji,jj)*vt_i(ji,jj)*aire(ji,jj) * 1.e-12_wp !mean age
+               vinfor(30) = vinfor(30) + smt_i(ji,jj)*vt_i(ji,jj)*aire(ji,jj) * 1.e-12_wp !mean salinity
                ! this diagnostic is not well computed (weighted by vol instead
                ! of area)
                vinfor(32) = vinfor(32) + vt_i(ji,jj)*( u_ice(ji,jj)*u_ice(ji,jj) + & 
                   v_ice(ji,jj)*v_ice(ji,jj) )*aire(ji,jj)/1.0e12 !ice vel
-               vinfor(54) = vinfor(54) + at_i(ji,jj)*emps(ji,jj)*aire(ji,jj) / 1.0e12 ! Total salt flux
-               vinfor(56) = vinfor(56) + at_i(ji,jj)*fsbri(ji,jj)*aire(ji,jj) / 1.0e12 ! Brine drainage salt flux
-               vinfor(58) = vinfor(58) + at_i(ji,jj)*fseqv(ji,jj)*aire(ji,jj) / 1.0e12 ! Equivalent salt flux
-               vinfor(60) = vinfor(60) +(sst_m(ji,jj)+rt0)*at_i(ji,jj)*aire(ji,jj) / 1.0e12  !SST
-               vinfor(62) = vinfor(62) + sss_m(ji,jj)*at_i(ji,jj)*aire(ji,jj) / 1.0e12  !SSS
-               vinfor(66) = vinfor(66) + et_s(ji,jj)/1.0e9*aire(ji,jj) / 1.0e12 ! snow temperature
-               vinfor(68) = vinfor(68) + et_i(ji,jj)/1.0e9*aire(ji,jj) / 1.0e12 ! ice enthalpy
-               vinfor(70) = vinfor(70) + v_i(ji,jj,1)*aire(ji,jj) / 1.0e12 !ice volume
-               vinfor(72) = vinfor(72) + v_i(ji,jj,2)*aire(ji,jj) / 1.0e12 !ice volume
-               vinfor(74) = vinfor(74) + v_i(ji,jj,3)*aire(ji,jj) / 1.0e12 !ice volume
-               vinfor(76) = vinfor(76) + v_i(ji,jj,4)*aire(ji,jj) / 1.0e12 !ice volume
-               vinfor(78) = vinfor(78) + v_i(ji,jj,5)*aire(ji,jj) / 1.0e12 !ice volume
+               vinfor(54) = vinfor(54) + at_i(ji,jj)*emps(ji,jj)*aire(ji,jj) * 1.e-12_wp ! Total salt flux
+               vinfor(56) = vinfor(56) + at_i(ji,jj)*fsbri(ji,jj)*aire(ji,jj) * 1.e-12_wp ! Brine drainage salt flux
+               vinfor(58) = vinfor(58) + at_i(ji,jj)*fseqv(ji,jj)*aire(ji,jj) * 1.e-12_wp ! Equivalent salt flux
+               vinfor(60) = vinfor(60) +(sst_m(ji,jj)+rt0)*at_i(ji,jj)*aire(ji,jj) * 1.e-12_wp  !SST
+               vinfor(62) = vinfor(62) + sss_m(ji,jj)*at_i(ji,jj)*aire(ji,jj) * 1.e-12_wp  !SSS
+               vinfor(66) = vinfor(66) + et_s(ji,jj)/1.0e9*aire(ji,jj) * 1.e-12_wp ! snow temperature
+               vinfor(68) = vinfor(68) + et_i(ji,jj)/1.0e9*aire(ji,jj) * 1.e-12_wp ! ice enthalpy
+               vinfor(70) = vinfor(70) + v_i(ji,jj,1)*aire(ji,jj) * 1.e-12_wp !ice volume
+               vinfor(72) = vinfor(72) + v_i(ji,jj,2)*aire(ji,jj) * 1.e-12_wp !ice volume
+               vinfor(74) = vinfor(74) + v_i(ji,jj,3)*aire(ji,jj) * 1.e-12_wp !ice volume
+               vinfor(76) = vinfor(76) + v_i(ji,jj,4)*aire(ji,jj) * 1.e-12_wp !ice volume
+               vinfor(78) = vinfor(78) + v_i(ji,jj,5)*aire(ji,jj) * 1.e-12_wp !ice volume
                vinfor(80) = 0.0
-               vinfor(82) = vinfor(82) + emp(ji,jj)*aire(ji,jj) / 1.0e12 ! mass flux
+               vinfor(82) = vinfor(82) + emp(ji,jj)*aire(ji,jj) * 1.e-12_wp ! mass flux
             ENDIF
          END DO
@@ -318,5 +288,5 @@
          DO jj = 2, njeqm1
             DO ji = fs_2, fs_jpim1   ! vector opt.
-               vinfor(12) = vinfor(12) + v_i(ji,jj,jl)*aire(ji,jj) / 1.0e12 !undef def ice volume
+               vinfor(12) = vinfor(12) + v_i(ji,jj,jl)*aire(ji,jj) * 1.e-12_wp !undef def ice volume
             END DO
          END DO
@@ -349,6 +319,6 @@
             DO ji = fs_2, fs_jpim1   ! vector opt.
                IF( tms(ji,jj) == 1 ) THEN
-                  vinfor(34) = vinfor(34) + d_v_i_trp(ji,jj,jl)*aire(ji,jj) / 1.0e12 !ice volume
-                  vinfor(36) = vinfor(36) + d_v_i_thd(ji,jj,jl)*aire(ji,jj) / 1.0e12 !ice volume
+                  vinfor(34) = vinfor(34) + d_v_i_trp(ji,jj,jl)*aire(ji,jj) * 1.e-12_wp !ice volume
+                  vinfor(36) = vinfor(36) + d_v_i_thd(ji,jj,jl)*aire(ji,jj) * 1.e-12_wp !ice volume
                ENDIF
             END DO
@@ -359,14 +329,13 @@
          DO ji = fs_2, fs_jpim1   ! vector opt.
             IF( tms(ji,jj) == 1 ) THEN
-               vinfor(38) = vinfor(38) + diag_sni_gr(ji,jj)*aire(ji,jj) / 1.0e12 !th growth rates
-               vinfor(40) = vinfor(40) + diag_lat_gr(ji,jj)*aire(ji,jj) / 1.0e12 
-               vinfor(42) = vinfor(42) + diag_bot_gr(ji,jj)*aire(ji,jj) / 1.0e12
-               vinfor(44) = vinfor(44) + diag_dyn_gr(ji,jj)*aire(ji,jj) / 1.0e12 
-               vinfor(46) = vinfor(46) + dv_dt_thd(ji,jj,5)*aire(ji,jj) / 1.0e12
-               vinfor(48) = vinfor(48) + v_newice(ji,jj) *aire(ji,jj) / 1.0e12 / rdt_ice ! volume acc in OW
+               vinfor(38) = vinfor(38) + diag_sni_gr(ji,jj)*aire(ji,jj) * 1.e-12_wp !th growth rates
+               vinfor(40) = vinfor(40) + diag_lat_gr(ji,jj)*aire(ji,jj) * 1.e-12_wp 
+               vinfor(42) = vinfor(42) + diag_bot_gr(ji,jj)*aire(ji,jj) * 1.e-12_wp
+               vinfor(44) = vinfor(44) + diag_dyn_gr(ji,jj)*aire(ji,jj) * 1.e-12_wp 
+               vinfor(46) = vinfor(46) + dv_dt_thd(ji,jj,5)*aire(ji,jj) * 1.e-12_wp
+               vinfor(48) = vinfor(48) + v_newice(ji,jj) *aire(ji,jj) * 1.e-12_wp / rdt_ice ! volume acc in OW
             ENDIF
          END DO
       END DO
-
 
       DO jl = 1, jpl
@@ -374,5 +343,5 @@
             DO ji = fs_2, fs_jpim1   ! vector opt.
                IF( tms(ji,jj) == 1 ) THEN
-                  vinfor(64) = vinfor(64) + t_su(ji,jj,jl)*a_i(ji,jj,jl)*aire(ji,jj) / 1.0e12
+                  vinfor(64) = vinfor(64) + t_su(ji,jj,jl)*a_i(ji,jj,jl)*aire(ji,jj) * 1.e-12_wp
                ENDIF
             END DO
@@ -386,25 +355,25 @@
          DO ji = fs_2, fs_jpim1   ! vector opt.
             IF( tms(ji,jj) == 1 ) THEN
-               zafy = 0.0
-               zamy = 0.0
+               zafy = 0._wp
+               zamy = 0._wp
                DO jl = 1, jpl
-                  IF ((o_i(ji,jj,jl) - zshift_date).LT.0.0) THEN
-                     vinfor(18) = vinfor(18) + a_i(ji,jj,jl)*aire(ji,jj) / 1.0e12 ! FY ice area
-                     vinfor(26) = vinfor(26) + v_i(ji,jj,jl)*aire(ji,jj) / 1.0e12 ! FY ice volume
+                  IF( (o_i(ji,jj,jl) - zshift_date) < 0._wp ) THEN
+                     vinfor(18) = vinfor(18) + a_i(ji,jj,jl)*aire(ji,jj) * 1.e-12_wp ! FY ice area
+                     vinfor(26) = vinfor(26) + v_i(ji,jj,jl)*aire(ji,jj) * 1.e-12_wp ! FY ice volume
                      zafy = zafy + a_i(ji,jj,jl)
-                     vinfor(50) = vinfor(50) + sm_i(ji,jj,jl)*v_i(ji,jj,jl)*aire(ji,jj) / 1.0e12 !FY ice salinity
+                     vinfor(50) = vinfor(50) + sm_i(ji,jj,jl)*v_i(ji,jj,jl)*aire(ji,jj) * 1.e-12_wp !FY ice salinity
                   ENDIF
-                  IF ((o_i(ji,jj,jl) - zshift_date).GT.0.0) THEN
-                     vinfor(20) = vinfor(20) + a_i(ji,jj,jl)*aire(ji,jj) / 1.0e12    ! MY ice area
-                     vinfor(28) = vinfor(28) + v_i(ji,jj,jl)*aire(ji,jj) / 1.0e12
-                     vinfor(52) = vinfor(52) + sm_i(ji,jj,jl)*v_i(ji,jj,jl)*aire(ji,jj) / 1.0e12 !FY ice salinity
+                  IF( (o_i(ji,jj,jl) - zshift_date) > 0._wp ) THEN
+                     vinfor(20) = vinfor(20) + a_i(ji,jj,jl)*aire(ji,jj) * 1.e-12_wp    ! MY ice area
+                     vinfor(28) = vinfor(28) + v_i(ji,jj,jl)*aire(ji,jj) * 1.e-12_wp
+                     vinfor(52) = vinfor(52) + sm_i(ji,jj,jl)*v_i(ji,jj,jl)*aire(ji,jj) * 1.e-12_wp !FY ice salinity
                      zamy = zamy + a_i(ji,jj,jl)
                   ENDIF
                END DO ! jl
                IF ((at_i(ji,jj).GT.0.15).AND.(zafy.GT.zamy)) THEN
-                  vinfor(22) = vinfor(22) + aire(ji,jj) / 1.0e12 ! Seasonal ice extent
+                  vinfor(22) = vinfor(22) + aire(ji,jj) * 1.e-12_wp ! Seasonal ice extent
                ENDIF
                IF ((at_i(ji,jj).GT.0.15).AND.(zafy.LE.zamy)) THEN
-                  vinfor(24) = vinfor(24) + aire(ji,jj) / 1.0e12 ! Perennial ice extent
+                  vinfor(24) = vinfor(24) + aire(ji,jj) * 1.e-12_wp ! Perennial ice extent
                ENDIF
             ENDIF ! tms
@@ -427,9 +396,10 @@
       naveg = 0
       DO jv = 1, nvinfo
-         vinfom(jv)=0.0
+         vinfom(jv) = 0._wp
       END DO
       !MV      ENDIF
-
+      !
    END SUBROUTINE lim_dia
+
 
    SUBROUTINE lim_dia_init
@@ -441,27 +411,18 @@
       !!
       !! ** input   : Namelist namicedia
-      !!
-      !! history :
-      !!  8.5  ! 03-08 (C. Ethe) original code
-      !!  9.0  ! 08-03 (M. Vancoppenolle) LIM3
-      !!-------------------------------------------------------------------
-      NAMELIST/namicedia/fmtinf, nfrinf, ninfo, ntmoy
-
-      INTEGER  ::   jv   ,     &  ! dummy loop indice
-         &          ntot ,     &
-         &          ndeb ,     &
-         &          irecl
-
-      REAL(wp) ::   zxx0, zxx1    ! temporary scalars
-
+      !!-------------------------------------------------------------------
+      INTEGER  ::   jv    ! dummy loop indice
+      INTEGER  ::   ntot , ndeb , irecl   ! local integers
+      REAL(wp) ::   zxx0, zxx1    ! local scalars
       CHARACTER(len=jpchinf) ::   titinf
       CHARACTER(len=50)      ::   clname
-      !!-------------------------------------------------------------------
-
-
-      ! Read Namelist namicedia
-      REWIND ( numnam_ice )
-      READ   ( numnam_ice  , namicedia )
-      IF(lwp) THEN
+      !!
+      NAMELIST/namicedia/fmtinf, nfrinf, ninfo, ntmoy
+      !!-------------------------------------------------------------------
+      !
+      REWIND( numnam_ice )             ! read namicedia namelist
+      READ  ( numnam_ice, namicedia )
+      !
+      IF(lwp) THEN                     ! control print
          WRITE(numout,*)
          WRITE(numout,*) 'lim_dia_init : ice parameters for ice diagnostics '
@@ -473,6 +434,5 @@
       ENDIF
 
-      ! masked grid cell area
-      aire(:,:) = area(:,:) * tms(:,:)
+      aire(:,:) = area(:,:) * tms(:,:) * tmask_i(:,:)      ! masked grid cell area (interior domain only)
 
       ! Titles of ice key variables :
@@ -580,12 +540,10 @@
 
       ! opening  "ice_evolu" file
-      IF ( lk_mpp ) THEN
-         WRITE(clname,FMT="('ice.evolu_',I4.4)") narea-1
-      ELSE
-         clname = 'ice.evolu'
+      IF( lk_mpp ) THEN   ;   WRITE(clname,FMT="('ice.evolu_',I4.4)") narea-1
+      ELSE                ;   clname = 'ice.evolu'
       END IF
       irecl = ( jpchinf + 1 ) * nvinfo 
       CALL ctl_opn( numevo_ice, clname, 'UNKNOWN', 'FORMATTED', 'SEQUENTIAL',    &
-         &         irecl, numout, .TRUE., narea )
+         &         irecl, numout, lwp, narea )
 
       !- ecriture de 2 lignes d''entete :
@@ -596,12 +554,11 @@
 
       !- ecriture de 2 lignes de titre :
-      WRITE(numevo_ice,'(A,I8,A,I8,A,I5)')                                      &
+      WRITE(numevo_ice,'(A,I8,A,I8,A,I5)')                 &
          'Evolution chronologique - Experience '//cexper   &
          //'   de', ndeb, ' a', nferme, ' pas', ninfo
       WRITE(numevo_ice,fmtitr) ( titvar(jv), jv = 1, nvinfo )
 
-
       !--preparation de "titvar" pour l''ecriture parmi les valeurs numeriques :
-      DO  jv = 2 , nvinfo
+      DO jv = 2 , nvinfo
          titinf     = titvar(jv)(:jpchinf)
          titvar(jv) = '  '//titinf
@@ -609,6 +566,6 @@
 
       !--Initialisation of the arrays for the accumulation
-      DO  jv = 1, nvinfo
-         vinfom(jv) = 0.
+      DO jv = 1, nvinfo
+         vinfom(jv) = 0._wp
       END DO
       naveg = 0
@@ -616,10 +573,10 @@
 1000  FORMAT( 3(A20),4(1x,I6) )
 1111  FORMAT( 3(F7.1,1X,F7.3,1X),I3,A )  
-
+      !
    END SUBROUTINE lim_dia_init
 
 #else
    !!----------------------------------------------------------------------
-   !!   Default option :                               NO LIM sea-ice model
+   !!   Default option :                             NO LIM-3 sea-ice model
    !!----------------------------------------------------------------------
 CONTAINS
