Index: trunk/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_lim.F90
===================================================================
--- trunk/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_lim.F90	(revision 4528)
+++ trunk/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_lim.F90	(revision 4689)
@@ -59,4 +59,5 @@
    USE prtctl          ! Print control
    USE lib_fortran     ! 
+   USE cpl_oasis3, ONLY : lk_cpl
 
 #if defined key_bdy 
@@ -68,4 +69,5 @@
 
    PUBLIC sbc_ice_lim  ! routine called by sbcmod.F90
+   PUBLIC lim_prt_state
    
    !! * Substitutions
@@ -133,5 +135,5 @@
       INTEGER, INTENT(in) ::   kblk    ! type of bulk (=3 CLIO, =4 CORE)
       !!
-      INTEGER  ::   jl      ! dummy loop index
+      INTEGER  ::   ji, jj, jl, jk      ! dummy loop index
       REAL(wp) ::   zcoef   ! local scalar
       REAL(wp), POINTER, DIMENSION(:,:,:)   ::   zalb_ice_os, zalb_ice_cs  ! albedo of the ice under overcast/clear sky
@@ -146,4 +148,6 @@
       REAL(wp), POINTER, DIMENSION(:,:) :: z_dqns_ice_all  ! Mean d(qns)/dT over all categories
       REAL(wp), POINTER, DIMENSION(:,:) :: z_dqla_ice_all  ! Mean d(qla)/dT over all categories
+      REAL(wp) ::   ztmelts           ! clem 2014: for HC diags
+      REAL(wp) ::   epsi20 = 1.e-20   !
       !!----------------------------------------------------------------------
 
@@ -152,11 +156,10 @@
       IF( nn_timing == 1 )  CALL timing_start('sbc_ice_lim')
 
-      CALL wrk_alloc( jpi,jpj,jpl, zalb_ice_os, zalb_ice_cs )
-
-#if defined key_coupled
-      IF ( ln_cpl .OR. ln_iceflx_ave .OR. ln_iceflx_linear ) CALL wrk_alloc( jpi,jpj,jpl, zalb_ice)
-      IF ( ln_iceflx_ave .OR. ln_iceflx_linear ) &
-         &   CALL wrk_alloc( jpi,jpj, ztem_ice_all, zalb_ice_all, z_qsr_ice_all, z_qns_ice_all, z_qla_ice_all, z_dqns_ice_all, z_dqla_ice_all)
-#endif
+      CALL wrk_alloc( jpi,jpj,jpl, zalb_ice_os, zalb_ice_cs, zalb_ice )
+
+      IF( lk_cpl ) THEN
+         IF ( ln_iceflx_ave .OR. ln_iceflx_linear ) &
+            &   CALL wrk_alloc( jpi,jpj, ztem_ice_all, zalb_ice_all, z_qsr_ice_all, z_qns_ice_all, z_qla_ice_all, z_dqns_ice_all, z_dqla_ice_all)
+      ENDIF
 
       IF( kt == nit000 ) THEN
@@ -168,5 +171,5 @@
          !
          IF( ln_nicep ) THEN      ! control print at a given point
-            jiindx = 177   ;   jjindx = 112
+            jiindx = 15    ;   jjindx =  44
             IF(lwp) WRITE(numout,*) ' The debugging point is : jiindx : ',jiindx, ' jjindx : ',jjindx
          ENDIF
@@ -176,12 +179,13 @@
       IF( MOD( kt-1, nn_fsbc ) == 0 ) THEN     !  Ice time-step only  !
          !                                     !----------------------!
-         !                                           !  Bulk Formulea !
+         !                                           !  Bulk Formulae !
          !                                           !----------------!
          !
          u_oce(:,:) = ssu_m(:,:)                     ! mean surface ocean current at ice velocity point
          v_oce(:,:) = ssv_m(:,:)                     ! (C-grid dynamics :  U- & V-points as the ocean)
-         !
-         t_bo(:,:) = tfreez( sss_m ) +  rt0          ! masked sea surface freezing temperature [Kelvin]
-         !                                           ! (set to rt0 over land)
+
+         ! masked sea surface freezing temperature [Kelvin]
+         t_bo(:,:) = ( tfreez( sss_m ) +  rt0 ) * tmask(:,:,1) + rt0 * ( 1. - tmask(:,:,1) )
+
          CALL albedo_ice( t_su, ht_i, ht_s, zalb_ice_cs, zalb_ice_os )  ! ... ice albedo
 
