Index: branches/UKMO/2015_CO6_CO5_shelfdiagnostic/NEMOGCM/NEMO/LIM_SRC_3/limthd.F90
===================================================================
--- branches/UKMO/2015_CO6_CO5_shelfdiagnostic/NEMOGCM/NEMO/LIM_SRC_3/limthd.F90	(revision 5422)
+++ branches/UKMO/2015_CO6_CO5_shelfdiagnostic/NEMOGCM/NEMO/LIM_SRC_3/limthd.F90	(revision 5666)
@@ -22,5 +22,4 @@
    USE phycst         ! physical constants
    USE dom_oce        ! ocean space and time domain variables
-   USE oce     , ONLY : fraqsr_1lev 
    USE ice            ! LIM: sea-ice variables
    USE sbc_oce        ! Surface boundary condition: ocean fields
@@ -28,5 +27,4 @@
    USE thd_ice        ! LIM thermodynamic sea-ice variables
    USE dom_ice        ! LIM sea-ice domain
-   USE domvvl         ! domain: variable volume level
    USE limthd_dif     ! LIM: thermodynamics, vertical diffusion
    USE limthd_dh      ! LIM: thermodynamics, ice and snow thickness variation
@@ -50,6 +48,6 @@
    PRIVATE
 
-   PUBLIC   lim_thd        ! called by limstp module
-   PUBLIC   lim_thd_init   ! called by sbc_lim_init
+   PUBLIC   lim_thd         ! called by limstp module
+   PUBLIC   lim_thd_init    ! called by sbc_lim_init
 
    !! * Substitutions
@@ -92,7 +90,5 @@
       REAL(wp), PARAMETER :: zch        = 0.0057_wp       ! heat transfer coefficient
       !
-      REAL(wp), POINTER, DIMENSION(:,:) ::  zqsr, zqns
       !!-------------------------------------------------------------------
-      CALL wrk_alloc( jpi,jpj, zqsr, zqns )
 
       IF( nn_timing == 1 )  CALL timing_start('limthd')
@@ -136,25 +132,4 @@
       ! 2) Partial computation of forcing for the thermodynamic sea ice model.      !
       !-----------------------------------------------------------------------------!
-
-      !--- Ocean solar and non solar fluxes to be used in zqld
-      IF ( .NOT. lk_cpl ) THEN   ! --- forced case, fluxes to the lead are the same as over the ocean
-         !
-         zqsr(:,:) = qsr(:,:)      ; zqns(:,:) = qns(:,:)
-         !
-      ELSE                       ! --- coupled case, fluxes to the lead are total - intercepted
-         !
-         zqsr(:,:) = qsr_tot(:,:)  ; zqns(:,:) = qns_tot(:,:)
-         !
-         DO jl = 1, jpl
-            DO jj = 1, jpj
-               DO ji = 1, jpi
-                  zqsr(ji,jj) = zqsr(ji,jj) - qsr_ice(ji,jj,jl) * a_i_b(ji,jj,jl)
-                  zqns(ji,jj) = zqns(ji,jj) - qns_ice(ji,jj,jl) * a_i_b(ji,jj,jl)
-               END DO
-            END DO
-         END DO
-         !
-      ENDIF
-
       DO jj = 1, jpj
          DO ji = 1, jpi
@@ -167,21 +142,7 @@
             !           !  temperature and turbulent mixing (McPhee, 1992)
             !
-
             ! --- Energy received in the lead, zqld is defined everywhere (J.m-2) --- !
