Index: branches/2015/dev_r5021_UKMO1_CICE_coupling/NEMOGCM/NEMO/LIM_SRC_3/limistate.F90
===================================================================
--- branches/2015/dev_r5021_UKMO1_CICE_coupling/NEMOGCM/NEMO/LIM_SRC_3/limistate.F90	(revision 5234)
+++ branches/2015/dev_r5021_UKMO1_CICE_coupling/NEMOGCM/NEMO/LIM_SRC_3/limistate.F90	(revision 5443)
@@ -22,5 +22,4 @@
    USE eosbn2           ! equation of state
    USE ice              ! sea-ice variables
-   USE par_ice          ! ice parameters
    USE par_oce          ! ocean parameters
    USE dom_ice          ! sea-ice domain
@@ -36,17 +35,17 @@
 
    !                          !!** init namelist (namiceini) **
-   REAL(wp) ::   thres_sst   ! threshold water temperature for initial sea ice
-   REAL(wp) ::   hts_ini_n   ! initial snow thickness in the north
-   REAL(wp) ::   hts_ini_s   ! initial snow thickness in the south
-   REAL(wp) ::   hti_ini_n   ! initial ice thickness in the north
-   REAL(wp) ::   hti_ini_s   ! initial ice thickness in the south
-   REAL(wp) ::   ati_ini_n   ! initial leads area in the north
-   REAL(wp) ::   ati_ini_s   ! initial leads area in the south
-   REAL(wp) ::   smi_ini_n   ! initial salinity 
-   REAL(wp) ::   smi_ini_s   ! initial salinity
-   REAL(wp) ::   tmi_ini_n   ! initial temperature
-   REAL(wp) ::   tmi_ini_s   ! initial temperature
-
-   LOGICAL  ::  ln_limini    ! initialization or not
+   REAL(wp) ::   rn_thres_sst   ! threshold water temperature for initial sea ice
+   REAL(wp) ::   rn_hts_ini_n   ! initial snow thickness in the north
+   REAL(wp) ::   rn_hts_ini_s   ! initial snow thickness in the south
+   REAL(wp) ::   rn_hti_ini_n   ! initial ice thickness in the north
+   REAL(wp) ::   rn_hti_ini_s   ! initial ice thickness in the south
+   REAL(wp) ::   rn_ati_ini_n   ! initial leads area in the north
+   REAL(wp) ::   rn_ati_ini_s   ! initial leads area in the south
+   REAL(wp) ::   rn_smi_ini_n   ! initial salinity 
+   REAL(wp) ::   rn_smi_ini_s   ! initial salinity
+   REAL(wp) ::   rn_tmi_ini_n   ! initial temperature
+   REAL(wp) ::   rn_tmi_ini_s   ! initial temperature
+
+   LOGICAL  ::  ln_iceini    ! initialization or not
    !!----------------------------------------------------------------------
    !!   LIM 3.0,  UCL-LOCEAN-IPSL (2008)
@@ -87,5 +86,5 @@
       !! * Local variables
       INTEGER    :: ji, jj, jk, jl             ! dummy loop indices
-      REAL(wp)   :: epsi20, ztmelts, zdh
+      REAL(wp)   :: ztmelts, zdh
       INTEGER    :: i_hemis, i_fill, jl0  
       REAL(wp)   :: ztest_1, ztest_2, ztest_3, ztest_4, ztests, zsigma, zarg, zA, zV, zA_cons, zV_cons, zconv
@@ -101,6 +100,4 @@
       CALL wrk_alloc(   2,      zht_i_ini, zat_i_ini, zvt_i_ini, zht_s_ini, zsm_i_ini, ztm_i_ini )
 
-      epsi20   = 1.e-20_wp
-
       IF(lwp) WRITE(numout,*)
       IF(lwp) WRITE(numout,*) 'lim_istate : Ice initialization '
@@ -115,12 +112,12 @@
       ! surface temperature
       DO jl = 1, jpl ! loop over categories
-         t_su  (:,:,jl) = rtt * tms(:,:)
-         tn_ice(:,:,jl) = rtt * tms(:,:)
+         t_su  (:,:,jl) = rt0 * tmask(:,:,1)
+         tn_ice(:,:,jl) = rt0 * tmask(:,:,1)
       END DO
 
       ! basal temperature (considered at freezing point)
