Index: branches/2013/dev_LOCEAN_2013/NEMOGCM/NEMO/LIM_SRC_3/limthd.F90
===================================================================
--- branches/2013/dev_LOCEAN_2013/NEMOGCM/NEMO/LIM_SRC_3/limthd.F90	(revision 4147)
+++ branches/2013/dev_LOCEAN_2013/NEMOGCM/NEMO/LIM_SRC_3/limthd.F90	(revision 4161)
@@ -11,4 +11,5 @@
    !!            3.3  ! 2010-11 (G. Madec) corrected snow melting heat (due to factor betas)
    !!            4.0  ! 2011-02 (G. Madec) dynamical allocation
+   !!             -   ! 2012-05 (C. Rousset) add penetration solar flux
    !!----------------------------------------------------------------------
 #if defined key_lim3
@@ -40,4 +41,5 @@
    USE prtctl         ! Print control
    USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)  
+   USE timing         ! Timing
 
    IMPLICIT NONE
@@ -92,8 +94,22 @@
       REAL(wp) ::   zfntlat, zpareff, zareamin, zcoef   !    -         -
       REAL(wp), POINTER, DIMENSION(:,:) ::   zqlbsbq   ! link with lead energy budget qldif
+      REAL(wp) :: zchk_v_i, zchk_smv, zchk_fs, zchk_fw, zchk_v_i_b, zchk_smv_b, zchk_fs_b, zchk_fw_b ! Check conservation (C Rousset)
+      REAL(wp) :: zchk_vmin, zchk_amin, zchk_amax ! Check errors (C Rousset)
       !!-------------------------------------------------------------------
+      IF( nn_timing == 1 )  CALL timing_start('limthd')
 
       CALL wrk_alloc( jpi, jpj, zqlbsbq )
    
+      ! -------------------------------
+      !- check conservation (C Rousset)
+      IF (ln_limdiahsb) THEN
+         zchk_v_i_b = glob_sum( SUM(   v_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) )
+         zchk_smv_b = glob_sum( SUM( smv_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) )
+         zchk_fw_b  = glob_sum( rdm_ice(:,:) * area(:,:) * tms(:,:) )
+         zchk_fs_b  = glob_sum( ( sfx_bri(:,:) + sfx_thd(:,:) + sfx_res(:,:) + sfx_mec(:,:) ) * area(:,:) * tms(:,:) )
+      ENDIF
+      !- check conservation (C Rousset)
+      ! -------------------------------
+
       !------------------------------------------------------------------------------!
       ! 1) Initialization of diagnostic variables                                    !
@@ -109,5 +125,5 @@
                DO ji = 1, jpi
                   !Energy of melting q(S,T) [J.m-3]
-                  e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) / ( area(ji,jj) * MAX( v_i(ji,jj,jl) , epsi06 ) ) * nlay_i
+                  e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) / ( area(ji,jj) * MAX( v_i(ji,jj,jl) , epsi06 ) ) * REAL( nlay_i )
                   !0 if no ice and 1 if yes
                   zindb = 1.0 - MAX(  0.0 , SIGN( 1.0 , - ht_i(ji,jj,jl) )  ) 
@@ -121,5 +137,5 @@
                DO ji = 1, jpi
                   !Energy of melting q(S,T) [J.m-3]
-                  e_s(ji,jj,jk,jl) = e_s(ji,jj,jk,jl) / ( area(ji,jj) * MAX( v_s(ji,jj,jl) , epsi06 ) ) * nlay_s
+                  e_s(ji,jj,jk,jl) = e_s(ji,jj,jk,jl) / ( area(ji,jj) * MAX( v_s(ji,jj,jl) , epsi06 ) ) * REAL( nlay_s )
                   !0 if no ice and 1 if yes
                   zindb = 1.0 - MAX(  0.0 , SIGN( 1.0 , - ht_s(ji,jj,jl) )  ) 
@@ -134,17 +150,17 @@
       ! 1.3) Set some dummies to 0
       !-----------------------------
