Index: /branches/UKMO/dev_merge_2017_CICE_interface/NEMOGCM/NEMO/OPA_SRC/SBC/sbc_ice.F90
===================================================================
--- /branches/UKMO/dev_merge_2017_CICE_interface/NEMOGCM/NEMO/OPA_SRC/SBC/sbc_ice.F90	(revision 9818)
+++ /branches/UKMO/dev_merge_2017_CICE_interface/NEMOGCM/NEMO/OPA_SRC/SBC/sbc_ice.F90	(revision 9819)
@@ -38,4 +38,10 @@
    LOGICAL         , PUBLIC, PARAMETER ::   lk_cice    = .TRUE.   !: CICE ice model
    CHARACTER(len=1), PUBLIC            ::   cp_ice_msh = 'F'      !: 'F'-grid ice-velocity
+   
+   ! Parameters imported from LIM to get CICE to work with NEMO4
+   INTEGER , PUBLIC, PARAMETER ::   np_jules_OFF    = 0  !: no Jules coupling (ice thermodynamics forced via qsr and qns)
+   INTEGER , PUBLIC, PARAMETER ::   np_jules_EMULE  = 1  !: emulated Jules coupling via icethd_zdf.F90 (BL99) (1st round compute qcn and qsr_tr, 2nd round use it)
+   INTEGER , PUBLIC, PARAMETER ::   np_jules_ACTIVE = 2  !: active Jules coupling                      (SM0L) (compute qcn and qsr_tr via sbcblk.F90 or sbccpl.F90)
+   INTEGER , PUBLIC, PARAMETER ::   nice_jules = np_jules_ACTIVE   !: Choice of jules coupling 
 # endif
 
@@ -56,5 +62,5 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   emp_ice        !: sublimation - precip over sea ice          [kg/m2/s]
 
-#if defined  key_lim3
+#if defined  key_lim3 || defined key_cice
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   evap_ice       !: sublimation                              [kg/m2/s]
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   devap_ice      !: sublimation sensitivity                [kg/m2/s/K]
Index: /branches/UKMO/dev_merge_2017_CICE_interface/NEMOGCM/NEMO/OPA_SRC/SBC/sbccpl.F90
===================================================================
--- /branches/UKMO/dev_merge_2017_CICE_interface/NEMOGCM/NEMO/OPA_SRC/SBC/sbccpl.F90	(revision 9818)
+++ /branches/UKMO/dev_merge_2017_CICE_interface/NEMOGCM/NEMO/OPA_SRC/SBC/sbccpl.F90	(revision 9819)
@@ -1246,5 +1246,5 @@
       !                                                      !   ice skin temp.   !
       !                                                      ! ================== !
-#if defined key_lim3
+#if defined key_lim3 || defined key_cice
       ! needed by Met Office
       IF( srcv(jpr_ts_ice)%laction ) THEN 
@@ -1715,5 +1715,6 @@
       REAL(wp), DIMENSION(jpi,jpj)     ::   zemp_tot, zemp_ice, zemp_oce, ztprecip, zsprecip  , zevap_oce, zdevap_ice
       REAL(wp), DIMENSION(jpi,jpj)     ::   zqns_tot, zqns_oce, zqsr_tot, zqsr_oce, zqprec_ice, zqemp_oce, zqemp_ice
-      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zqns_ice, zqsr_ice, zdqns_ice, zqevap_ice, zevap_ice    !!gm , zfrqsr_tr_i
+      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zqns_ice, zqsr_ice, zdqns_ice, zqevap_ice    !!gm , zfrqsr_tr_i
+      REAL(wp), DIMENSION(jpi,jpj)     ::   ztmp
       !!----------------------------------------------------------------------
       !
@@ -1735,5 +1736,35 @@
          ztprecip(:,:) =   frcv(jpr_rain)%z3(:,:,1) + zsprecip(:,:)  ! May need to ensure positive here
          zemp_tot(:,:) =   frcv(jpr_tevp)%z3(:,:,1) - ztprecip(:,:)