-      t_bo(:,:) = ( eos_fzp( tsn(:,:,1,jp_sal) ) + rt0 ) * tms(:,:) 
-
-      IF( ln_limini ) THEN
+      t_bo(:,:) = ( eos_fzp( sss_m(:,:) ) + rt0 ) * tmask(:,:,1) 
+
+      IF( ln_iceini ) THEN
 
       !--------------------------------------------------------------------
@@ -130,8 +127,8 @@
       DO jj = 1, jpj                                       ! ice if sst <= t-freez + ttest
          DO ji = 1, jpi
-            IF( ( tsn(ji,jj,1,jp_tem)  - ( t_bo(ji,jj) - rt0 ) ) * tms(ji,jj) >= thres_sst ) THEN 
-               zswitch(ji,jj) = 0._wp * tms(ji,jj)    ! no ice
+            IF( ( sst_m(ji,jj)  - ( t_bo(ji,jj) - rt0 ) ) * tmask(ji,jj,1) >= rn_thres_sst ) THEN 
+               zswitch(ji,jj) = 0._wp * tmask(ji,jj,1)    ! no ice
             ELSE                                                                                   
-               zswitch(ji,jj) = 1._wp * tms(ji,jj)    !    ice
+               zswitch(ji,jj) = 1._wp * tmask(ji,jj,1)    !    ice
             ENDIF
          END DO
@@ -158,9 +155,9 @@
       !-----------------------------
       ! assign initial thickness, concentration, snow depth and salinity to an hemisphere-dependent array
-      zht_i_ini(1) = hti_ini_n ; zht_i_ini(2) = hti_ini_s  ! ice thickness
-      zht_s_ini(1) = hts_ini_n ; zht_s_ini(2) = hts_ini_s  ! snow depth
-      zat_i_ini(1) = ati_ini_n ; zat_i_ini(2) = ati_ini_s  ! ice concentration
-      zsm_i_ini(1) = smi_ini_n ; zsm_i_ini(2) = smi_ini_s  ! bulk ice salinity
-      ztm_i_ini(1) = tmi_ini_n ; ztm_i_ini(2) = tmi_ini_s  ! temperature (ice and snow)
+      zht_i_ini(1) = rn_hti_ini_n ; zht_i_ini(2) = rn_hti_ini_s  ! ice thickness
+      zht_s_ini(1) = rn_hts_ini_n ; zht_s_ini(2) = rn_hts_ini_s  ! snow depth
+      zat_i_ini(1) = rn_ati_ini_n ; zat_i_ini(2) = rn_ati_ini_s  ! ice concentration
+      zsm_i_ini(1) = rn_smi_ini_n ; zsm_i_ini(2) = rn_smi_ini_s  ! bulk ice salinity
+      ztm_i_ini(1) = rn_tmi_ini_n ; ztm_i_ini(2) = rn_tmi_ini_s  ! temperature (ice and snow)
 
       zvt_i_ini(:) = zht_i_ini(:) * zat_i_ini(:)   ! ice volume
@@ -197,11 +194,11 @@
                !--- Ice thicknesses in the i_fill - 1 first categories
                DO jl = 1, i_fill - 1
-                  zh_i_ini(jl,i_hemis)    = 0.5 * ( hi_max(jl) + hi_max(jl-1) )
+                  zh_i_ini(jl,i_hemis) = hi_mean(jl)
                END DO
                
                !--- jl0: most likely index where cc will be maximum
                DO jl = 1, jpl
-                  IF ( ( zht_i_ini(i_hemis) .GT. hi_max(jl-1) ) .AND. &
-                     ( zht_i_ini(i_hemis) .LE. hi_max(jl)   ) ) THEN
+                  IF ( ( zht_i_ini(i_hemis) >  hi_max(jl-1) ) .AND. &
+                     & ( zht_i_ini(i_hemis) <= hi_max(jl)   ) ) THEN
                      jl0 = jl
                   ENDIF