-      rdvosif(:,:) = 0.e0   ! variation of ice volume at surface
-      rdvobif(:,:) = 0.e0   ! variation of ice volume at bottom
-      fdvolif(:,:) = 0.e0   ! total variation of ice volume
-      rdvonif(:,:) = 0.e0   ! lateral variation of ice volume
-      fstric (:,:) = 0.e0   ! part of solar radiation transmitted through the ice
-      ffltbif(:,:) = 0.e0   ! linked with fstric
-      qfvbq  (:,:) = 0.e0   ! linked with fstric
-      rdm_snw(:,:) = 0.e0   ! variation of snow mass per unit area
-      rdm_ice(:,:) = 0.e0   ! variation of ice mass per unit area
-      hicifp (:,:) = 0.e0   ! daily thermodynamic ice production. 
-      sfx_bri(:,:) = 0.e0   ! brine flux contribution to salt flux to the ocean
-      fhbri  (:,:) = 0.e0   ! brine flux contribution to heat flux to the ocean
-      sfx_thd(:,:) = 0.e0   ! equivalent salt flux to the ocean due to ice/growth decay
+      !clem rdvosif(:,:) = 0.e0   ! variation of ice volume at surface
+      !clem rdvobif(:,:) = 0.e0   ! variation of ice volume at bottom
+      !clem fdvolif(:,:) = 0.e0   ! total variation of ice volume
+      !clem rdvonif(:,:) = 0.e0   ! lateral variation of ice volume
+      !clem fstric (:,:) = 0.e0   ! part of solar radiation transmitted through the ice
+      !clem ffltbif(:,:) = 0.e0   ! linked with fstric
+      !clem qfvbq  (:,:) = 0.e0   ! linked with fstric
+      !clem rdm_snw(:,:) = 0.e0   ! variation of snow mass per unit area
+      !clem rdm_ice(:,:) = 0.e0   ! variation of ice mass per unit area
+      !clem hicifp (:,:) = 0.e0   ! daily thermodynamic ice production. 
+      !clem sfx_bri(:,:) = 0.e0   ! brine flux contribution to salt flux to the ocean
+      !clem fhbri  (:,:) = 0.e0   ! brine flux contribution to heat flux to the ocean
+      !clem sfx_thd(:,:) = 0.e0   ! equivalent salt flux to the ocean due to ice/growth decay
 
       !-----------------------------------
@@ -165,9 +181,7 @@
 !CDIR NOVERRCHK
          DO ji = 1, jpi
-            zthsnice       = SUM( ht_s(ji,jj,1:jpl) ) + SUM( ht_i(ji,jj,1:jpl) )
-            zindb          = tms(ji,jj) * ( 1.0 - MAX( zzero , SIGN( zone , - zthsnice ) ) ) 
             phicif(ji,jj)  = vt_i(ji,jj)
             pfrld(ji,jj)   = 1.0 - at_i(ji,jj)
-            zinda          = 1.0 - MAX( zzero , SIGN( zone , - ( 1.0 - pfrld(ji,jj) ) ) )
+            zinda          = tms(ji,jj) * ( 1.0 - MAX( zzero , SIGN( zone , - at_i(ji,jj) ) ) )
             !
             !           !  solar irradiance transmission at the mixed layer bottom and used in the lead heat budget
@@ -180,16 +194,16 @@
 
             ! here the drag will depend on ice thickness and type (0.006)
-            fdtcn(ji,jj)  = zindb * rau0 * rcp * 0.006  * zfric_u * ( (sst_m(ji,jj) + rt0) - t_bo(ji,jj) ) 
+            fdtcn(ji,jj)  = zinda * rau0 * rcp * 0.006  * zfric_u * ( (sst_m(ji,jj) + rt0) - t_bo(ji,jj) ) 
             ! also category dependent
             !           !-- Energy from the turbulent oceanic heat flux heat flux coming in the lead 
-            qdtcn(ji,jj)  = zindb * fdtcn(ji,jj) * (1.0 - at_i(ji,jj)) * rdt_ice
+            qdtcn(ji,jj)  = zinda * fdtcn(ji,jj) * (1.0 - at_i(ji,jj)) * rdt_ice
             !                       
             !           !-- Lead heat budget, qldif (part 1, next one is in limthd_dh) 
             !           !   caution: exponent betas used as more snow can fallinto leads
             qldif(ji,jj) =  tms(ji,jj) * rdt_ice  * (                             &
-               &   pfrld(ji,jj)        * (  qsr(ji,jj)                            &   ! solar heat
+               &   pfrld(ji,jj)        * (  qsr(ji,jj) * oatte(ji,jj)             &   ! solar heat + clem modif
                &                            + qns(ji,jj)                          &   ! non solar heat
                &                            + fdtcn(ji,jj)                        &   ! turbulent ice-ocean heat
-               &                            + fsbbq(ji,jj) * ( 1.0 - zindb )  )   &   ! residual heat from previous step
+               &                            + fsbbq(ji,jj) * ( 1.0 - zinda )  )   &   ! residual heat from previous step
                & - pfrld(ji,jj)**betas * sprecip(ji,jj) * lfus                    )   ! latent heat of sprecip melting
             !