+         
+#if defined key_cice
+         IF ( TRIM(sn_rcv_emp%clcat) == 'yes' ) THEN
+            ! zemp_ice is the sum of frcv(jpr_ievp)%z3(:,:,1) over all layers - snow
+            zemp_ice(:,:) = - frcv(jpr_snow)%z3(:,:,1)
+            DO jl=1,jpl
+               zemp_ice(:,:   ) = zemp_ice(:,:) + frcv(jpr_ievp)%z3(:,:,jl)
+            ENDDO
+            ! latent heat coupled for each category in CICE
+            zevap_ice(:,:,1:jpl) = - frcv(jpr_ievp)%z3(:,:,1:jpl)
+         ELSE
+            ! If CICE has multicategories it still expects coupling fields for
+            ! each even if we treat as a single field
+            ! The latent heat flux is split between the ice categories according
+            ! to the fraction of the ice in each category
+            zemp_ice(:,:) = frcv(jpr_ievp)%z3(:,:,1) - frcv(jpr_snow)%z3(:,:,1)
+            WHERE ( ziceld(:,:) /= 0._wp ) 
+               ztmp(:,:) = 1./ziceld(:,:)
+            ELSEWHERE 
+               ztmp(:,:) = 0.e0
+            END WHERE  
+            DO jl=1,jpl
+               zevap_ice(:,:,jl) = - a_i(:,:,jl) * ztmp(:,:) * frcv(jpr_ievp)%z3(:,:,1)
+            END DO
+            WHERE ( ziceld(:,:) == 0._wp )  zevap_ice(:,:,1) = -frcv(jpr_ievp)%z3(:,:,1)
+         ENDIF
+         
+#else             
          zemp_ice(:,:) = ( frcv(jpr_ievp)%z3(:,:,1) - frcv(jpr_snow)%z3(:,:,1) ) * picefr(:,:)
+#endif
+
       CASE( 'oce and ice'   )   ! received fields: jpr_sbpr, jpr_semp, jpr_oemp, jpr_ievp
          zemp_tot(:,:) = ziceld(:,:) * frcv(jpr_oemp)%z3(:,:,1) + picefr(:,:) * frcv(jpr_sbpr)%z3(:,:,1)
@@ -1961,5 +1992,5 @@
       ELSE
          qns_tot(:,:  ) = zqns_tot(:,:  )
-         qns_ice(:,:,:) = zqns_ice(:,:,:)
+         IF ( TRIM( sn_rcv_qsr%cldes ) /= 'oce only' ) qns_ice(:,:,:) = zqns_ice(:,:,:)
       ENDIF
 
@@ -2043,7 +2074,5 @@
       ELSE
          qsr_tot(:,:  ) = zqsr_tot(:,:  )
-         IF ( TRIM( sn_rcv_qsr%cldes ) /= 'oce only' ) THEN
-            qsr_ice(:,:,:) = zqsr_ice(:,:,:)
-         ENDIF
+         IF ( TRIM( sn_rcv_qsr%cldes ) /= 'oce only' ) qsr_ice(:,:,:) = zqsr_ice(:,:,:)
       ENDIF
 
@@ -2070,5 +2099,5 @@
       ENDIF
 
-#if defined key_lim3      
+#if defined key_lim3 || defined key_cice
       !                                                      ! ========================= !
       SELECT CASE( TRIM( sn_rcv_iceflx%cldes ) )             !  ice topmelt and botmelt  !
Index: /branches/UKMO/dev_merge_2017_CICE_interface/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_cice.F90
===================================================================
--- /branches/UKMO/dev_merge_2017_CICE_interface/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_cice.F90	(revision 9818)
+++ /branches/UKMO/dev_merge_2017_CICE_interface/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_cice.F90	(revision 9819)
@@ -209,5 +209,5 @@
          ! there is no restart file.
          ! Values from a CICE restart file would overwrite this
-         CALL nemo2cice( tsn(:,:,1,jp_tem) , sst , 'T' , 1.) 
+         CALL nemo2cice( tsn(:,:,1,jp_tem) , sst , 'T' , 1.,'sst1') 
 #endif
 