@@ -267,5 +264,5 @@
 
             ! Test 3: thickness of the last category is in-bounds ?
-            IF ( zh_i_ini(i_fill, i_hemis) .GT. hi_max(i_fill-1) ) THEN
+            IF ( zh_i_ini(i_fill, i_hemis) > hi_max(i_fill-1) ) THEN
                ztest_3 = 1
             ELSE
@@ -317,9 +314,9 @@
             DO ji = 1, jpi
                a_i(ji,jj,jl)   = zswitch(ji,jj) * za_i_ini (jl,zhemis(ji,jj))  ! concentration
-               ht_i(ji,jj,jl)  = zswitch(ji,jj) * zh_i_ini(jl,zhemis(ji,jj))  ! ice thickness
+               ht_i(ji,jj,jl)  = zswitch(ji,jj) * zh_i_ini(jl,zhemis(ji,jj))   ! ice thickness
                ht_s(ji,jj,jl)  = ht_i(ji,jj,jl) * ( zht_s_ini( zhemis(ji,jj) ) / zht_i_ini( zhemis(ji,jj) ) )  ! snow depth
-               sm_i(ji,jj,jl)  = zswitch(ji,jj) * zsm_i_ini(zhemis(ji,jj)) !+ ( 1._wp - zswitch(ji,jj) ) * s_i_min ! salinity
-               o_i(ji,jj,jl)   = zswitch(ji,jj) * 1._wp + ( 1._wp - zswitch(ji,jj) ) ! age
-               t_su(ji,jj,jl)  = zswitch(ji,jj) * ztm_i_ini(zhemis(ji,jj)) + ( 1._wp - zswitch(ji,jj) ) * rtt ! surf temp
+               sm_i(ji,jj,jl)  = zswitch(ji,jj) * zsm_i_ini(zhemis(ji,jj))     ! salinity
+               o_i(ji,jj,jl)   = zswitch(ji,jj) * 1._wp                        ! age (1 day)
+               t_su(ji,jj,jl)  = zswitch(ji,jj) * ztm_i_ini(zhemis(ji,jj)) + ( 1._wp - zswitch(ji,jj) ) * rt0 ! surf temp
 
                ! This case below should not be used if (ht_s/ht_i) is ok in namelist
@@ -329,5 +326,5 @@
                ! recompute ht_i, ht_s avoiding out of bounds values
                ht_i(ji,jj,jl) = MIN( hi_max(jl), ht_i(ji,jj,jl) + zdh )
-               ht_s(ji,jj,jl) = MAX( 0._wp, ht_s(ji,jj,jl) - zdh * rhoic / rhosn )
+               ht_s(ji,jj,jl) = MAX( 0._wp, ht_s(ji,jj,jl) - zdh * rhoic * r1_rhosn )
 
                ! ice volume, salt content, age content
@@ -336,7 +333,7 @@
                smv_i(ji,jj,jl) = MIN( sm_i(ji,jj,jl) , sss_m(ji,jj) ) * v_i(ji,jj,jl) ! salt content
                oa_i(ji,jj,jl)  = o_i(ji,jj,jl) * a_i(ji,jj,jl)               ! age content
-            END DO ! ji
-         END DO ! jj
-      END DO ! jl
+            END DO
+         END DO
+      END DO
 
       ! Snow temperature and heat content
@@ -345,15 +342,14 @@
             DO jj = 1, jpj
                DO ji = 1, jpi
-                   t_s(ji,jj,jk,jl) = zswitch(ji,jj) * ztm_i_ini(zhemis(ji,jj)) + ( 1._wp - zswitch(ji,jj) ) * rtt
+                   t_s(ji,jj,jk,jl) = zswitch(ji,jj) * ztm_i_ini(zhemis(ji,jj)) + ( 1._wp - zswitch(ji,jj) ) * rt0
                    ! Snow energy of melting