@@ -206,8 +220,8 @@
             !
             ! Energy needed to bring ocean surface layer until its freezing (qcmif, limflx)
-            qcmif  (ji,jj) =  rau0 * rcp * fse3t(ji,jj,1) * ( t_bo(ji,jj) - (sst_m(ji,jj) + rt0) ) * ( 1. - zinda )
+            qcmif  (ji,jj) =  rau0 * rcp * fse3t(ji,jj,1) * ( t_bo(ji,jj) - (sst_m(ji,jj) + rt0) )
             !
             ! oceanic heat flux (limthd_dh)
-            fbif   (ji,jj) = zindb * (  fsbbq(ji,jj) / MAX( at_i(ji,jj) , epsi20 ) + fdtcn(ji,jj) )
+            fbif   (ji,jj) = zinda * (  fsbbq(ji,jj) / MAX( at_i(ji,jj) , epsi20 ) + fdtcn(ji,jj) )
             !
          END DO
@@ -294,4 +308,6 @@
             CALL tab_2d_1d( nbpb, qfvbq_1d   (1:nbpb), qfvbq           , jpi, jpj, npb(1:nbpb) )
 
+            CALL tab_2d_1d( nbpb, iatte_1d   (1:nbpb), iatte           , jpi, jpj, npb(1:nbpb) ) ! clem modif
+            CALL tab_2d_1d( nbpb, oatte_1d   (1:nbpb), oatte           , jpi, jpj, npb(1:nbpb) ) ! clem modif
             !--------------------------------
             ! 4.3) Thermodynamic processes
@@ -411,6 +427,6 @@
       ! 5.4) Diagnostic thermodynamic growth rates
       !--------------------------------------------
-      d_v_i_thd(:,:,:) = v_i      (:,:,:) - old_v_i(:,:,:)    ! ice volumes 
-      dv_dt_thd(:,:,:) = d_v_i_thd(:,:,:) * r1_rdtice * rday
+!clem@useless      d_v_i_thd(:,:,:) = v_i      (:,:,:) - old_v_i(:,:,:)    ! ice volumes 
+!clem@mv-to-itd    dv_dt_thd(:,:,:) = d_v_i_thd(:,:,:) * r1_rdtice * rday
 
       IF( con_i )   fbif(:,:) = fbif(:,:) + zqlbsbq(:,:)
@@ -448,6 +464,31 @@
       ENDIF
       !
+      ! -------------------------------
+      !- check conservation (C Rousset)
+      IF (ln_limdiahsb) THEN
+         zchk_fs  = glob_sum( ( sfx_bri(:,:) + sfx_thd(:,:) + sfx_res(:,:) + sfx_mec(:,:) ) * area(:,:) * tms(:,:) ) - zchk_fs_b
+         zchk_fw  = glob_sum( rdm_ice(:,:) * area(:,:) * tms(:,:) ) - zchk_fw_b
+ 
+         zchk_v_i = ( glob_sum( SUM(   v_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) ) - zchk_v_i_b - ( zchk_fw / rhoic ) ) * r1_rdtice
+         zchk_smv = ( glob_sum( SUM( smv_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) ) - zchk_smv_b ) * r1_rdtice + ( zchk_fs / rhoic )
+
+         zchk_vmin = glob_min(v_i)
+         zchk_amax = glob_max(SUM(a_i,dim=3))
+         zchk_amin = glob_min(a_i)
+       
+         IF(lwp) THEN
+            IF ( ABS( zchk_v_i   ) >  1.e-5 ) WRITE(numout,*) 'violation volume [m3/day]     (limthd) = ',(zchk_v_i * rday)
+            IF ( ABS( zchk_smv   ) >  1.e-4 ) WRITE(numout,*) 'violation saline [psu*m3/day] (limthd) = ',(zchk_smv * rday)
+            IF ( zchk_vmin <  0.            ) WRITE(numout,*) 'violation v_i<0  [mm]         (limthd) = ',(zchk_vmin * 1.e-3)
+            IF ( zchk_amax >  amax+epsi10   ) WRITE(numout,*) 'violation a_i>amax            (limthd) = ',zchk_amax
+            IF ( zchk_amin <  0.            ) WRITE(numout,*) 'violation a_i<0               (limthd) = ',zchk_amin
+         ENDIF
+      ENDIF
+      !- check conservation (C Rousset)
+      ! -------------------------------
+      !
       CALL wrk_dealloc( jpi, jpj, zqlbsbq )
       !
