Index: branches/2012/dev_3352_UKMO8_CICE/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_cice.F90
===================================================================
--- branches/2012/dev_3352_UKMO8_CICE/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_cice.F90	(revision 3354)
+++ branches/2012/dev_3352_UKMO8_CICE/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_cice.F90	(revision 3355)
@@ -146,5 +146,5 @@
       !!---------------------------------------------------------------------
 
-      INTEGER  ::   ji, jj, jpl                        ! dummy loop indices
+      INTEGER  ::   ji, jj, jl                        ! dummy loop indices
 
       IF( nn_timing == 1 )  CALL timing_start('cice_sbc_init')
@@ -182,6 +182,6 @@
       CALL cice2nemo(aice,fr_i, 'T', 1. )
       IF ( (nsbc == 2).OR.(nsbc == 5) ) THEN
-         DO jpl=1,ncat
-            CALL cice2nemo(aicen(:,:,jpl,:),a_i(:,:,jpl), 'T', 1. )
+         DO jl=1,ncat
+            CALL cice2nemo(aicen(:,:,jl,:),a_i(:,:,jl), 'T', 1. )
          ENDDO
       ENDIF
@@ -212,5 +212,5 @@
       INTEGER, INTENT(in   ) ::   nsbc ! surface forcing type
 
-      INTEGER  ::   ji, jj, jpl                   ! dummy loop indices      
+      INTEGER  ::   ji, jj, jl                   ! dummy loop indices      
       REAL(wp), DIMENSION(:,:), POINTER :: ztmp
       REAL(wp), DIMENSION(:,:,:), POINTER :: ztmpn
@@ -259,6 +259,6 @@
 ! Surface downward latent heat flux (CI_5)
          IF (nsbc == 2) THEN
-            DO jpl=1,ncat
-               ztmpn(:,:,jpl)=qla_ice(:,:,1)*a_i(:,:,jpl)
+            DO jl=1,ncat
+               ztmpn(:,:,jl)=qla_ice(:,:,1)*a_i(:,:,jl)
             ENDDO
          ELSE
@@ -269,12 +269,12 @@
                DO ji=1,jpi
                   IF (fr_i(ji,jj).eq.0.0) THEN
-                     DO jpl=1,ncat
-                        ztmpn(ji,jj,jpl)=0.0
+                     DO jl=1,ncat
+                        ztmpn(ji,jj,jl)=0.0
                      ENDDO
                      ! This will then be conserved in CICE
                      ztmpn(ji,jj,1)=qla_ice(ji,jj,1)
                   ELSE
-                     DO jpl=1,ncat
-                        ztmpn(ji,jj,jpl)=qla_ice(ji,jj,1)*a_i(ji,jj,jpl)/fr_i(ji,jj)
+                     DO jl=1,ncat
+                        ztmpn(ji,jj,jl)=qla_ice(ji,jj,1)*a_i(ji,jj,jl)/fr_i(ji,jj)
                      ENDDO
                   ENDIF
@@ -282,24 +282,24 @@
             ENDDO
          ENDIF
-         DO jpl=1,ncat
-            CALL nemo2cice(ztmpn(:,:,jpl),flatn_f(:,:,jpl,:),'T', 1. )
+         DO jl=1,ncat
+            CALL nemo2cice(ztmpn(:,:,jl),flatn_f(:,:,jl,:),'T', 1. )
 
 ! GBM conductive flux through ice (CI_6)
 !  Convert to GBM
             IF (nsbc == 2) THEN
-               ztmp(:,:) = botmelt(:,:,jpl)*a_i(:,:,jpl)
+               ztmp(:,:) = botmelt(:,:,jl)*a_i(:,:,jl)
             ELSE
-               ztmp(:,:) = botmelt(:,:,jpl)
+               ztmp(:,:) = botmelt(:,:,jl)
             ENDIF
-            CALL nemo2cice(ztmp,fcondtopn_f(:,:,jpl,:),'T', 1. )
+            CALL nemo2cice(ztmp,fcondtopn_f(:,:,jl,:),'T', 1. )
 
 ! GBM surface heat flux (CI_7)
 !  Convert to GBM
             IF (nsbc == 2) THEN
-               ztmp(:,:) = (topmelt(:,:,jpl)+botmelt(:,:,jpl))*a_i(:,:,jpl) 
+               ztmp(:,:) = (topmelt(:,:,jl)+botmelt(:,:,jl))*a_i(:,:,jl) 
             ELSE