-            ! REMARK valid at least in forced mode from clem
-            ! precip is included in qns but not in qns_ice
-            IF ( lk_cpl ) THEN
-               zqld =  tmask(ji,jj,1) * rdt_ice *  &
-                  &    (   zqsr(ji,jj) * fraqsr_1lev(ji,jj) + zqns(ji,jj)               &   ! pfrld already included in coupled mode
-                  &    + ( pfrld(ji,jj)**rn_betas - pfrld(ji,jj) ) * sprecip(ji,jj)  *     &   ! heat content of precip
-                  &      ( cpic * ( MIN( tatm_ice(ji,jj), rt0_snow ) - rt0 ) - lfus )   &
-                  &    + ( 1._wp - pfrld(ji,jj) ) * ( tprecip(ji,jj) - sprecip(ji,jj) ) * rcp * ( tatm_ice(ji,jj) - rt0 ) )
-            ELSE
-               zqld =  tmask(ji,jj,1) * rdt_ice *  &
-                  &      ( pfrld(ji,jj) * ( zqsr(ji,jj) * fraqsr_1lev(ji,jj) + zqns(ji,jj) )    &
-                  &    + ( pfrld(ji,jj)**rn_betas - pfrld(ji,jj) ) * sprecip(ji,jj)  *             &  ! heat content of precip
-                  &      ( cpic * ( MIN( tatm_ice(ji,jj), rt0_snow ) - rt0 ) - lfus )           &
-                  &    + ( 1._wp - pfrld(ji,jj) ) * ( tprecip(ji,jj) - sprecip(ji,jj) ) * rcp * ( tatm_ice(ji,jj) - rt0 ) )
-            ENDIF
+            zqld =  tmask(ji,jj,1) * rdt_ice *  &
+               &    ( pfrld(ji,jj) * qsr_oce(ji,jj) * frq_m(ji,jj) + pfrld(ji,jj) * qns_oce(ji,jj) + qemp_oce(ji,jj) )
 
             ! --- Energy needed to bring ocean surface layer until its freezing (<0, J.m-2) --- !
@@ -210,34 +171,17 @@
             ! Net heat flux on top of ice-ocean [W.m-2]
             ! -----------------------------------------
-            !     heat flux at the ocean surface + precip
-            !   + heat flux at the ice   surface 
-            hfx_in(ji,jj) = hfx_in(ji,jj)                                                                                         & 
-               ! heat flux above the ocean
-               &    +             pfrld(ji,jj)   * ( zqns(ji,jj) + zqsr(ji,jj) )                                                  &
-               ! latent heat of precip (note that precip is included in qns but not in qns_ice)
-               &    +   ( 1._wp - pfrld(ji,jj) ) * sprecip(ji,jj) * ( cpic * ( MIN( tatm_ice(ji,jj), rt0_snow ) - rt0 ) - lfus )  &
-               &    +   ( 1._wp - pfrld(ji,jj) ) * ( tprecip(ji,jj) - sprecip(ji,jj) ) * rcp * ( tatm_ice(ji,jj) - rt0 )          &
-               ! heat flux above the ice
-               &    +   SUM(    a_i_b(ji,jj,:)   * ( qns_ice(ji,jj,:) + qsr_ice(ji,jj,:) ) )
+            hfx_in(ji,jj) = qns_tot(ji,jj) + qsr_tot(ji,jj) 
 
             ! -----------------------------------------------------------------------------
-            ! Net heat flux that is retroceded to the ocean or taken from the ocean [W.m-2]
+            ! Net heat flux on top of the ocean after ice thermo (1st step) [W.m-2]
             ! -----------------------------------------------------------------------------
             !     First  step here              :  non solar + precip - qlead - qturb
             !     Second step in limthd_dh      :  heat remaining if total melt (zq_rema) 
             !     Third  step in limsbc         :  heat from ice-ocean mass exchange (zf_mass) + solar