-                   e_s(ji,jj,jk,jl) = zswitch(ji,jj) * rhosn * ( cpic * ( rtt - t_s(ji,jj,jk,jl) ) + lfus )
-                   ! Change dimensions
-                   e_s(ji,jj,jk,jl) = e_s(ji,jj,jk,jl) / unit_fac
-                   ! Multiply by volume, so that heat content in Joules
-                   e_s(ji,jj,jk,jl) = e_s(ji,jj,jk,jl) * area(ji,jj) * v_s(ji,jj,jl) / nlay_s
-               END DO ! ji
-            END DO ! jj
-         END DO ! jl
-      END DO ! jk
+                   e_s(ji,jj,jk,jl) = zswitch(ji,jj) * rhosn * ( cpic * ( rt0 - t_s(ji,jj,jk,jl) ) + lfus )
+
+                   ! Mutliply by volume, and divide by number of layers to get heat content in J/m2
+                   e_s(ji,jj,jk,jl) = e_s(ji,jj,jk,jl) * v_s(ji,jj,jl) * r1_nlay_s
+               END DO
+            END DO
+         END DO
+      END DO
 
       ! Ice salinity, temperature and heat content
@@ -362,27 +358,24 @@
             DO jj = 1, jpj
                DO ji = 1, jpi
-                   t_i(ji,jj,jk,jl) = zswitch(ji,jj) * ztm_i_ini(zhemis(ji,jj)) + ( 1._wp - zswitch(ji,jj) ) * rtt 
-                   s_i(ji,jj,jk,jl) = zswitch(ji,jj) * zsm_i_ini(zhemis(ji,jj)) !+ ( 1._wp - zswitch(ji,jj) ) * s_i_min
-                   ztmelts          = - tmut * s_i(ji,jj,jk,jl) + rtt !Melting temperature in K
+                   t_i(ji,jj,jk,jl) = zswitch(ji,jj) * ztm_i_ini(zhemis(ji,jj)) + ( 1._wp - zswitch(ji,jj) ) * rt0 
+                   s_i(ji,jj,jk,jl) = zswitch(ji,jj) * zsm_i_ini(zhemis(ji,jj)) !+ ( 1._wp - zswitch(ji,jj) ) * rn_simin
+                   ztmelts          = - tmut * s_i(ji,jj,jk,jl) + rt0 !Melting temperature in K
 
                    ! heat content per unit volume
                    e_i(ji,jj,jk,jl) = zswitch(ji,jj) * rhoic * (   cpic    * ( ztmelts - t_i(ji,jj,jk,jl) ) &
-                      +   lfus    * ( 1._wp - (ztmelts-rtt) / MIN((t_i(ji,jj,jk,jl)-rtt),-epsi20) ) &
-                      -   rcp     * ( ztmelts - rtt ) )
-
-                   ! Correct dimensions to avoid big values
-                   e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) / unit_fac 
-
-                   ! Mutliply by ice volume, and divide by number of layers to get heat content in J
-                   e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * area(ji,jj) * v_i(ji,jj,jl) / nlay_i
-               END DO ! ji
-            END DO ! jj
-         END DO ! jl
-      END DO ! jk
+                      +   lfus    * ( 1._wp - (ztmelts-rt0) / MIN((t_i(ji,jj,jk,jl)-rt0),-epsi20) ) &
+                      -   rcp     * ( ztmelts - rt0 ) )
+
+                   ! Mutliply by ice volume, and divide by number of layers to get heat content in J/m2
+                   e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * v_i(ji,jj,jl) * r1_nlay_i
+               END DO
+            END DO
+         END DO
+      END DO
 
       tn_ice (:,:,:) = t_su (:,:,:)
 
       ELSE 
-         ! if ln_limini=false
+         ! if ln_iceini=false
          a_i  (:,:,:) = 0._wp
          v_i  (:,:,:) = 0._wp
@@ -400,12 +393,12 @@
          DO jl = 1, jpl
             DO jk = 1, nlay_i
-               t_i(:,:,jk,jl) = rtt * tms(:,:)
+               t_i(:,:,jk,jl) = rt0 * tmask(:,:,1)
             END DO
             DO jk = 1, nlay_s
-               t_s(:,:,jk,jl) = rtt * tms(:,:)
+               t_s(:,:,jk,jl) = rt0 * tmask(:,:,1)
             END DO
          END DO
       
-      ENDIF ! ln_limini
+      ENDIF ! ln_iceini
       
       at_i (:,:) = 0.0_wp