@@ -192,13 +196,13 @@
          IF ( ln_cpl ) zalb_ice (:,:,:) = 0.5 * ( zalb_ice_cs (:,:,:) +  zalb_ice_os (:,:,:) )
          
-#if defined key_coupled
-         IF ( ln_iceflx_ave .OR. ln_iceflx_linear ) THEN
-            !
-            ! Compute mean albedo and temperature
-            zalb_ice_all (:,:) = fice_ice_ave ( zalb_ice (:,:,:) ) 
-            ztem_ice_all (:,:) = fice_ice_ave ( tn_ice   (:,:,:) ) 
-            !
+         IF( lk_cpl ) THEN
+            IF ( ln_iceflx_ave .OR. ln_iceflx_linear ) THEN
+               !
+               ! Compute mean albedo and temperature
+               zalb_ice_all (:,:) = fice_ice_ave ( zalb_ice (:,:,:) ) 
+               ztem_ice_all (:,:) = fice_ice_ave ( tn_ice   (:,:,:) ) 
+               !
+            ENDIF
          ENDIF
-#endif
                                                ! Bulk formulea - provides the following fields:
          ! utau_ice, vtau_ice : surface ice stress                     (U- & V-points)   [N/m2]
@@ -218,5 +222,13 @@
             !         
          CASE( 4 )                                       ! CORE bulk formulation