@@ -216,5 +216,5 @@
       ! calculate surface freezing temperature and send to CICE
       CALL  eos_fzp(sss_m(:,:), sstfrz(:,:), gdept_n(:,:,1)) 
-      CALL nemo2cice(sstfrz,Tf, 'T', 1. )
+      CALL nemo2cice(sstfrz,Tf, 'T', 1.,'tf1' )
 
       CALL cice2nemo(aice,fr_i, 'T', 1. )
@@ -320,4 +320,5 @@
          IF(lwp) WRITE(numout,*)'cice_sbc_in'
       ENDIF
+      WRITE(numout,*)'zevap_ice 7: min, max = ', MINVAL(zevap_ice(:,:,:)), MAXVAL(zevap_ice(:,:,:))
 
       ztmp(:,:)=0.0
@@ -340,5 +341,6 @@
             ENDDO
          ENDDO
-         CALL nemo2cice(ztmp,strax,'F', -1. )
+         
+         CALL nemo2cice(ztmp,strax,'F', -1.,'strax' )         
 
 ! y comp of wind stress (CI_2)
@@ -350,5 +352,5 @@
             ENDDO
          ENDDO
-         CALL nemo2cice(ztmp,stray,'F', -1. )
+         CALL nemo2cice(ztmp,stray,'F', -1.,'stray' )
 
 
@@ -379,5 +381,5 @@
 
          DO jl=1,ncat
-            CALL nemo2cice(ztmpn(:,:,jl),flatn_f(:,:,jl,:),'T', 1. )
+            CALL nemo2cice(ztmpn(:,:,jl),flatn_f(:,:,jl,:),'T', 1.,'flatn' )
 
 ! GBM conductive flux through ice (CI_6)
@@ -388,5 +390,5 @@
                ztmp(:,:) = qcn_ice(:,:,jl)
             ENDIF
-            CALL nemo2cice(ztmp,fcondtopn_f(:,:,jl,:),'T', 1. )
+            CALL nemo2cice(ztmp,fcondtopn_f(:,:,jl,:),'T', 1.,'fcondtopn' )
 
 ! GBM surface heat flux (CI_7)
@@ -397,5 +399,5 @@
                ztmp(:,:) = (qml_ice(:,:,jl)+qcn_ice(:,:,jl))
             ENDIF
-            CALL nemo2cice(ztmp,fsurfn_f(:,:,jl,:),'T', 1. )
+            CALL nemo2cice(ztmp,fsurfn_f(:,:,jl,:),'T', 1. ,'fsurfn_f')
          ENDDO
 
@@ -405,24 +407,24 @@
 ! x comp and y comp of atmosphere surface wind (CICE expects on T points)
          ztmp(:,:) = wndi_ice(:,:)
-         CALL nemo2cice(ztmp,uatm,'T', -1. )
+         CALL nemo2cice(ztmp,uatm,'T', -1. ,'uatm')
          ztmp(:,:) = wndj_ice(:,:)
-         CALL nemo2cice(ztmp,vatm,'T', -1. )
+         CALL nemo2cice(ztmp,vatm,'T', -1.,'vatm' )
          ztmp(:,:) = SQRT ( wndi_ice(:,:)**2 + wndj_ice(:,:)**2 )
-         CALL nemo2cice(ztmp,wind,'T', 1. )    ! Wind speed (m/s)
+         CALL nemo2cice(ztmp,wind,'T', 1.,'wind' )    ! Wind speed (m/s)
          ztmp(:,:) = qsr_ice(:,:,1)
-         CALL nemo2cice(ztmp,fsw,'T', 1. )     ! Incoming short-wave (W/m^2)
+         CALL nemo2cice(ztmp,fsw,'T', 1.,'fsw' )     ! Incoming short-wave (W/m^2)
          ztmp(:,:) = qlw_ice(:,:,1)
-         CALL nemo2cice(ztmp,flw,'T', 1. )     ! Incoming long-wave (W/m^2)
+         CALL nemo2cice(ztmp,flw,'T', 1.,'flw' )     ! Incoming long-wave (W/m^2)
          ztmp(:,:) = tatm_ice(:,:)