+      IF( nn_timing == 1 )  CALL timing_stop('limthd')
    END SUBROUTINE lim_thd
 
@@ -472,10 +513,10 @@
       DO jk = 1, nlay_i                ! total q over all layers, ice [J.m-2]
          DO ji = kideb, kiut
-            etilayer(ji,jk) = q_i_b(ji,jk) * ht_i_b(ji) / nlay_i
+            etilayer(ji,jk) = q_i_b(ji,jk) * ht_i_b(ji) / REAL( nlay_i )
             eti     (ji,jl) = eti(ji,jl) + etilayer(ji,jk) 
          END DO
       END DO
       DO ji = kideb, kiut              ! total q over all layers, snow [J.m-2]
-         ets(ji,jl) = ets(ji,jl) + q_s_b(ji,1) * ht_s_b(ji) / nlay_s
+         ets(ji,jl) = ets(ji,jl) + q_s_b(ji,1) * ht_s_b(ji) / REAL( nlay_s )
       END DO
       !
@@ -498,5 +539,5 @@
 
       INTEGER  ::   ji, jk         ! loop indices
-      INTEGER  ::   zji, zjj
+      INTEGER  ::   ii, ij
       INTEGER  ::   numce          ! number of points for which conservation is violated
       REAL(wp) ::   meance         ! mean conservation error
@@ -521,8 +562,8 @@
       !----------------------------------------
       DO ji = kideb, kiut
-         zji = MOD( npb(ji) - 1 , jpi ) + 1
-         zjj =    ( npb(ji) - 1 ) / jpi + 1
+         ii = MOD( npb(ji) - 1 , jpi ) + 1
+         ij =    ( npb(ji) - 1 ) / jpi + 1
          fatm     (ji,jl) = qnsr_ice_1d(ji) + ( 1._wp - i0(ji) ) * qsr_ice_1d(ji)
-         sum_fluxq(ji,jl) = fc_su(ji) - fc_bo_i(ji) + qsr_ice_1d(ji) * i0(ji) - fstroc(zji,zjj,jl)
+         sum_fluxq(ji,jl) = fc_su(ji) - fc_bo_i(ji) + qsr_ice_1d(ji) * i0(ji) - fstroc(ii,ij,jl)
       END DO
 
@@ -579,6 +620,6 @@
          IF ( ( ( t_su_b(ji) .LT. rtt ) .AND. ( surf_error(ji,jl) .GT. max_surf_err ) ) .OR. &
             ( cons_error(ji,jl) .GT. max_cons_err  ) ) THEN
-            zji                 = MOD( npb(ji) - 1, jpi ) + 1
-            zjj                 = ( npb(ji) - 1 ) / jpi + 1
+            ii                 = MOD( npb(ji) - 1, jpi ) + 1
+            ij                 = ( npb(ji) - 1 ) / jpi + 1
             !
             WRITE(numout,*) ' alerte 1     '
@@ -586,6 +627,6 @@
             WRITE(numout,*) ' heat diffusion in the ice '
             WRITE(numout,*) ' Category   : ', jl
-            WRITE(numout,*) ' zji , zjj  : ', zji, zjj
-            WRITE(numout,*) ' lat, lon   : ', gphit(zji,zjj), glamt(zji,zjj)
+            WRITE(numout,*) ' ii , ij  : ', ii, ij
+            WRITE(numout,*) ' lat, lon   : ', gphit(ii,ij), glamt(ii,ij)
             WRITE(numout,*) ' cons_error : ', cons_error(ji,jl)
             WRITE(numout,*) ' surf_error : ', surf_error(ji,jl)
@@ -615,5 +656,5 @@
             WRITE(numout,*) ' fc_bo      : ', - fc_bo_i  (ji)
             WRITE(numout,*) ' foc        : ', fbif_1d(ji)