@@ -481,6 +474,6 @@
       !!  8.5  ! 07-11 (M. Vancoppenolle) rewritten initialization
       !!-----------------------------------------------------------------------------
-      NAMELIST/namiceini/ ln_limini, thres_sst, hts_ini_n, hts_ini_s, hti_ini_n, hti_ini_s,  &
-         &                                      ati_ini_n, ati_ini_s, smi_ini_n, smi_ini_s, tmi_ini_n, tmi_ini_s
+      NAMELIST/namiceini/ ln_iceini, rn_thres_sst, rn_hts_ini_n, rn_hts_ini_s, rn_hti_ini_n, rn_hti_ini_s,  &
+         &                                      rn_ati_ini_n, rn_ati_ini_s, rn_smi_ini_n, rn_smi_ini_s, rn_tmi_ini_n, rn_tmi_ini_s
       INTEGER :: ios                 ! Local integer output status for namelist read
       !!-----------------------------------------------------------------------------
@@ -502,16 +495,16 @@
          WRITE(numout,*) 'lim_istate_init : ice parameters inititialisation '
          WRITE(numout,*) '~~~~~~~~~~~~~~~'
-         WRITE(numout,*) '   initialization with ice (T) or not (F)       ln_limini   = ', ln_limini
-         WRITE(numout,*) '   threshold water temp. for initial sea-ice    thres_sst  = ', thres_sst
-         WRITE(numout,*) '   initial snow thickness in the north          hts_ini_n  = ', hts_ini_n
-         WRITE(numout,*) '   initial snow thickness in the south          hts_ini_s  = ', hts_ini_s 
-         WRITE(numout,*) '   initial ice thickness  in the north          hti_ini_n  = ', hti_ini_n
-         WRITE(numout,*) '   initial ice thickness  in the south          hti_ini_s  = ', hti_ini_s
-         WRITE(numout,*) '   initial ice concentr.  in the north          ati_ini_n  = ', ati_ini_n
-         WRITE(numout,*) '   initial ice concentr.  in the north          ati_ini_s  = ', ati_ini_s
-         WRITE(numout,*) '   initial  ice salinity  in the north          smi_ini_n  = ', smi_ini_n
-         WRITE(numout,*) '   initial  ice salinity  in the south          smi_ini_s  = ', smi_ini_s
-         WRITE(numout,*) '   initial  ice/snw temp  in the north          tmi_ini_n  = ', tmi_ini_n
-         WRITE(numout,*) '   initial  ice/snw temp  in the south          tmi_ini_s  = ', tmi_ini_s
+         WRITE(numout,*) '   initialization with ice (T) or not (F)       ln_iceini     = ', ln_iceini
+         WRITE(numout,*) '   threshold water temp. for initial sea-ice    rn_thres_sst  = ', rn_thres_sst
+         WRITE(numout,*) '   initial snow thickness in the north          rn_hts_ini_n  = ', rn_hts_ini_n
+         WRITE(numout,*) '   initial snow thickness in the south          rn_hts_ini_s  = ', rn_hts_ini_s 
+         WRITE(numout,*) '   initial ice thickness  in the north          rn_hti_ini_n  = ', rn_hti_ini_n
+         WRITE(numout,*) '   initial ice thickness  in the south          rn_hti_ini_s  = ', rn_hti_ini_s
+         WRITE(numout,*) '   initial ice concentr.  in the north          rn_ati_ini_n  = ', rn_ati_ini_n
+         WRITE(numout,*) '   initial ice concentr.  in the north          rn_ati_ini_s  = ', rn_ati_ini_s
+         WRITE(numout,*) '   initial  ice salinity  in the north          rn_smi_ini_n  = ', rn_smi_ini_n
+         WRITE(numout,*) '   initial  ice salinity  in the south          rn_smi_ini_s  = ', rn_smi_ini_s
+         WRITE(numout,*) '   initial  ice/snw temp  in the north          rn_tmi_ini_n  = ', rn_tmi_ini_n
+         WRITE(numout,*) '   initial  ice/snw temp  in the south          rn_tmi_ini_s  = ', rn_tmi_ini_s
       ENDIF
 