-               ztmp(:,:) = (topmelt(:,:,jpl)+botmelt(:,:,jpl))
+               ztmp(:,:) = (topmelt(:,:,jl)+botmelt(:,:,jl))
             ENDIF
-            CALL nemo2cice(ztmp,fsurfn_f(:,:,jpl,:),'T', 1. )
+            CALL nemo2cice(ztmp,fsurfn_f(:,:,jl,:),'T', 1. )
          ENDDO
 
@@ -420,11 +420,11 @@
       INTEGER, INTENT( in  ) ::   nsbc ! surface forcing type
       
-      INTEGER  ::   ji, jj, jpl                 ! dummy loop indices
-      REAL(wp), DIMENSION(:,:), POINTER :: ztmp
+      INTEGER  ::   ji, jj, jl                 ! dummy loop indices
+      REAL(wp), DIMENSION(:,:), POINTER :: ztmp1, ztmp2
       !!---------------------------------------------------------------------
 
       IF( nn_timing == 1 )  CALL timing_start('cice_sbc_out')
       !
-      CALL wrk_alloc( jpi,jpj, ztmp )
+      CALL wrk_alloc( jpi,jpj, ztmp1, ztmp2 )
       
       IF( kt == nit000 )  THEN
@@ -433,10 +433,10 @@
       
 ! x comp of ocean-ice stress 
-      CALL cice2nemo(strocnx,ztmp,'F', -1. )
+      CALL cice2nemo(strocnx,ztmp1,'F', -1. )
       ss_iou(:,:)=0.0
 ! F point to U point
       DO jj=2,jpjm1
          DO ji=2,jpim1
-            ss_iou(ji,jj) = 0.5 * ( ztmp(ji,jj-1) + ztmp(ji,jj) ) * umask(ji,jj,1)
+            ss_iou(ji,jj) = 0.5 * ( ztmp1(ji,jj-1) + ztmp1(ji,jj) ) * umask(ji,jj,1)
          ENDDO
       ENDDO
@@ -444,5 +444,5 @@
 
 ! y comp of ocean-ice stress 
-      CALL cice2nemo(strocny,ztmp,'F', -1. )
+      CALL cice2nemo(strocny,ztmp1,'F', -1. )
       ss_iov(:,:)=0.0
 ! F point to V point
@@ -450,5 +450,5 @@
       DO jj=1,jpjm1
          DO ji=2,jpim1
-            ss_iov(ji,jj) = 0.5 * ( ztmp(ji-1,jj) + ztmp(ji,jj) ) * vmask(ji,jj,1)
+            ss_iov(ji,jj) = 0.5 * ( ztmp1(ji-1,jj) + ztmp1(ji,jj) ) * vmask(ji,jj,1)
          ENDDO
       ENDDO
@@ -477,30 +477,31 @@
       ENDIF
 
-! Subtract fluxes from CICE to get freshwater equivalent flux used in 
-! salinity calculation
-      CALL cice2nemo(fresh_gbm,ztmp,'T', 1. )
-      emps(:,:)=emp(:,:)-ztmp(:,:)
-! Note the 1000.0 is to convert from kg salt to g salt (needed for PSU)
-      CALL cice2nemo(fsalt_gbm,ztmp,'T', 1. )
-      DO jj=1,jpj
-         DO ji=1,jpi
-            IF (sss_m(ji,jj).gt.0.0) THEN
-               emps(ji,jj)=emps(ji,jj)+ztmp(ji,jj)*1000.0/sss_m(ji,jj)
-            ENDIF
-         ENDDO
-      ENDDO
-
-! No longer remove precip over ice from free surface calculation on basis that the
+      CALL cice2nemo(fresh_gbm,ztmp1,'T', 1. )
+      CALL cice2nemo(fsalt_gbm,ztmp2,'T', 1. )
+! Check to avoid unphysical expression when ice is forming (ztmp1 negative)
+! Otherwise we are effectively allowing ice of higher salinity than the ocean to form
+! which has to be compensated for by the ocean salinity potentially going negative
+! This check breaks conservation but seems reasonable until we have prognostic ice salinity
+! The lines below ensure that when ice is forming emps will not be negative
+! Note the 1000.0 below is to convert from kg salt to g salt (needed for PSU)
+      emps(:,:)=0.0
+      WHERE (ztmp1(:,:).lt.0.0) ztmp2(:,:)=MAX(ztmp2(:,:),ztmp1(:,:)*sss_m(:,:)/1000.0)
+      WHERE (sss_m(:,:).gt.0.0) emps(:,:)=-ztmp1(:,:)+ztmp2(:,:)*1000.0/sss_m(:,:)
+! Now add other non-ice freshwater contributions into emps
+      emps(:,:)=emp(:,:)+emps(:,:)
+ 
+      IF (lk_vvl) THEN
+
+! The relevant quantity for calculating the salinity change is emps-emp
+! This is in fact just the ztmp2*1000.0/sss_m term above
+         emp(:,:)=emp(:,:)-ztmp1(:,:)
+
+      ELSE
+
+! Don't remove precip over ice from free surface calculation on basis that the
 ! weight of the precip will affect the free surface even if it falls on the ice