-            WRITE(numout,*) ' fstroc     : ', fstroc   (zji,zjj,jl)
+            WRITE(numout,*) ' fstroc     : ', fstroc   (ii,ij,jl)
             WRITE(numout,*) ' i0         : ', i0(ji)
             WRITE(numout,*) ' qsr_ice    : ', (1.0-i0(ji))*qsr_ice_1d(ji)
@@ -651,5 +692,5 @@
       !
       INTEGER  ::   ji                ! loop indices
-      INTEGER  ::   zji, zjj, numce         ! local integers
+      INTEGER  ::   ii, ij, numce         ! local integers
       REAL(wp) ::   meance, max_cons_err    !local scalar
       !!---------------------------------------------------------------------
@@ -669,9 +710,9 @@
       !----------------------------------------
       DO ji = kideb, kiut
-         zji = MOD( npb(ji) - 1 , jpi ) + 1
-         zjj =    ( npb(ji) - 1 ) / jpi + 1
+         ii = MOD( npb(ji) - 1 , jpi ) + 1
+         ij =    ( npb(ji) - 1 ) / jpi + 1
 
          fatm      (ji,jl) = qnsr_ice_1d(ji) + qsr_ice_1d(ji)                       ! total heat flux
-         sum_fluxq (ji,jl) = fatm(ji,jl) + fbif_1d(ji) - ftotal_fin(ji) - fstroc(zji,zjj,jl) 
+         sum_fluxq (ji,jl) = fatm(ji,jl) + fbif_1d(ji) - ftotal_fin(ji) - fstroc(ii,ij,jl) 
          cons_error(ji,jl) = ABS( dq_i(ji,jl) * r1_rdtice + sum_fluxq(ji,jl) )
       END DO
@@ -706,11 +747,11 @@
       DO ji = kideb, kiut
          IF ( cons_error(ji,jl) .GT. max_cons_err  ) THEN
-            zji = MOD( npb(ji) - 1, jpi ) + 1
-            zjj =    ( npb(ji) - 1 ) / jpi + 1
+            ii = MOD( npb(ji) - 1, jpi ) + 1
+            ij =    ( npb(ji) - 1 ) / jpi + 1
             !
             WRITE(numout,*) ' alerte 1 - category : ', jl
             WRITE(numout,*) ' Untolerated conservation error after limthd_ent '
-            WRITE(numout,*) ' zji , zjj  : ', zji, zjj
-            WRITE(numout,*) ' lat, lon   : ', gphit(zji,zjj), glamt(zji,zjj)
+            WRITE(numout,*) ' ii , ij  : ', ii, ij
+            WRITE(numout,*) ' lat, lon   : ', gphit(ii,ij), glamt(ii,ij)
             WRITE(numout,*) ' * '
             WRITE(numout,*) ' Ftotal     : ', sum_fluxq(ji,jl)
@@ -724,5 +765,5 @@
             WRITE(numout,*) ' foce       : ', fbif_1d(ji)
             WRITE(numout,*) ' fres       : ', ftotal_fin(ji)
-            WRITE(numout,*) ' fhbri      : ', fhbricat(zji,zjj,jl)
+            WRITE(numout,*) ' fhbri      : ', fhbricat(ii,ij,jl)
             WRITE(numout,*) ' * '
             WRITE(numout,*) ' Heat contents --- : '
@@ -793,5 +834,5 @@
       INTEGER  ::   ios                 ! Local integer output status for namelist read
       NAMELIST/namicethd/ hmelt , hiccrit, fraz_swi, maxfrazb, vfrazb, Cfrazb,   &
-         &                hicmin, hiclim, amax  ,                                &
+         &                hicmin, hiclim,                                        &
          &                sbeta  , parlat, hakspl, hibspl, exld,                 &
          &                hakdif, hnzst  , thth  , parsub, alphs, betas,         & 
@@ -825,5 +866,4 @@
          WRITE(numout,*)'      ice thick. corr. to max. energy stored in brine pocket  hicmin       = ', hicmin  
          WRITE(numout,*)'      minimum ice thickness                                   hiclim       = ', hiclim 
-         WRITE(numout,*)'      maximum lead fraction                                   amax         = ', amax 
          WRITE(numout,*)'      numerical carac. of the scheme for diffusion in ice '
          WRITE(numout,*)'      Cranck-Nicholson (=0.5), implicit (=1), explicit (=0)   sbeta        = ', sbeta