-         CALL nemo2cice(ztmp,Tair,'T', 1. )    ! Air temperature (K)
-         CALL nemo2cice(ztmp,potT,'T', 1. )    ! Potential temp (K)
+         CALL nemo2cice(ztmp,Tair,'T', 1.,'tair' )    ! Air temperature (K)
+         CALL nemo2cice(ztmp,potT,'T', 1.,'pott' )    ! Potential temp (K)
 ! Following line uses MAX(....) to avoid problems if tatm_ice has unset halo rows  
          ztmp(:,:) = 101000. / ( 287.04 * MAX(1.0,tatm_ice(:,:)) )    
                                                ! Constant (101000.) atm pressure assumed
-         CALL nemo2cice(ztmp,rhoa,'T', 1. )    ! Air density (kg/m^3)
+         CALL nemo2cice(ztmp,rhoa,'T', 1.,'rhoa' )    ! Air density (kg/m^3)
          ztmp(:,:) = qatm_ice(:,:)
-         CALL nemo2cice(ztmp,Qa,'T', 1. )      ! Specific humidity (kg/kg)
+         CALL nemo2cice(ztmp,Qa,'T', 1.,'qa' )      ! Specific humidity (kg/kg)
          ztmp(:,:)=10.0
-         CALL nemo2cice(ztmp,zlvl,'T', 1. )    ! Atmos level height (m)
+         CALL nemo2cice(ztmp,zlvl,'T', 1.,'zlvl' )    ! Atmos level height (m)
 
 ! May want to check all values are physically realistic (as in CICE routine 
@@ -431,11 +433,11 @@
 ! Divide shortwave into spectral bands (as in prepare_forcing)
          ztmp(:,:)=qsr_ice(:,:,1)*frcvdr       ! visible direct
-         CALL nemo2cice(ztmp,swvdr,'T', 1. )             
+         CALL nemo2cice(ztmp,swvdr,'T', 1.,'swvdr' )             
          ztmp(:,:)=qsr_ice(:,:,1)*frcvdf       ! visible diffuse
-         CALL nemo2cice(ztmp,swvdf,'T', 1. )              
+         CALL nemo2cice(ztmp,swvdf,'T', 1.,'swvdf' )              
          ztmp(:,:)=qsr_ice(:,:,1)*frcidr       ! near IR direct
-         CALL nemo2cice(ztmp,swidr,'T', 1. )
+         CALL nemo2cice(ztmp,swidr,'T', 1.,'swidr' )
          ztmp(:,:)=qsr_ice(:,:,1)*frcidf       ! near IR diffuse
-         CALL nemo2cice(ztmp,swidf,'T', 1. )
+         CALL nemo2cice(ztmp,swidf,'T', 1. ,'swidf')
 
       ENDIF
@@ -444,5 +446,5 @@
 !Ice concentration change (from assimilation)
       ztmp(:,:)=ndaice_da(:,:)*tmask(:,:,1)
-      Call nemo2cice(ztmp,daice_da,'T', 1. )
+      Call nemo2cice(ztmp,daice_da,'T', 1. ,'daice_da')
 #endif 
 
@@ -451,29 +453,29 @@
       IF( iom_use('snowpre') )   CALL iom_put('snowpre',MAX( (1.0-fr_i(:,:))*sprecip(:,:) ,0.0)) !!Joakim edit  
       ztmp(:,:)=MAX(fr_i(:,:)*sprecip(:,:),0.0)  
-      CALL nemo2cice(ztmp,fsnow,'T', 1. ) 
+      CALL nemo2cice(ztmp,fsnow,'T', 1.,'fsnow' ) 
 
 ! Rainfall
       IF( iom_use('precip') )   CALL iom_put('precip', (1.0-fr_i(:,:))*(tprecip(:,:)-sprecip(:,:)) ) !!Joakim edit
       ztmp(:,:)=fr_i(:,:)*(tprecip(:,:)-sprecip(:,:))
-      CALL nemo2cice(ztmp,frain,'T', 1. ) 
+      CALL nemo2cice(ztmp,frain,'T', 1.,'frain' ) 
 
 ! Recalculate freezing temperature and send to CICE 
       CALL eos_fzp(sss_m(:,:), sstfrz(:,:), gdept_n(:,:,1)) 
-      CALL nemo2cice(sstfrz,Tf,'T', 1. )
+      CALL nemo2cice(sstfrz,Tf,'T', 1. ,'Tf')
 
 ! Freezing/melting potential
 ! Calculated over NEMO leapfrog timestep (hence 2*dt)
       nfrzmlt(:,:)=rau0*rcp*e3t_m(:,:)*(sstfrz(:,:)-sst_m(:,:))/(2.0*dt) 
-      CALL nemo2cice(nfrzmlt,frzmlt,'T', 1. )
+      CALL nemo2cice(nfrzmlt,frzmlt,'T', 1.,'frzmlt' )
 
 ! SST  and SSS
 
-      CALL nemo2cice(sst_m,sst,'T', 1. )
-      CALL nemo2cice(sss_m,sss,'T', 1. )
+      CALL nemo2cice(sst_m,sst,'T', 1.,'sst' )
+      CALL nemo2cice(sss_m,sss,'T', 1.,'sss' )
 
       IF( ksbc == jp_purecpl ) THEN
 ! Sea ice surface skin temperature
          DO jl=1,ncat
-           CALL nemo2cice(tsfc_ice(:,:,jl), trcrn(:,:,nt_tsfc,jl,:),'T',1.)
+           CALL nemo2cice(tsfc_ice(:,:,jl), trcrn(:,:,nt_tsfc,jl,:),'T',1., 'trcrn')
          ENDDO 
       ENDIF
@@ -486,5 +488,5 @@
          ENDDO
       ENDDO
-      CALL nemo2cice(ztmp,uocn,'F', -1. )
+      CALL nemo2cice(ztmp,uocn,'F', -1., 'uocn' )
 
 ! V point to F point
@@ -494,5 +496,5 @@
          ENDDO
       ENDDO
-      CALL nemo2cice(ztmp,vocn,'F', -1. )
+      CALL nemo2cice(ztmp,vocn,'F', -1. ,'vocn')
 
       IF( ln_ice_embd ) THEN             !== embedded sea ice: compute representative ice top surface ==!
@@ -521,5 +523,5 @@
          END DO
       END DO
-      CALL nemo2cice( ztmp,ss_tltx,'F', -1. )
+      CALL nemo2cice( ztmp,ss_tltx,'F', -1. , 'ss_tltx')
 
 ! T point to F point
@@ -530,5 +532,5 @@
          END DO
       END DO
-      CALL nemo2cice(ztmp,ss_tlty,'F', -1. )
+      CALL nemo2cice(ztmp,ss_tlty,'F', -1. ,'ss_tlty')
       !
    END SUBROUTINE cice_sbc_in
@@ -900,5 +902,5 @@
    END SUBROUTINE cice_sbc_force
 
-   SUBROUTINE nemo2cice( pn, pc, cd_type, psgn)
+   SUBROUTINE nemo2cice( pn, pc, cd_type, psgn, varname)
       !!---------------------------------------------------------------------
       !!                    ***  ROUTINE nemo2cice  ***
@@ -931,4 +933,5 @@
 #endif
       REAL (kind=dbl_kind), dimension(nx_block,ny_block,max_blocks) :: pc
+      CHARACTER(len=10), INTENT( in ) :: varname
       INTEGER (int_kind) :: &
          field_type,        &! id for type of field (scalar, vector, angle)
@@ -1010,4 +1013,8 @@
 #endif
 
+      IF ( ln_ctl ) THEN
+         WRITE(numout,*)'nemo2cice: ',varname,' min,max = ',MINVAL(pc(:,:,:)), MAXVAL(pc(:,:,:))
+         CALL flush(numout)
+      ENDIF
    END SUBROUTINE nemo2cice
 