-            CALL blk_ice_core( t_su , u_ice     , v_ice     , zalb_ice_cs,               &
+            ! MV 2014
+            ! We must account for cloud fraction in the computation of the albedo
+            ! The present ref just uses the clear sky value
+            ! The overcast sky value is 0.06 higher, and polar skies are mostly overcast
+            ! CORE has no cloud fraction, hence we must prescribe it
+            ! Mean summer cloud fraction computed from CLIO = 0.81
+            zalb_ice(:,:,:) = 0.19 * zalb_ice_cs(:,:,:) + 0.81 * zalb_ice_os(:,:,:)
+            ! Following line, we replace zalb_ice_cs by simply zalb_ice
+            CALL blk_ice_core( t_su , u_ice     , v_ice     , zalb_ice   ,               &
                &                      utau_ice  , vtau_ice  , qns_ice    , qsr_ice   ,   &
                &                      qla_ice   , dqns_ice  , dqla_ice   ,               &
@@ -239,5 +251,5 @@
 
          ! Average over all categories
-#if defined key_coupled
+         IF( lk_cpl ) THEN
          IF ( ln_iceflx_ave .OR. ln_iceflx_linear ) THEN
 
@@ -269,5 +281,5 @@
             END IF
          END IF
-#endif
+         ENDIF
          !                                           !----------------------!
          !                                           ! LIM-3  time-stepping !
@@ -285,8 +297,8 @@
          old_smv_i(:,:,:)   = smv_i(:,:,:)     ! salt content
          old_oa_i (:,:,:)   = oa_i (:,:,:)     ! areal age content
-         !
-         old_u_ice(:,:) = u_ice(:,:)
-         old_v_ice(:,:) = v_ice(:,:)
-         !                                           ! intialisation to zero    !!gm is it truly necessary ???
+         old_u_ice(:,:)     = u_ice(:,:)
+         old_v_ice(:,:)     = v_ice(:,:)
+
+         ! trends    !!gm is it truly necessary ???
          d_a_i_thd  (:,:,:)   = 0._wp   ;   d_a_i_trp  (:,:,:)   = 0._wp
          d_v_i_thd  (:,:,:)   = 0._wp   ;   d_v_i_trp  (:,:,:)   = 0._wp
@@ -296,33 +308,40 @@
          d_smv_i_thd(:,:,:)   = 0._wp   ;   d_smv_i_trp(:,:,:)   = 0._wp
          d_oa_i_thd (:,:,:)   = 0._wp   ;   d_oa_i_trp (:,:,:)   = 0._wp
-         !
-         d_u_ice_dyn(:,:) = 0._wp
-         d_v_ice_dyn(:,:) = 0._wp
-         !
-         sfx    (:,:) = 0._wp   ;   sfx_thd  (:,:) = 0._wp
-         sfx_bri(:,:) = 0._wp   ;   sfx_mec  (:,:) = 0._wp   ;   sfx_res  (:,:) = 0._wp
-         fhbri  (:,:) = 0._wp   ;   fheat_mec(:,:) = 0._wp   ;   fheat_res(:,:) = 0._wp
-         fhmec  (:,:) = 0._wp   ;   
-         fmmec  (:,:) = 0._wp
-         fmmflx (:,:) = 0._wp     
-         focea2D(:,:) = 0._wp
-         fsup2D (:,:) = 0._wp
-
-         ! used in limthd.F90
-         rdvosif(:,:) = 0._wp   ! variation of ice volume at surface
-         rdvobif(:,:) = 0._wp   ! variation of ice volume at bottom
-         fdvolif(:,:) = 0._wp   ! total variation of ice volume
-         rdvonif(:,:) = 0._wp   ! lateral variation of ice volume
-         fstric (:,:) = 0._wp   ! part of solar radiation transmitted through the ice
-         ffltbif(:,:) = 0._wp   ! linked with fstric
-         qfvbq  (:,:) = 0._wp   ! linked with fstric
-         rdm_snw(:,:) = 0._wp   ! variation of snow mass per unit area
-         rdm_ice(:,:) = 0._wp   ! variation of ice mass per unit area
-         hicifp (:,:) = 0._wp   ! daily thermodynamic ice production. 
-         !
-         diag_sni_gr(:,:) = 0._wp   ;   diag_lat_gr(:,:) = 0._wp
-         diag_bot_gr(:,:) = 0._wp   ;   diag_dyn_gr(:,:) = 0._wp
-         diag_bot_me(:,:) = 0._wp   ;   diag_sur_me(:,:) = 0._wp
-         diag_res_pr(:,:) = 0._wp   ;   diag_trp_vi(:,:) = 0._wp
+         d_u_ice_dyn(:,:)     = 0._wp   ;   d_v_ice_dyn(:,:)     = 0._wp
+
+         ! salt, heat and mass fluxes
+         sfx    (:,:) = 0._wp   ;
+         sfx_bri(:,:) = 0._wp   ;   sfx_dyn(:,:) = 0._wp 
+         sfx_sni(:,:) = 0._wp   ;   sfx_opw(:,:) = 0._wp
+         sfx_bog(:,:) = 0._wp   ;   sfx_dyn(:,:) = 0._wp
+         sfx_bom(:,:) = 0._wp   ;   sfx_sum(:,:) = 0._wp
+         sfx_res(:,:) = 0._wp
+
+         wfx_snw(:,:) = 0._wp   ;   wfx_ice(:,:) = 0._wp
+         wfx_sni(:,:) = 0._wp   ;   wfx_opw(:,:) = 0._wp
+         wfx_bog(:,:) = 0._wp   ;   wfx_dyn(:,:) = 0._wp
+         wfx_bom(:,:) = 0._wp   ;   wfx_sum(:,:) = 0._wp
+         wfx_res(:,:) = 0._wp   ;   wfx_sub(:,:) = 0._wp
+         wfx_spr(:,:) = 0._wp   ;   
+
+         hfx_in (:,:) = 0._wp   ;   hfx_out(:,:) = 0._wp
+         hfx_thd(:,:) = 0._wp   ;   
+         hfx_snw(:,:) = 0._wp   ;   hfx_opw(:,:) = 0._wp
+         hfx_bog(:,:) = 0._wp   ;   hfx_dyn(:,:) = 0._wp
+         hfx_bom(:,:) = 0._wp   ;   hfx_sum(:,:) = 0._wp
+         hfx_res(:,:) = 0._wp   ;   hfx_sub(:,:) = 0._wp
+         hfx_spr(:,:) = 0._wp   ;   hfx_dif(:,:) = 0._wp 
+         hfx_err(:,:) = 0._wp   ;   hfx_err_rem(:,:) = 0._wp
+
+         !
+         fhld  (:,:)    = 0._wp 
+         fmmflx(:,:)    = 0._wp     
+         ! part of solar radiation transmitted through the ice
+         ftr_ice(:,:,:) = 0._wp
+
+         ! diags
+         diag_trp_vi  (:,:) = 0._wp  ; diag_trp_vs(:,:) = 0._wp  ;  diag_trp_ei(:,:) = 0._wp  ;  diag_trp_es(:,:) = 0._wp
+         diag_heat_dhc(:,:) = 0._wp  
+
          ! dynamical invariants
          delta_i(:,:) = 0._wp       ;   divu_i(:,:) = 0._wp       ;   shear_i(:,:) = 0._wp
@@ -375,5 +394,4 @@
                           zcoef = rdt_ice /rday           !  Ice natural aging
                           oa_i(:,:,:) = oa_i(:,:,:) + a_i(:,:,:) * zcoef
-                          CALL lim_var_glo2eqv            ! this CALL is maybe not necessary (Martin)
          IF( ln_nicep )   CALL lim_prt_state( kt, jiindx, jjindx, 1, ' - ice thermodyn. - ' )   ! control print
                           CALL lim_itd_th( kt )           !  Remap ice categories, lateral accretion  !
@@ -391,7 +409,7 @@
          !                                           ! Diagnostics and outputs 
          IF (ln_limdiaout) CALL lim_diahsb
-!clem # if ! defined key_iomput
+
                           CALL lim_wri( 1  )              ! Ice outputs 
-!clem # endif
+
          IF( kt == nit000 .AND. ln_rstart )   &
             &             CALL iom_close( numrir )        ! clem: close input ice restart file
@@ -413,12 +431,10 @@
       
 !!gm   remark, the ocean-ice stress is not saved in ice diag call above .....  find a solution!!!
-      !
-      CALL wrk_dealloc( jpi,jpj,jpl, zalb_ice_os, zalb_ice_cs )
-
-#if defined key_coupled
-      IF ( ln_cpl .OR. ln_iceflx_ave .OR. ln_iceflx_linear ) CALL wrk_dealloc( jpi,jpj,jpl, zalb_ice)
-      IF ( ln_iceflx_ave .OR. ln_iceflx_linear ) &
-         &    CALL wrk_dealloc( jpi,jpj, ztem_ice_all, zalb_ice_all, z_qsr_ice_all, z_qns_ice_all, z_qla_ice_all, z_dqns_ice_all, z_dqla_ice_all)
-#endif
+      CALL wrk_dealloc( jpi,jpj,jpl, zalb_ice_os, zalb_ice_cs, zalb_ice )
+
+      IF( lk_cpl ) THEN
+         IF ( ln_iceflx_ave .OR. ln_iceflx_linear ) &
+            &    CALL wrk_dealloc( jpi,jpj, ztem_ice_all, zalb_ice_all, z_qsr_ice_all, z_qns_ice_all, z_qla_ice_all, z_dqns_ice_all, z_dqla_ice_all)
+      ENDIF
       !
       IF( nn_timing == 1 )  CALL timing_stop('sbc_ice_lim')
@@ -534,5 +550,4 @@
 !                 WRITE(numout,*) ' sst                  : ', sst_m(ji,jj)
 !                 WRITE(numout,*) ' sss                  : ', sss_m(ji,jj)
-!                 WRITE(numout,*) ' s_i_newice           : ', s_i_newice(ji,jj,1:jpl)
 !                 WRITE(numout,*) 
                   inb_alp(ialert_id) = inb_alp(ialert_id) + 1
@@ -591,17 +606,4 @@
                !WRITE(numout,*) ' sst       : ', sst_m(ji,jj)
                !WRITE(numout,*) ' sss       : ', sss_m(ji,jj)
-               !WRITE(numout,*) ' qcmif     : ', qcmif(ji,jj)
-               !WRITE(numout,*) ' qldif     : ', qldif(ji,jj)
-               !WRITE(numout,*) ' qcmif     : ', qcmif(ji,jj) / rdt_ice
-               !WRITE(numout,*) ' qldif     : ', qldif(ji,jj) / rdt_ice
-               !WRITE(numout,*) ' qfvbq     : ', qfvbq(ji,jj)
-               !WRITE(numout,*) ' qdtcn     : ', qdtcn(ji,jj)
-               !WRITE(numout,*) ' qfvbq / dt: ', qfvbq(ji,jj) / rdt_ice
-               !WRITE(numout,*) ' qdtcn / dt: ', qdtcn(ji,jj) / rdt_ice
-               !WRITE(numout,*) ' fdtcn     : ', fdtcn(ji,jj) 
-               !WRITE(numout,*) ' fhmec     : ', fhmec(ji,jj) 
-               !WRITE(numout,*) ' fheat_mec : ', fheat_mec(ji,jj) 
-               !WRITE(numout,*) ' fheat_res : ', fheat_res(ji,jj) 
-               !WRITE(numout,*) ' fhbri     : ', fhbri(ji,jj) 
                !
                !CALL lim_prt_state( kt, ji, jj, 2, '   ')
@@ -790,14 +792,8 @@
                WRITE(numout,*) ' - Heat / FW fluxes '
                WRITE(numout,*) '   ~~~~~~~~~~~~~~~~ '
-               WRITE(numout,*) ' emp        : ', emp      (ji,jj)
-               WRITE(numout,*) ' sfx        : ', sfx      (ji,jj)
-               WRITE(numout,*) ' sfx_thd    : ', sfx_thd(ji,jj)
-               WRITE(numout,*) ' sfx_bri    : ', sfx_bri  (ji,jj)
-               WRITE(numout,*) ' sfx_mec    : ', sfx_mec  (ji,jj)
-               WRITE(numout,*) ' sfx_res    : ', sfx_res(ji,jj)
-               WRITE(numout,*) ' fmmec      : ', fmmec    (ji,jj)
-               WRITE(numout,*) ' fhmec      : ', fhmec    (ji,jj)
-               WRITE(numout,*) ' fhbri      : ', fhbri    (ji,jj)
-               WRITE(numout,*) ' fheat_mec  : ', fheat_mec(ji,jj)
+               WRITE(numout,*) ' - Heat fluxes in and out the ice ***'
+               WRITE(numout,*) ' qsr_ini       : ', pfrld(ji,jj) * qsr(ji,jj) + SUM( old_a_i(ji,jj,:) * qsr_ice(ji,jj,:) )
+               WRITE(numout,*) ' qns_ini       : ', pfrld(ji,jj) * qns(ji,jj) + SUM( old_a_i(ji,jj,:) * qns_ice(ji,jj,:) )
+               WRITE(numout,*)
                WRITE(numout,*) 
                WRITE(numout,*) ' sst        : ', sst_m(ji,jj)  
@@ -829,21 +825,27 @@
                WRITE(numout,*) ' qsr       : ', qsr(ji,jj)
                WRITE(numout,*) ' qns       : ', qns(ji,jj)
-               WRITE(numout,*) ' fdtcn     : ', fdtcn(ji,jj)
-               WRITE(numout,*) ' qcmif     : ', qcmif(ji,jj) * r1_rdtice
-               WRITE(numout,*) ' qldif     : ', qldif(ji,jj) * r1_rdtice
+               WRITE(numout,*)
+               WRITE(numout,*) ' hfx_mass     : ', hfx_thd(ji,jj) + hfx_dyn(ji,jj) + hfx_snw(ji,jj) + hfx_res(ji,jj)
+               WRITE(numout,*) ' hfx_in       : ', hfx_in(ji,jj)
+               WRITE(numout,*) ' hfx_out      : ', hfx_out(ji,jj)
+               WRITE(numout,*) ' dhc          : ', diag_heat_dhc(ji,jj)              
+               WRITE(numout,*)
+               WRITE(numout,*) ' hfx_dyn      : ', hfx_dyn(ji,jj)
+               WRITE(numout,*) ' hfx_thd      : ', hfx_thd(ji,jj)
+               WRITE(numout,*) ' hfx_res      : ', hfx_res(ji,jj)
+               WRITE(numout,*) ' fhtur        : ', fhtur(ji,jj) 
+               WRITE(numout,*) ' qlead        : ', qlead(ji,jj) * r1_rdtice
                WRITE(numout,*)
                WRITE(numout,*) ' - Salt fluxes at bottom interface ***'
                WRITE(numout,*) ' emp       : ', emp    (ji,jj)
-               WRITE(numout,*) ' sfx_bri   : ', sfx_bri(ji,jj)
                WRITE(numout,*) ' sfx       : ', sfx    (ji,jj)
                WRITE(numout,*) ' sfx_res   : ', sfx_res(ji,jj)
-               WRITE(numout,*) ' sfx_mec   : ', sfx_mec(ji,jj)
-               WRITE(numout,*) ' - Heat fluxes at bottom interface ***'
-               WRITE(numout,*) ' fheat_res : ', fheat_res(ji,jj)
+               WRITE(numout,*) ' sfx_bri   : ', sfx_bri(ji,jj)
+               WRITE(numout,*) ' sfx_dyn   : ', sfx_dyn(ji,jj)
                WRITE(numout,*)
                WRITE(numout,*) ' - Momentum fluxes '
                WRITE(numout,*) ' utau      : ', utau(ji,jj) 
                WRITE(numout,*) ' vtau      : ', vtau(ji,jj)
-            ENDIF
+            ENDIF 
             WRITE(numout,*) ' '
             !