-            hfx_out(ji,jj) = hfx_out(ji,jj)                                                                                       & 
-               ! Non solar heat flux received by the ocean
-               &    +        pfrld(ji,jj) * zqns(ji,jj)                                                                            &
-               ! latent heat of precip (note that precip is included in qns but not in qns_ice)
-               &    +      ( pfrld(ji,jj)**rn_betas - pfrld(ji,jj) ) * sprecip(ji,jj)       &
-               &         * ( cpic * ( MIN( tatm_ice(ji,jj), rt0_snow ) - rt0 ) - lfus )  &
-               &    +      ( 1._wp - pfrld(ji,jj) ) * ( tprecip(ji,jj) - sprecip(ji,jj) ) * rcp * ( tatm_ice(ji,jj) - rt0 )       &
-               ! heat flux taken from the ocean where there is open water ice formation
-               &    -      qlead(ji,jj) * r1_rdtice                                                                               &
-               ! heat flux taken from the ocean during bottom growth/melt (fhld should be 0 while bott growth)
-               &    -      at_i(ji,jj) * fhtur(ji,jj)                                                                             &
-               &    -      at_i(ji,jj) *  fhld(ji,jj)
-
+            hfx_out(ji,jj) =   pfrld(ji,jj) * qns_oce(ji,jj) + qemp_oce(ji,jj)  &  ! Non solar heat flux received by the ocean               
+               &             - qlead(ji,jj) * r1_rdtice                         &  ! heat flux taken from the ocean where there is open water ice formation
+               &             - at_i(ji,jj) * fhtur(ji,jj)                       &  ! heat flux taken by turbulence
+               &             - at_i(ji,jj) *  fhld(ji,jj)                          ! heat flux taken during bottom growth/melt 
+                                                                                   !    (fhld should be 0 while bott growth)
          END DO
       END DO
@@ -412,6 +356,4 @@
       IF( ln_limdiahsb ) CALL lim_cons_hsm(1, 'limthd', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b)
 
-      CALL wrk_dealloc( jpi,jpj, zqsr, zqns )
-
       !------------------------------------------------------------------------------|
       !  6) Transport of ice between thickness categories.                           |
@@ -472,4 +414,5 @@
    END SUBROUTINE lim_thd 
 
+ 
    SUBROUTINE lim_thd_temp( kideb, kiut )
       !!-----------------------------------------------------------------------
@@ -570,5 +513,5 @@
          END DO
          
-         CALL tab_2d_1d( nbpb, tatm_ice_1d(1:nbpb), tatm_ice(:,:)   , jpi, jpj, npb(1:nbpb) )
+         CALL tab_2d_1d( nbpb, qprec_ice_1d(1:nbpb), qprec_ice(:,:) , jpi, jpj, npb(1:nbpb) )
          CALL tab_2d_1d( nbpb, qsr_ice_1d (1:nbpb), qsr_ice(:,:,jl) , jpi, jpj, npb(1:nbpb) )
          CALL tab_2d_1d( nbpb, fr1_i0_1d  (1:nbpb), fr1_i0          , jpi, jpj, npb(1:nbpb) )
@@ -576,8 +519,5 @@
          CALL tab_2d_1d( nbpb, qns_ice_1d (1:nbpb), qns_ice(:,:,jl) , jpi, jpj, npb(1:nbpb) )
          CALL tab_2d_1d( nbpb, ftr_ice_1d (1:nbpb), ftr_ice(:,:,jl) , jpi, jpj, npb(1:nbpb) )
-         IF( .NOT. lk_cpl ) THEN
-            CALL tab_2d_1d( nbpb, qla_ice_1d (1:nbpb), qla_ice(:,:,jl) , jpi, jpj, npb(1:nbpb) )
-            CALL tab_2d_1d( nbpb, dqla_ice_1d(1:nbpb), dqla_ice(:,:,jl), jpi, jpj, npb(1:nbpb) )
-         ENDIF
+         CALL tab_2d_1d( nbpb, evap_ice_1d (1:nbpb), evap_ice(:,:,jl), jpi, jpj, npb(1:nbpb) )
          CALL tab_2d_1d( nbpb, dqns_ice_1d(1:nbpb), dqns_ice(:,:,jl), jpi, jpj, npb(1:nbpb) )
          CALL tab_2d_1d( nbpb, t_bo_1d     (1:nbpb), t_bo            , jpi, jpj, npb(1:nbpb) )
@@ -670,5 +610,5 @@
          CALL tab_1d_2d( nbpb, qns_ice(:,:,jl), npb, qns_ice_1d(1:nbpb) , jpi, jpj)
          CALL tab_1d_2d( nbpb, ftr_ice(:,:,jl), npb, ftr_ice_1d(1:nbpb) , jpi, jpj )
-                  
+         !         
       END SELECT
 