-! (same to the argument that freezing / melting of ice doesn't change the free surface) 
-! Sublimation from the ice is treated in a similar way (included in emp but not emps)  
-!
-! This should not be done in the variable volume case
-
-      IF (.NOT. lk_vvl) THEN
-
+! (same argument that freezing / melting of ice doesn't change the free surface) 
          emp(:,:)  = emp(:,:) - tprecip(:,:)*fr_i(:,:)
-
-! Take sublimation into account
+! Sublimation from the ice is treated in a similar way (included in emp but not emps)
          IF (nsbc == 5 ) THEN 
             emp(:,:) = emp(:,:) + ( emp_ice(:,:) + sprecip(:,:) )
@@ -532,6 +533,6 @@
 ! Now add in ice / snow related terms
 ! [fswthru will be zero unless running with calc_Tsfc=T in CICE]
-      CALL cice2nemo(fswthru_gbm,ztmp,'T', 1. )
-      qsr(:,:)=qsr(:,:)+ztmp(:,:)
+      CALL cice2nemo(fswthru_gbm,ztmp1,'T', 1. )
+      qsr(:,:)=qsr(:,:)+ztmp1(:,:)
       CALL lbc_lnk( qsr , 'T', 1. )
 
@@ -542,6 +543,6 @@
       ENDDO
 
-      CALL cice2nemo(fhocn_gbm,ztmp,'T', 1. )
-      qns(:,:)=qns(:,:)+nfrzmlt(:,:)+ztmp(:,:)
+      CALL cice2nemo(fhocn_gbm,ztmp1,'T', 1. )
+      qns(:,:)=qns(:,:)+nfrzmlt(:,:)+ztmp1(:,:)
 
       CALL lbc_lnk( qns , 'T', 1. )
@@ -551,6 +552,6 @@
       CALL cice2nemo(aice,fr_i,'T', 1. )
       IF ( (nsbc == 2).OR.(nsbc == 5) ) THEN
-         DO jpl=1,ncat
-            CALL cice2nemo(aicen(:,:,jpl,:),a_i(:,:,jpl), 'T', 1. )
+         DO jl=1,ncat
+            CALL cice2nemo(aicen(:,:,jl,:),a_i(:,:,jl), 'T', 1. )
          ENDDO
       ENDIF
@@ -570,5 +571,5 @@
 ! Release work space
 
-      CALL wrk_dealloc( jpi,jpj, ztmp )
+      CALL wrk_dealloc( jpi,jpj, ztmp1, ztmp2 )
       !
       IF( nn_timing == 1 )  CALL timing_stop('cice_sbc_out')
@@ -587,5 +588,5 @@
       !!---------------------------------------------------------------------
 
-      INTEGER  ::   jpl                        ! dummy loop index
+      INTEGER  ::   jl                        ! dummy loop index
       INTEGER  ::   ierror
 
@@ -610,7 +611,7 @@
 ! Snow and ice thicknesses (CO_2 and CO_3)
 
-      DO jpl = 1,ncat
-         CALL cice2nemo(vsnon(:,:,jpl,:),ht_s(:,:,jpl),'T', 1. )
-         CALL cice2nemo(vicen(:,:,jpl,:),ht_i(:,:,jpl),'T', 1. )
+      DO jl = 1,ncat
+         CALL cice2nemo(vsnon(:,:,jl,:),ht_s(:,:,jl),'T', 1. )
+         CALL cice2nemo(vicen(:,:,jl,:),ht_i(:,:,jl),'T', 1. )
       ENDDO
       !
