Index: NEMO/branches/UKMO/NEMO_4.0.2_ENHANCE-02_ISF_nemo/src/OCE/SBC/cpl_oasis3.F90
===================================================================
--- NEMO/branches/UKMO/NEMO_4.0.2_ENHANCE-02_ISF_nemo/src/OCE/SBC/cpl_oasis3.F90	(revision 12705)
+++ NEMO/branches/UKMO/NEMO_4.0.2_ENHANCE-02_ISF_nemo/src/OCE/SBC/cpl_oasis3.F90	(revision 12706)
@@ -306,7 +306,13 @@
       ! End of definition phase
       !------------------------------------------------------------------
-      
+      !     
+#if defined key_agrif
+      IF( agrif_fixed() == Agrif_Nb_Fine_Grids() ) THEN
+#endif
       CALL oasis_enddef(nerror)
       IF( nerror /= OASIS_Ok )   CALL oasis_abort ( ncomp_id, 'cpl_define', 'Failure in oasis_enddef')
+#if defined key_agrif
+      ENDIF
+#endif
       !
       IF ( ltmp_wapatch ) THEN
@@ -357,7 +363,7 @@
                      WRITE(numout,*) 'oasis_put:  kstep ', kstep
                      WRITE(numout,*) 'oasis_put:   info ', kinfo
-                     WRITE(numout,*) '     - Minimum value is ', MINVAL(pdata(:,:,jc))
-                     WRITE(numout,*) '     - Maximum value is ', MAXVAL(pdata(:,:,jc))
-                     WRITE(numout,*) '     -     Sum value is ', SUM(pdata(:,:,jc))
+                     WRITE(numout,*) '     - Minimum value is ', MINVAL(pdata(nldi:nlei,nldj:nlej,jc))
+                     WRITE(numout,*) '     - Maximum value is ', MAXVAL(pdata(nldi:nlei,nldj:nlej,jc))
+                     WRITE(numout,*) '     -     Sum value is ', SUM(pdata(nldi:nlei,nldj:nlej,jc))
                      WRITE(numout,*) '****************'
                   ENDIF
Index: NEMO/branches/UKMO/NEMO_4.0.2_ENHANCE-02_ISF_nemo/src/OCE/SBC/fldread.F90
===================================================================
--- NEMO/branches/UKMO/NEMO_4.0.2_ENHANCE-02_ISF_nemo/src/OCE/SBC/fldread.F90	(revision 12705)
+++ NEMO/branches/UKMO/NEMO_4.0.2_ENHANCE-02_ISF_nemo/src/OCE/SBC/fldread.F90	(revision 12706)
@@ -833,6 +833,6 @@
 
       REAL(wp), DIMENSION(:,:,:), INTENT(in   ) ::   pdta_read       ! data read in bdy file
-      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   pdta_read_z     ! depth of the data read in bdy file
-      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   pdta_read_dz    ! thickness of the levels in bdy file
+      REAL(wp), DIMENSION(:,:,:), INTENT(in   ) ::   pdta_read_z     ! depth of the data read in bdy file
+      REAL(wp), DIMENSION(:,:,:), INTENT(in   ) ::   pdta_read_dz    ! thickness of the levels in bdy file
       REAL(wp), DIMENSION(:,:,:), INTENT(  out) ::   pdta            ! output field on model grid (2 dimensional)
       REAL(wp)                  , INTENT(in   ) ::   pfv             ! fillvalue of the data read in bdy file
@@ -841,11 +841,11 @@
       INTEGER                   , INTENT(in   ) ::   kbdy            ! bdy number
       !!
-      INTEGER                                   ::   ipi             ! length of boundary data on local process
-      INTEGER                                   ::   ipkb            ! number of vertical levels in boundary data file
-      INTEGER                                   ::   jb, ji, jj, jk, jkb   ! loop counters
-      REAL(wp)                                  ::   zcoef
-      REAL(wp)                                  ::   zl, zi, zh      ! tmp variable for current depth and interpolation factor
-      REAL(wp)                                  ::   zfv_alt, ztrans, ztrans_new ! fillvalue and alternative -ABS(pfv)
-      REAL(wp), DIMENSION(jpk)                  ::   zdepth, zdhalf  ! level and half-level depth
+      INTEGER                  ::   ipi                 ! length of boundary data on local process
+      INTEGER                  ::   ipkb                ! number of vertical levels in boundary data file
+      INTEGER                  ::   ipkmax              ! number of vertical levels in boundary data file where no mask
+      INTEGER                  ::   jb, ji, jj, jk, jkb ! loop counters
+      REAL(wp)                 ::   zcoef, zi           ! 
+      REAL(wp)                 ::   ztrans, ztrans_new  ! transports
+      REAL(wp), DIMENSION(jpk) ::   zdepth, zdhalf      ! level and half-level depth
       !!---------------------------------------------------------------------
      
@@ -853,45 +853,20 @@
       ipkb = SIZE( pdta_read, 3 )
       
-      zfv_alt = -ABS(pfv)  ! set _FillValue < 0 as we make use of MAXVAL and MAXLOC later
-      !
-      WHERE( pdta_read == pfv )
-         pdta_read_z  = zfv_alt ! safety: put fillvalue into external depth field so consistent with data
-         pdta_read_dz = 0._wp   ! safety: put 0._wp into external thickness factors to ensure transport is correct
-      ENDWHERE
-      
       DO jb = 1, ipi
          ji = idx_bdy(kbdy)%nbi(jb,kgrd)
          jj = idx_bdy(kbdy)%nbj(jb,kgrd)
-         zh  = SUM(pdta_read_dz(jb,1,:) )
-         !
-         ! Warnings to flag differences in the input and model topgraphy - is this useful/necessary?
+         !
+         ! --- max jk where input data /= FillValue --- !
+         ipkmax = 1
+         DO jkb = 2, ipkb
+            IF( pdta_read(jb,1,jkb) /= pfv )   ipkmax = MAX( ipkmax, jkb )
+         END DO
+         !
+         ! --- calculate depth at t,u,v points --- !
          SELECT CASE( kgrd )                         
-         CASE(1)
-            IF( ABS( (zh - ht_n(ji,jj)) / ht_n(ji,jj)) * tmask(ji,jj,1) > 0.01_wp ) THEN
-               WRITE(ctmp1,"(I10.10)") jb 
-               CALL ctl_warn('fld_bdy_interp: T depths differ between grids at BDY point '//TRIM(ctmp1)//' by more than 1%')
-               !   IF(lwp) WRITE(numout,*) 'DEPTHT', zh, sum(e3t_n(ji,jj,:), mask=tmask(ji,jj,:)==1),  ht_n(ji,jj), jb, jb, ji, jj
-            ENDIF
-         CASE(2)
-            IF( ABS( (zh - hu_n(ji,jj)) * r1_hu_n(ji,jj)) * umask(ji,jj,1) > 0.01_wp ) THEN
-               WRITE(ctmp1,"(I10.10)") jb 
-               CALL ctl_warn('fld_bdy_interp: U depths differ between grids at BDY point '//TRIM(ctmp1)//' by more than 1%')
-               !   IF(lwp) WRITE(numout,*) 'DEPTHU', zh, SUM(e3u_n(ji,jj,:), mask=umask(ji,jj,:)==1),  SUM(umask(ji,jj,:)), &
-               !      &                hu_n(ji,jj), jb, jb, ji, jj, narea-1, pdta_read(jb,1,:)
-            ENDIF
-         CASE(3)
-            IF( ABS( (zh - hv_n(ji,jj)) * r1_hv_n(ji,jj)) * vmask(ji,jj,1) > 0.01_wp ) THEN
-               WRITE(ctmp1,"(I10.10)") jb
-               CALL ctl_warn('fld_bdy_interp: V depths differ between grids at BDY point '//TRIM(ctmp1)//' by more than 1%')
-            ENDIF
-         END SELECT
-         !
-         SELECT CASE( kgrd )                         
-         CASE(1)
-            ! depth of T points:
+         CASE(1)            ! depth of T points:
             zdepth(:) = gdept_n(ji,jj,:)
-         CASE(2)
-            ! depth of U points: we must not use gdept_n as we don't want to do a communication
-            !   --> copy what is done for gdept_n in domvvl...
+         CASE(2)            ! depth of U points: we must not use gdept_n as we don't want to do a communication
+            !                 --> copy what is done for gdept_n in domvvl...
             zdhalf(1) = 0.0_wp
             zdepth(1) = 0.5_wp * e3uw_n(ji,jj,1)
@@ -903,10 +878,9 @@
                zcoef = ( umask(ji,jj,jk) - wumask(ji,jj,jk) )
                zdhalf(jk) = zdhalf(jk-1) + e3u_n(ji,jj,jk-1)
-               zdepth(jk) =      zcoef  * ( zdhalf(jk  ) + 0.5 * e3uw_n(ji,jj,jk))  &
-                  &         + (1-zcoef) * ( zdepth(jk-1) +       e3uw_n(ji,jj,jk))
+               zdepth(jk) =       zcoef  * ( zdhalf(jk  ) + 0.5 * e3uw_n(ji,jj,jk))  &
+                  &         + (1.-zcoef) * ( zdepth(jk-1) +       e3uw_n(ji,jj,jk))
             END DO
-         CASE(3)
-            ! depth of V points: we must not use gdept_n as we don't want to do a communication
-            !   --> copy what is done for gdept_n in domvvl...
+         CASE(3)            ! depth of V points: we must not use gdept_n as we don't want to do a communication
+            !                 --> copy what is done for gdept_n in domvvl...
             zdhalf(1) = 0.0_wp
             zdepth(1) = 0.5_wp * e3vw_n(ji,jj,1)
@@ -918,76 +892,85 @@
                zcoef = ( vmask(ji,jj,jk) - wvmask(ji,jj,jk) )
                zdhalf(jk) = zdhalf(jk-1) + e3v_n(ji,jj,jk-1)
-               zdepth(jk) =      zcoef  * ( zdhalf(jk  ) + 0.5 * e3vw_n(ji,jj,jk))  &
-                  &         + (1-zcoef) * ( zdepth(jk-1) +       e3vw_n(ji,jj,jk))
+               zdepth(jk) =       zcoef  * ( zdhalf(jk  ) + 0.5 * e3vw_n(ji,jj,jk))  &
+                  &         + (1.-zcoef) * ( zdepth(jk-1) +       e3vw_n(ji,jj,jk))
             END DO
          END SELECT
-         !
-         DO jk = 1, jpk                      
-            IF(     zdepth(jk) < pdta_read_z(jb,1,          1) ) THEN                ! above the first level of external data
-               pdta(jb,1,jk) =  pdta_read(jb,1,1)
-            ELSEIF( zdepth(jk) > pdta_read_z(jb,1,ipkb) ) THEN                       ! below the last level of external data 
-               pdta(jb,1,jk) =  pdta_read(jb,1,MAXLOC(pdta_read_z(jb,1,:),1))
-            ELSE                                                             ! inbetween: vertical interpolation between jkb & jkb+1
-               DO jkb = 1, ipkb-1                                            ! when  gdept_n(jkb) < zdepth(jk) < gdept_n(jkb+1)
-                  IF( ( ( zdepth(jk) - pdta_read_z(jb,1,jkb) ) * ( zdepth(jk) - pdta_read_z(jb,1,jkb+1) ) <= 0._wp ) &
-                     &    .AND. ( pdta_read_z(jb,1,jkb+1) /= zfv_alt) ) THEN   ! linear interpolation between 2 levels
+         !         
+         ! --- interpolate bdy data on the model grid --- !
+         DO jk = 1, jpk
+            IF(     zdepth(jk) <= pdta_read_z(jb,1,1)      ) THEN   ! above the first level of external data
+               pdta(jb,1,jk) = pdta_read(jb,1,1)
+            ELSEIF( zdepth(jk) >  pdta_read_z(jb,1,ipkmax) ) THEN   ! below the last level of external data /= FillValue
+               pdta(jb,1,jk) = pdta_read(jb,1,ipkmax)
+            ELSE                                                    ! inbetween: vertical interpolation between jkb & jkb+1
+               DO jkb = 1, ipkmax-1
+                  IF( ( ( zdepth(jk) - pdta_read_z(jb,1,jkb) ) * ( zdepth(jk) - pdta_read_z(jb,1,jkb+1) ) ) <= 0._wp ) THEN ! linear interpolation between 2 levels
                      zi = ( zdepth(jk) - pdta_read_z(jb,1,jkb) ) / ( pdta_read_z(jb,1,jkb+1) - pdta_read_z(jb,1,jkb) )
-                     pdta(jb,1,jk) = pdta_read(jb,1,jkb) + ( pdta_read  (jb,1,jkb+1) - pdta_read  (jb,1,jkb) ) * zi
+                     pdta(jb,1,jk) = pdta_read(jb,1,jkb) + zi * ( pdta_read(jb,1,jkb+1) - pdta_read(jb,1,jkb) )
                   ENDIF
                END DO
             ENDIF
-         END DO   ! jpk
+         END DO
          !
       END DO   ! ipi
-      
-      IF(kgrd == 2) THEN                                  ! do we need to adjust the transport term?
+
+      ! --- mask data and adjust transport --- !
+      SELECT CASE( kgrd )                         
+
+      CASE(1)                                 ! mask data (probably unecessary)
          DO jb = 1, ipi
             ji = idx_bdy(kbdy)%nbi(jb,kgrd)
             jj = idx_bdy(kbdy)%nbj(jb,kgrd)
-            zh  = SUM(pdta_read_dz(jb,1,:) )
+            DO jk = 1, jpk                      
+               pdta(jb,1,jk) = pdta(jb,1,jk) * tmask(ji,jj,jk)
+            END DO
+         END DO
+         
+      CASE(2)                                 ! adjust the U-transport term
+         DO jb = 1, ipi
+            ji = idx_bdy(kbdy)%nbi(jb,kgrd)
+            jj = idx_bdy(kbdy)%nbj(jb,kgrd)
             ztrans = 0._wp
+            DO jkb = 1, ipkb                              ! calculate transport on input grid
+               IF( pdta_read(jb,1,jkb) /= pfv )   ztrans = ztrans + pdta_read(jb,1,jkb) * pdta_read_dz(jb,1,jkb)
+            ENDDO
             ztrans_new = 0._wp
-            DO jkb = 1, ipkb                              ! calculate transport on input grid
-               ztrans     = ztrans     + pdta_read(jb,1,jkb) * pdta_read_dz(jb, 1,jkb)
-            ENDDO
             DO jk = 1, jpk                                ! calculate transport on model grid
-               ztrans_new = ztrans_new +      pdta(jb,1,jk ) *        e3u_n(ji,jj,jk ) * umask(ji,jj,jk)
+               ztrans_new = ztrans_new + pdta(jb,1,jk ) * e3u_n(ji,jj,jk) * umask(ji,jj,jk)
             ENDDO
             DO jk = 1, jpk                                ! make transport correction
                IF(ldtotvel) THEN ! bdy data are total velocity so adjust bt transport term to match input data
                   pdta(jb,1,jk) = ( pdta(jb,1,jk) + ( ztrans - ztrans_new ) * r1_hu_n(ji,jj) ) * umask(ji,jj,jk)
-               ELSE ! we're just dealing with bc velocity so bt transport term should sum to zero
-                  pdta(jb,1,jk) =   pdta(jb,1,jk) + (  0._wp - ztrans_new ) * r1_hu_n(ji,jj)   * umask(ji,jj,jk)
+               ELSE              ! we're just dealing with bc velocity so bt transport term should sum to zero
+                  pdta(jb,1,jk) = ( pdta(jb,1,jk) + (  0._wp - ztrans_new ) * r1_hu_n(ji,jj) ) * umask(ji,jj,jk)
                ENDIF
             ENDDO
          ENDDO
-      ENDIF
-      
-      IF(kgrd == 3) THEN                                  ! do we need to adjust the transport term?
+
+      CASE(3)                                 ! adjust the V-transport term
          DO jb = 1, ipi
             ji = idx_bdy(kbdy)%nbi(jb,kgrd)
             jj = idx_bdy(kbdy)%nbj(jb,kgrd)
-            zh  = SUM(pdta_read_dz(jb,1,:) )
             ztrans = 0._wp
+            DO jkb = 1, ipkb                              ! calculate transport on input grid
+               IF( pdta_read(jb,1,jkb) /= pfv )   ztrans = ztrans + pdta_read(jb,1,jkb) * pdta_read_dz(jb,1,jkb)
+            ENDDO
             ztrans_new = 0._wp
-            DO jkb = 1, ipkb                              ! calculate transport on input grid
-               ztrans     = ztrans     + pdta_read(jb,1,jkb) * pdta_read_dz(jb, 1,jkb)
-            ENDDO
             DO jk = 1, jpk                                ! calculate transport on model grid
-               ztrans_new = ztrans_new +      pdta(jb,1,jk ) *        e3v_n(ji,jj,jk ) * vmask(ji,jj,jk)
+               ztrans_new = ztrans_new + pdta(jb,1,jk ) * e3v_n(ji,jj,jk) * vmask(ji,jj,jk)
             ENDDO
             DO jk = 1, jpk                                ! make transport correction
                IF(ldtotvel) THEN ! bdy data are total velocity so adjust bt transport term to match input data
                   pdta(jb,1,jk) = ( pdta(jb,1,jk) + ( ztrans - ztrans_new ) * r1_hv_n(ji,jj) ) * vmask(ji,jj,jk)
-               ELSE ! we're just dealing with bc velocity so bt transport term should sum to zero
-                  pdta(jb,1,jk) =   pdta(jb,1,jk) + (  0._wp - ztrans_new ) * r1_hv_n(ji,jj)   * vmask(ji,jj,jk)
+               ELSE              ! we're just dealing with bc velocity so bt transport term should sum to zero
+                  pdta(jb,1,jk) = ( pdta(jb,1,jk) + (  0._wp - ztrans_new ) * r1_hv_n(ji,jj) ) * vmask(ji,jj,jk)
                ENDIF
             ENDDO
          ENDDO
-      ENDIF
-      
+      END SELECT
+
    END SUBROUTINE fld_bdy_interp
 
-
+   
    SUBROUTINE fld_rot( kt, sd )
       !!---------------------------------------------------------------------
Index: NEMO/branches/UKMO/NEMO_4.0.2_ENHANCE-02_ISF_nemo/src/OCE/SBC/sbc_oce.F90
===================================================================
--- NEMO/branches/UKMO/NEMO_4.0.2_ENHANCE-02_ISF_nemo/src/OCE/SBC/sbc_oce.F90	(revision 12705)
+++ NEMO/branches/UKMO/NEMO_4.0.2_ENHANCE-02_ISF_nemo/src/OCE/SBC/sbc_oce.F90	(revision 12706)
@@ -104,5 +104,5 @@
    !!              Ocean Surface Boundary Condition fields
    !!----------------------------------------------------------------------
-   INTEGER , PUBLIC ::  ncpl_qsr_freq            !: qsr coupling frequency per days from atmosphere
+   INTEGER , PUBLIC ::  ncpl_qsr_freq = 0        !: qsr coupling frequency per days from atmosphere (used by top)
    !
    LOGICAL , PUBLIC ::   lhftau = .FALSE.        !: HF tau used in TKE: mean(stress module) - module(mean stress)
Index: NEMO/branches/UKMO/NEMO_4.0.2_ENHANCE-02_ISF_nemo/src/OCE/SBC/sbcblk.F90
===================================================================
--- NEMO/branches/UKMO/NEMO_4.0.2_ENHANCE-02_ISF_nemo/src/OCE/SBC/sbcblk.F90	(revision 12705)
+++ NEMO/branches/UKMO/NEMO_4.0.2_ENHANCE-02_ISF_nemo/src/OCE/SBC/sbcblk.F90	(revision 12706)
@@ -801,4 +801,5 @@
       REAL(wp), DIMENSION(jpi,jpj)     ::   zevap, zsnw   ! evaporation and snw distribution after wind blowing (SI3)
       REAL(wp), DIMENSION(jpi,jpj)     ::   zrhoa
+      REAL(wp), DIMENSION(jpi,jpj)     ::   ztmp, ztmp2
       !!---------------------------------------------------------------------
       !
@@ -913,4 +914,24 @@
          qtr_ice_top(:,:,:) = 0._wp 
       END WHERE
+      !
+
+      IF( iom_use('evap_ao_cea') .OR. iom_use('hflx_evap_cea') ) THEN
+         ztmp(:,:) = zevap(:,:) * ( 1._wp - at_i_b(:,:) ) 
+         CALL iom_put( 'evap_ao_cea'  , ztmp(:,:) * tmask(:,:,1) )   ! ice-free oce evap (cell average)
+         CALL iom_put( 'hflx_evap_cea', ztmp(:,:) * sst_m(:,:) * rcp * tmask(:,:,1) )   ! heat flux from evap (cell average)
+      ENDIF
+      IF( iom_use('hflx_rain_cea') ) THEN
+         ztmp(:,:) = rcp * ( SUM( (ptsu-rt0) * a_i_b, dim=3 ) + sst_m(:,:) * ( 1._wp - at_i_b(:,:) ) )
+         CALL iom_put( 'hflx_rain_cea', ( tprecip(:,:) - sprecip(:,:) ) * ztmp(:,:) )   ! heat flux from rain (cell average)
+      ENDIF
+      IF( iom_use('hflx_snow_cea') .OR. iom_use('hflx_snow_ao_cea') .OR. iom_use('hflx_snow_ai_cea')  )  THEN
+          WHERE( SUM( a_i_b, dim=3 ) > 1.e-10 ) ;   ztmp(:,:) = rcpi * SUM( (ptsu-rt0) * a_i_b, dim=3 ) / SUM( a_i_b, dim=3 )
+          ELSEWHERE                             ;   ztmp(:,:) = rcp * sst_m(:,:)    
+          ENDWHERE
+          ztmp2(:,:) = sprecip(:,:) * ( ztmp(:,:) - rLfus ) 
+          CALL iom_put('hflx_snow_cea'   , ztmp2(:,:) ) ! heat flux from snow (cell average)
+          CALL iom_put('hflx_snow_ao_cea', ztmp2(:,:) * ( 1._wp - zsnw(:,:) ) ) ! heat flux from snow (over ocean)
+          CALL iom_put('hflx_snow_ai_cea', ztmp2(:,:) *           zsnw(:,:)   ) ! heat flux from snow (over ice)
+      ENDIF
       !
       IF(ln_ctl) THEN
Index: NEMO/branches/UKMO/NEMO_4.0.2_ENHANCE-02_ISF_nemo/src/OCE/SBC/sbccpl.F90
===================================================================
--- NEMO/branches/UKMO/NEMO_4.0.2_ENHANCE-02_ISF_nemo/src/OCE/SBC/sbccpl.F90	(revision 12705)
+++ NEMO/branches/UKMO/NEMO_4.0.2_ENHANCE-02_ISF_nemo/src/OCE/SBC/sbccpl.F90	(revision 12706)
@@ -574,4 +574,10 @@
       IF ( TRIM( sn_rcv_emp%clcat    ) == 'yes' )   srcv(jpr_ievp)%nct       = nn_cats_cpl
 
+#if defined key_si3
+      IF( ln_cndflx .AND. .NOT.ln_cndemulate ) THEN 
+         IF( .NOT.srcv(jpr_ts_ice)%laction )  &
+            &   CALL ctl_stop( 'sbc_cpl_init: srcv(jpr_ts_ice)%laction should be set to true when ln_cndflx=T' )     
+      ENDIF
+#endif
       !                                                      ! ------------------------- !
       !                                                      !      Wave breaking        !    
@@ -863,5 +869,4 @@
       ELSE IF( sn_snd_crt%clvgrd /= 'T' ) THEN  
          CALL ctl_stop( 'sn_snd_crt%clvgrd must be equal to T' )
-         ssnd(jps_ocx1:jps_ivz1)%clgrid  = 'T'      ! all oce and ice components on the same unique grid
       ENDIF
       ssnd(jps_ocx1:jps_ivz1)%laction = .TRUE.   ! default: all are send
@@ -1041,9 +1046,4 @@
       ENDIF
       xcplmask(:,:,0) = 1. - SUM( xcplmask(:,:,1:nn_cplmodel), dim = 3 )
-      !
-      ncpl_qsr_freq = cpl_freq( 'O_QsrOce' ) + cpl_freq( 'O_QsrMix' ) + cpl_freq( 'I_QsrOce' ) + cpl_freq( 'I_QsrMix' )
-      IF( ln_dm2dc .AND. ln_cpl .AND. ncpl_qsr_freq /= 86400 )   &
-         &   CALL ctl_stop( 'sbc_cpl_init: diurnal cycle reconstruction (ln_dm2dc) needs daily couping for solar radiation' )
-      IF( ln_dm2dc .AND. ln_cpl ) ncpl_qsr_freq = 86400 / ncpl_qsr_freq
       !
    END SUBROUTINE sbc_cpl_init
@@ -1111,4 +1111,12 @@
       REAL(wp), DIMENSION(jpi,jpj) ::   ztx, zty, zmsk, zemp, zqns, zqsr
       !!----------------------------------------------------------------------
+      !
+      IF( kt == nit000 ) THEN
+      !   cannot be done in the init phase when we use agrif as cpl_freq requires that oasis_enddef is done
+         ncpl_qsr_freq = cpl_freq( 'O_QsrOce' ) + cpl_freq( 'O_QsrMix' ) + cpl_freq( 'I_QsrOce' ) + cpl_freq( 'I_QsrMix' )
+         IF( ln_dm2dc .AND. ncpl_qsr_freq /= 86400 )   &
+            &   CALL ctl_stop( 'sbc_cpl_rcv: diurnal cycle reconstruction (ln_dm2dc) needs daily couping for solar radiation' )
+         ncpl_qsr_freq = 86400 / ncpl_qsr_freq   ! used by top
+      ENDIF
       !
       IF( ln_mixcpl )   zmsk(:,:) = 1. - xcplmask(:,:,0)
@@ -1244,16 +1252,4 @@
       IF( srcv(jpr_co2)%laction )   atm_co2(:,:) = frcv(jpr_co2)%z3(:,:,1)
       !
-      !                                                      ! ================== !
-      !                                                      !   ice skin temp.   !
-      !                                                      ! ================== !
-#if defined key_si3
-      ! needed by Met Office
-      IF( srcv(jpr_ts_ice)%laction ) THEN 
-         WHERE    ( frcv(jpr_ts_ice)%z3(:,:,:) > 0.0  )   ;   tsfc_ice(:,:,:) = 0.0 
-         ELSEWHERE( frcv(jpr_ts_ice)%z3(:,:,:) < -60. )   ;   tsfc_ice(:,:,:) = -60.
-         ELSEWHERE                                        ;   tsfc_ice(:,:,:) = frcv(jpr_ts_ice)%z3(:,:,:)
-         END WHERE
-      ENDIF 
-#endif
       !                                                      ! ========================= ! 
       !                                                      ! Mean Sea Level Pressure   !   (taum) 
@@ -1635,11 +1631,11 @@
       !!                   sprecip           solid precipitation over the ocean  
       !!----------------------------------------------------------------------
-      REAL(wp), INTENT(in), DIMENSION(:,:)             ::   picefr     ! ice fraction                [0 to 1]
-      !                                                !!           ! optional arguments, used only in 'mixed oce-ice' case
-      REAL(wp), INTENT(in), DIMENSION(:,:,:), OPTIONAL ::   palbi      ! all skies ice albedo 
-      REAL(wp), INTENT(in), DIMENSION(:,:  ), OPTIONAL ::   psst       ! sea surface temperature     [Celsius]
-      REAL(wp), INTENT(in), DIMENSION(:,:,:), OPTIONAL ::   pist       ! ice surface temperature     [Kelvin]
-      REAL(wp), INTENT(in), DIMENSION(:,:,:), OPTIONAL ::   phs        ! snow depth                  [m]
-      REAL(wp), INTENT(in), DIMENSION(:,:,:), OPTIONAL ::   phi        ! ice thickness               [m]
+      REAL(wp), INTENT(in)   , DIMENSION(:,:)             ::   picefr     ! ice fraction                [0 to 1]
+      !                                                   !!           ! optional arguments, used only in 'mixed oce-ice' case or for Met-Office coupling
+      REAL(wp), INTENT(in)   , DIMENSION(:,:,:), OPTIONAL ::   palbi      ! all skies ice albedo 
+      REAL(wp), INTENT(in)   , DIMENSION(:,:  ), OPTIONAL ::   psst       ! sea surface temperature     [Celsius]
+      REAL(wp), INTENT(inout), DIMENSION(:,:,:), OPTIONAL ::   pist       ! ice surface temperature     [Kelvin] => inout for Met-Office
+      REAL(wp), INTENT(in)   , DIMENSION(:,:,:), OPTIONAL ::   phs        ! snow depth                  [m]
+      REAL(wp), INTENT(in)   , DIMENSION(:,:,:), OPTIONAL ::   phi        ! ice thickness               [m]
       !
       INTEGER  ::   ji, jj, jl   ! dummy loop index
@@ -1648,5 +1644,5 @@
       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, zevap_ice, zqtr_ice_top, ztsu
       !!----------------------------------------------------------------------
       !
@@ -1774,11 +1770,12 @@
       IF( srcv(jpr_cal)%laction )   CALL iom_put( 'calving_cea' , frcv(jpr_cal)%z3(:,:,1) * tmask(:,:,1)                )  ! calving
       IF( srcv(jpr_icb)%laction )   CALL iom_put( 'iceberg_cea' , frcv(jpr_icb)%z3(:,:,1) * tmask(:,:,1)                )  ! icebergs
-      IF( iom_use('snowpre') )      CALL iom_put( 'snowpre'     , sprecip(:,:)                                          )  ! Snow
-      IF( iom_use('precip') )       CALL iom_put( 'precip'      , tprecip(:,:)                                          )  ! total  precipitation
-      IF( iom_use('rain') )         CALL iom_put( 'rain'        , tprecip(:,:) - sprecip(:,:)                           )  ! liquid precipitation 
-      IF( iom_use('snow_ao_cea') )  CALL iom_put( 'snow_ao_cea' , sprecip(:,:) * ( 1._wp - zsnw(:,:) )                  )  ! Snow over ice-free ocean  (cell average)
-      IF( iom_use('snow_ai_cea') )  CALL iom_put( 'snow_ai_cea' , sprecip(:,:) *           zsnw(:,:)                    )  ! Snow over sea-ice         (cell average)
-      IF( iom_use('subl_ai_cea') )  CALL iom_put( 'subl_ai_cea' , frcv(jpr_ievp)%z3(:,:,1) * picefr(:,:) * tmask(:,:,1) )  ! Sublimation over sea-ice (cell average)
-      IF( iom_use('evap_ao_cea') )  CALL iom_put( 'evap_ao_cea' , ( frcv(jpr_tevp)%z3(:,:,1)  &
+      CALL iom_put( 'snowpre'     , sprecip(:,:)                                          )  ! Snow
+      CALL iom_put( 'precip'      , tprecip(:,:)                                          )  ! total  precipitation
+      IF ( iom_use('rain') ) CALL iom_put( 'rain'        , tprecip(:,:) - sprecip(:,:)                           )  ! liquid precipitation 
+      IF ( iom_use('snow_ao_cea') ) CALL iom_put( 'snow_ao_cea' , sprecip(:,:) * ( 1._wp - zsnw(:,:) )                  )  ! Snow over ice-free ocean  (cell average)
+      IF ( iom_use('snow_ai_cea') ) CALL iom_put( 'snow_ai_cea' , sprecip(:,:) *           zsnw(:,:)                    )  ! Snow over sea-ice         (cell average)
+      IF ( iom_use('rain_ao_cea') ) CALL iom_put( 'rain_ao_cea' , ( tprecip(:,:) - sprecip(:,:) ) * picefr(:,:)         )  ! liquid precipitation over ocean (cell average)
+      IF ( iom_use('subl_ai_cea') ) CALL iom_put( 'subl_ai_cea' , frcv(jpr_ievp)%z3(:,:,1) * picefr(:,:) * tmask(:,:,1) )  ! Sublimation over sea-ice (cell average)
+      IF ( iom_use('evap_ao_cea') ) CALL iom_put( 'evap_ao_cea' , ( frcv(jpr_tevp)%z3(:,:,1)  &
          &                                                        - frcv(jpr_ievp)%z3(:,:,1) * picefr(:,:) ) * tmask(:,:,1) )  ! ice-free oce evap (cell average)
       ! note: runoff output is done in sbcrnf (which includes icebergs too) and iceshelf output is done in sbcisf
@@ -1815,7 +1812,17 @@
 ! ** NEED TO SORT OUT HOW THIS SHOULD WORK IN THE MULTI-CATEGORY CASE - CURRENTLY NOT ALLOWED WHEN INTERFACE INITIALISED **
          zqns_tot(:,:  ) = frcv(jpr_qnsmix)%z3(:,:,1)
-         zqns_ice(:,:,1) = frcv(jpr_qnsmix)%z3(:,:,1)    &
-            &            + frcv(jpr_dqnsdt)%z3(:,:,1) * ( pist(:,:,1) - ( (rt0 + psst(:,:  ) ) * ziceld(:,:)   &
-            &                                           + pist(:,:,1) * picefr(:,:) ) )
+         IF ( TRIM(sn_rcv_qsr%clcat) == 'yes' ) THEN
+            DO jl = 1, jpl
+               zqns_ice(:,:,jl) = frcv(jpr_qnsmix)%z3(:,:,jl)    &
+                  &             + frcv(jpr_dqnsdt)%z3(:,:,jl) * ( pist(:,:,jl) - ( ( rt0 + psst(:,:) ) * ziceld(:,:)   &
+                  &                                                                + pist(:,:,jl) * picefr(:,:) ) )
+            END DO
+         ELSE
+            DO jl = 1, jpl
+               zqns_ice(:,:,jl) = frcv(jpr_qnsmix)%z3(:,:, 1)    &
+                  &             + frcv(jpr_dqnsdt)%z3(:,:, 1) * ( pist(:,:,jl) - ( ( rt0 + psst(:,:) ) * ziceld(:,:)   &
+                  &                                                                + pist(:,:,jl) * picefr(:,:) ) )
+            END DO
+         ENDIF
       END SELECT
       !                                     
@@ -1902,14 +1909,14 @@
 #endif
       ! outputs
-      IF ( srcv(jpr_cal)%laction       ) CALL iom_put('hflx_cal_cea'    , - frcv(jpr_cal)%z3(:,:,1) * rLfus )                      ! latent heat from calving
-      IF ( srcv(jpr_icb)%laction       ) CALL iom_put('hflx_icb_cea'    , - frcv(jpr_icb)%z3(:,:,1) * rLfus )                      ! latent heat from icebergs melting
-      IF ( iom_use('hflx_rain_cea')    ) CALL iom_put('hflx_rain_cea'   , ( tprecip(:,:) - sprecip(:,:) ) * zcptrain(:,:) )        ! heat flux from rain (cell average)
-      IF ( iom_use('hflx_evap_cea')    ) CALL iom_put('hflx_evap_cea'   , ( frcv(jpr_tevp)%z3(:,:,1) - frcv(jpr_ievp)%z3(:,:,1) &
-           &                                                              * picefr(:,:) ) * zcptn(:,:) * tmask(:,:,1) )            ! heat flux from evap (cell average)
-      IF ( iom_use('hflx_snow_cea')    ) CALL iom_put('hflx_snow_cea'   , sprecip(:,:) * ( zcptsnw(:,:) - rLfus )  )               ! heat flux from snow (cell average)
-      IF ( iom_use('hflx_snow_ao_cea') ) CALL iom_put('hflx_snow_ao_cea', sprecip(:,:) * ( zcptsnw(:,:) - rLfus ) &
-           &                                                              * ( 1._wp - zsnw(:,:) )                  )               ! heat flux from snow (over ocean)
-      IF ( iom_use('hflx_snow_ai_cea') ) CALL iom_put('hflx_snow_ai_cea', sprecip(:,:) * ( zcptsnw(:,:) - rLfus ) & 
-           &                                                              *           zsnw(:,:)                    )               ! heat flux from snow (over ice)
+      IF ( srcv(jpr_cal)%laction ) CALL iom_put('hflx_cal_cea'    , - frcv(jpr_cal)%z3(:,:,1) * rLfus )   ! latent heat from calving
+      IF ( srcv(jpr_icb)%laction ) CALL iom_put('hflx_icb_cea'    , - frcv(jpr_icb)%z3(:,:,1) * rLfus )   ! latent heat from icebergs melting
+      IF ( iom_use(   'hflx_rain_cea') ) CALL iom_put('hflx_rain_cea' , ( tprecip(:,:) - sprecip(:,:) ) * zcptrain(:,:) )        ! heat flux from rain (cell average)
+      IF ( iom_use(   'hflx_evap_cea') ) CALL iom_put('hflx_evap_cea' , ( frcv(jpr_tevp)%z3(:,:,1) - frcv(jpr_ievp)%z3(:,:,1) * picefr(:,:) )  &
+           &                         * zcptn(:,:) * tmask(:,:,1) )            ! heat flux from evap (cell average)
+      IF ( iom_use(   'hflx_prec_cea') ) CALL iom_put('hflx_prec_cea' ,    sprecip(:,:) * ( zcptsnw(:,:) - rLfus )  &                    ! heat flux from all precip (cell avg)
+         &                          + ( tprecip(:,:) - sprecip(:,:) ) * zcptrain(:,:) )
+      IF ( iom_use(   'hflx_snow_cea') ) CALL iom_put('hflx_snow_cea'   , sprecip(:,:) * ( zcptsnw(:,:) - rLfus )  )               ! heat flux from snow (cell average)
+      IF ( iom_use('hflx_snow_ao_cea') ) CALL iom_put('hflx_snow_ao_cea', sprecip(:,:) * ( zcptsnw(:,:) - rLfus ) * ( 1._wp - zsnw(:,:) ) )   ! heat flux from snow (over ocean)
+      IF ( iom_use('hflx_snow_ai_cea') ) CALL iom_put('hflx_snow_ai_cea', sprecip(:,:) * ( zcptsnw(:,:) - rLfus ) *  zsnw(:,:) )              ! heat flux from snow (over ice)
       ! note: hflx for runoff and iceshelf are done in sbcrnf and sbcisf resp.
       !
@@ -1929,6 +1936,4 @@
             END DO
          ENDIF
-         zqsr_tot(:,:  ) = frcv(jpr_qsrmix)%z3(:,:,1)
-         zqsr_ice(:,:,1) = frcv(jpr_qsrice)%z3(:,:,1)
       CASE( 'oce and ice' )
          zqsr_tot(:,:  ) =  ziceld(:,:) * frcv(jpr_qsroce)%z3(:,:,1)
@@ -1950,7 +1955,17 @@
 !       Create solar heat flux over ice using incoming solar heat flux and albedos
 !       ( see OASIS3 user guide, 5th edition, p39 )
-         zqsr_ice(:,:,1) = frcv(jpr_qsrmix)%z3(:,:,1) * ( 1.- palbi(:,:,1) )   &
-            &            / (  1.- ( alb_oce_mix(:,:  ) * ziceld(:,:)       &
-            &                     + palbi      (:,:,1) * picefr(:,:) ) )
+         IF ( TRIM(sn_rcv_qsr%clcat) == 'yes' ) THEN
+            DO jl = 1, jpl
+               zqsr_ice(:,:,jl) = frcv(jpr_qsrmix)%z3(:,:,jl) * ( 1.- palbi(:,:,jl) )   &
+                  &            / (  1.- ( alb_oce_mix(:,:   ) * ziceld(:,:)       &
+                  &                     + palbi      (:,:,jl) * picefr(:,:) ) )
+            END DO
+         ELSE
+            DO jl = 1, jpl
+               zqsr_ice(:,:,jl) = frcv(jpr_qsrmix)%z3(:,:, 1) * ( 1.- palbi(:,:,jl) )   &
+                  &            / (  1.- ( alb_oce_mix(:,:   ) * ziceld(:,:)       &
+                  &                     + palbi      (:,:,jl) * picefr(:,:) ) )
+            END DO
+         ENDIF
       CASE( 'none'      )       ! Not available as for now: needs additional coding  
       !                         ! since fields received, here zqsr_tot,  are not defined with none option
@@ -2012,8 +2027,14 @@
       !                                                      ! ========================= !
       CASE ('coupled')
-         qml_ice(:,:,:) = frcv(jpr_topm)%z3(:,:,:)
-         qcn_ice(:,:,:) = frcv(jpr_botm)%z3(:,:,:)
+         IF( ln_mixcpl ) THEN
+            DO jl=1,jpl
+               qml_ice(:,:,jl) = qml_ice(:,:,jl) * xcplmask(:,:,0) + frcv(jpr_topm)%z3(:,:,jl) * zmsk(:,:)
+               qcn_ice(:,:,jl) = qcn_ice(:,:,jl) * xcplmask(:,:,0) + frcv(jpr_botm)%z3(:,:,jl) * zmsk(:,:)
+            ENDDO
+         ELSE
+            qml_ice(:,:,:) = frcv(jpr_topm)%z3(:,:,:)
+            qcn_ice(:,:,:) = frcv(jpr_botm)%z3(:,:,:)
+         ENDIF
       END SELECT
-      !
       !                                                      ! ========================= !
       !                                                      !      Transmitted Qsr      !   [W/m2]
@@ -2022,9 +2043,13 @@
          !
          !                    ! ===> used prescribed cloud fraction representative for polar oceans in summer (0.81)
-         ztri = 0.18 * ( 1.0 - cldf_ice ) + 0.35 * cldf_ice    ! surface transmission parameter (Grenfell Maykut 77)
+         ztri = 0.18 * ( 1.0 - cldf_ice ) + 0.35 * cldf_ice    ! surface transmission when hi>10cm (Grenfell Maykut 77)
          !
-         qtr_ice_top(:,:,:) = ztri * qsr_ice(:,:,:)
-         WHERE( phs(:,:,:) >= 0.0_wp )   qtr_ice_top(:,:,:) = 0._wp            ! snow fully opaque
-         WHERE( phi(:,:,:) <= 0.1_wp )   qtr_ice_top(:,:,:) = qsr_ice(:,:,:)   ! thin ice transmits all solar radiation
+         WHERE    ( phs(:,:,:) <= 0._wp .AND. phi(:,:,:) <  0.1_wp )       ! linear decrease from hi=0 to 10cm  
+            zqtr_ice_top(:,:,:) = qsr_ice(:,:,:) * ( ztri + ( 1._wp - ztri ) * ( 1._wp - phi(:,:,:) * 10._wp ) )
+         ELSEWHERE( phs(:,:,:) <= 0._wp .AND. phi(:,:,:) >= 0.1_wp )       ! constant (ztri) when hi>10cm
+            zqtr_ice_top(:,:,:) = qsr_ice(:,:,:) * ztri
+         ELSEWHERE                                                         ! zero when hs>0
+            zqtr_ice_top(:,:,:) = 0._wp
+         END WHERE
          !     
       ELSEIF( ln_cndflx .AND. .NOT.ln_cndemulate ) THEN      !==  conduction flux as surface forcing  ==!
@@ -2032,5 +2057,32 @@
          !                    ! ===> here we must receive the qtr_ice_top array from the coupler
          !                           for now just assume zero (fully opaque ice)
-         qtr_ice_top(:,:,:) = 0._wp
+         zqtr_ice_top(:,:,:) = 0._wp
+         !
+      ENDIF
+      !
+      IF( ln_mixcpl ) THEN
+         DO jl=1,jpl
+            qtr_ice_top(:,:,jl) = qtr_ice_top(:,:,jl) * xcplmask(:,:,0) + zqtr_ice_top(:,:,jl) * zmsk(:,:)
+         ENDDO
+      ELSE
+         qtr_ice_top(:,:,:) = zqtr_ice_top(:,:,:)
+      ENDIF
+      !                                                      ! ================== !
+      !                                                      !   ice skin temp.   !
+      !                                                      ! ================== !
+      ! needed by Met Office
+      IF( srcv(jpr_ts_ice)%laction ) THEN 
+         WHERE    ( frcv(jpr_ts_ice)%z3(:,:,:) > 0.0  )   ;   ztsu(:,:,:) =   0. + rt0 
+         ELSEWHERE( frcv(jpr_ts_ice)%z3(:,:,:) < -60. )   ;   ztsu(:,:,:) = -60. + rt0
+         ELSEWHERE                                        ;   ztsu(:,:,:) = frcv(jpr_ts_ice)%z3(:,:,:) + rt0
+         END WHERE
+         !
+         IF( ln_mixcpl ) THEN
+            DO jl=1,jpl
+               pist(:,:,jl) = pist(:,:,jl) * xcplmask(:,:,0) + ztsu(:,:,jl) * zmsk(:,:)
+            ENDDO
+         ELSE
+            pist(:,:,:) = ztsu(:,:,:)
+         ENDIF
          !
       ENDIF
@@ -2195,5 +2247,5 @@
          CASE default    ;   CALL ctl_stop( 'sbc_cpl_snd: wrong definition of sn_snd_thick%clcat' )
          END SELECT
-         IF( ssnd(jps_fice)%laction )   CALL cpl_snd( jps_fice, isec, ztmp3, info )
+         CALL cpl_snd( jps_fice, isec, ztmp3, info )
       ENDIF
 
@@ -2255,5 +2307,5 @@
       !                                                      !      Ice melt ponds       ! 
       !                                                      ! ------------------------- !
-      ! needed by Met Office
+      ! needed by Met Office: 1) fraction of ponded ice 2) local/actual pond depth 
       IF( ssnd(jps_a_p)%laction .OR. ssnd(jps_ht_p)%laction ) THEN 
          SELECT CASE( sn_snd_mpnd%cldes)  
@@ -2261,12 +2313,12 @@
             SELECT CASE( sn_snd_mpnd%clcat )  
             CASE( 'yes' )  
-               ztmp3(:,:,1:jpl) =  a_ip(:,:,1:jpl)
-               ztmp4(:,:,1:jpl) =  v_ip(:,:,1:jpl)  
+               ztmp3(:,:,1:jpl) =  a_ip_frac(:,:,1:jpl)
+               ztmp4(:,:,1:jpl) =  h_ip(:,:,1:jpl)  
             CASE( 'no' )  
                ztmp3(:,:,:) = 0.0  
                ztmp4(:,:,:) = 0.0  
                DO jl=1,jpl  
-                 ztmp3(:,:,1) = ztmp3(:,:,1) + a_ip(:,:,jpl)  
-                 ztmp4(:,:,1) = ztmp4(:,:,1) + v_ip(:,:,jpl) 
+                 ztmp3(:,:,1) = ztmp3(:,:,1) + a_ip_frac(:,:,jpl)  
+                 ztmp4(:,:,1) = ztmp4(:,:,1) + h_ip(:,:,jpl) 
                ENDDO  
             CASE default   ;   CALL ctl_stop( 'sbc_cpl_snd: wrong definition of sn_snd_mpnd%clcat' )  
@@ -2306,5 +2358,8 @@
       !                                                      !  CO2 flux from PISCES     ! 
       !                                                      ! ------------------------- !
-      IF( ssnd(jps_co2)%laction .AND. l_co2cpl )   CALL cpl_snd( jps_co2, isec, RESHAPE ( oce_co2, (/jpi,jpj,1/) ) , info )
+      IF( ssnd(jps_co2)%laction .AND. l_co2cpl )   THEN  
+         ztmp1(:,:) = oce_co2(:,:) * 1000.  ! conversion in molC/m2/s 
+         CALL cpl_snd( jps_co2, isec, RESHAPE ( ztmp1, (/jpi,jpj,1/) ) , info ) 
+      ENDIF 
       !
       !                                                      ! ------------------------- !
Index: NEMO/branches/UKMO/NEMO_4.0.2_ENHANCE-02_ISF_nemo/src/OCE/SBC/sbcmod.F90
===================================================================
--- NEMO/branches/UKMO/NEMO_4.0.2_ENHANCE-02_ISF_nemo/src/OCE/SBC/sbcmod.F90	(revision 12705)
+++ NEMO/branches/UKMO/NEMO_4.0.2_ENHANCE-02_ISF_nemo/src/OCE/SBC/sbcmod.F90	(revision 12706)
@@ -236,4 +236,12 @@
 #endif
       !
+      !
+      IF( sbc_ssr_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'sbc_init : unable to allocate sbc_ssr arrays' )
+      IF( .NOT.ln_ssr ) THEN               !* Initialize qrp and erp if no restoring 
+         qrp(:,:) = 0._wp
+         erp(:,:) = 0._wp
+      ENDIF
+      !
+
       IF( nn_ice == 0 ) THEN        !* No sea-ice in the domain : ice fraction is always zero
          IF( nn_components /= jp_iam_opa )   fr_i(:,:) = 0._wp    ! except for OPA in SAS-OPA coupled case
@@ -536,4 +544,6 @@
          CALL iom_put( "taum"  , taum       )                   ! wind stress module
          CALL iom_put( "wspd"  , wndm       )                   ! wind speed  module over free ocean or leads in presence of sea-ice
+         CALL iom_put( "qrp", qrp )                             ! heat flux damping
+         CALL iom_put( "erp", erp )                             ! freshwater flux damping
       ENDIF
       !
Index: NEMO/branches/UKMO/NEMO_4.0.2_ENHANCE-02_ISF_nemo/src/OCE/SBC/sbcrnf.F90
===================================================================
--- NEMO/branches/UKMO/NEMO_4.0.2_ENHANCE-02_ISF_nemo/src/OCE/SBC/sbcrnf.F90	(revision 12705)
+++ NEMO/branches/UKMO/NEMO_4.0.2_ENHANCE-02_ISF_nemo/src/OCE/SBC/sbcrnf.F90	(revision 12706)
@@ -20,5 +20,5 @@
    USE sbc_oce        ! surface boundary condition variables
    USE eosbn2         ! Equation Of State
-   USE closea         ! closed seas
+   USE closea, ONLY: l_clo_rnf, clo_rnf ! closed seas
    !
    USE in_out_manager ! I/O manager
@@ -42,8 +42,10 @@
    REAL(wp)                   ::      rn_dep_max        !: depth over which runoffs is spread       (ln_rnf_depth_ini =T)
    INTEGER                    ::      nn_rnf_depth_file !: create (=1) a runoff depth file or not (=0)
+   LOGICAL                    ::   ln_rnf_icb        !: iceberg flux is specified in a file
    LOGICAL                    ::   ln_rnf_tem        !: temperature river runoffs attribute specified in a file
    LOGICAL           , PUBLIC ::   ln_rnf_sal        !: salinity    river runoffs attribute specified in a file
    TYPE(FLD_N)       , PUBLIC ::   sn_rnf            !: information about the runoff file to be read
    TYPE(FLD_N)                ::   sn_cnf            !: information about the runoff mouth file to be read
+   TYPE(FLD_N)                ::   sn_i_rnf        !: information about the iceberg flux file to be read
    TYPE(FLD_N)                ::   sn_s_rnf          !: information about the salinities of runoff file to be read
    TYPE(FLD_N)                ::   sn_t_rnf          !: information about the temperatures of runoff file to be read
@@ -64,4 +66,5 @@
 
    TYPE(FLD),        ALLOCATABLE, DIMENSION(:) ::   sf_rnf       ! structure: river runoff (file information, fields read)
+   TYPE(FLD),        ALLOCATABLE, DIMENSION(:) ::   sf_i_rnf     ! structure: iceberg flux (file information, fields read)
    TYPE(FLD),        ALLOCATABLE, DIMENSION(:) ::   sf_s_rnf     ! structure: river runoff salinity (file information, fields read)  
    TYPE(FLD),        ALLOCATABLE, DIMENSION(:) ::   sf_t_rnf     ! structure: river runoff temperature (file information, fields read)  
@@ -111,5 +114,9 @@
       !                                            !-------------------!
       !
-      IF( .NOT. l_rnfcpl )   CALL fld_read ( kt, nn_fsbc, sf_rnf   )    ! Read Runoffs data and provide it at kt
+      !
+      IF( .NOT. l_rnfcpl )  THEN
+                            CALL fld_read ( kt, nn_fsbc, sf_rnf   )    ! Read Runoffs data and provide it at kt ( runoffs + iceberg )
+         IF( ln_rnf_icb )   CALL fld_read ( kt, nn_fsbc, sf_i_rnf )    ! idem for iceberg flux if required
+      ENDIF
       IF(   ln_rnf_tem   )   CALL fld_read ( kt, nn_fsbc, sf_t_rnf )    ! idem for runoffs temperature if required
       IF(   ln_rnf_sal   )   CALL fld_read ( kt, nn_fsbc, sf_s_rnf )    ! idem for runoffs salinity    if required
@@ -117,5 +124,12 @@
       IF( MOD( kt - 1, nn_fsbc ) == 0 ) THEN
          !
-         IF( .NOT. l_rnfcpl )   rnf(:,:) = rn_rfact * ( sf_rnf(1)%fnow(:,:,1) ) * tmask(:,:,1)       ! updated runoff value at time step kt
+         IF( .NOT. l_rnfcpl ) THEN
+             rnf(:,:) = rn_rfact * ( sf_rnf(1)%fnow(:,:,1) ) * tmask(:,:,1)  ! updated runoff value at time step kt
+             IF( ln_rnf_icb ) THEN
+                fwficb(:,:) = rn_rfact * ( sf_i_rnf(1)%fnow(:,:,1) ) * tmask(:,:,1)  ! updated runoff value at time step kt
+                CALL iom_put( 'iceberg_cea'  , fwficb(:,:)  )         ! output iceberg flux
+                CALL iom_put( 'hflx_icb_cea' , fwficb(:,:) * rLfus )   ! output Heat Flux into Sea Water due to Iceberg Thermodynamics -->
+             ENDIF
+         ENDIF
          !
          !                                                           ! set temperature & salinity content of runoffs
@@ -128,10 +142,10 @@
          ELSE                                                        ! use SST as runoffs temperature
             !CEOD River is fresh water so must at least be 0 unless we consider ice
-            rnf_tsc(:,:,jp_tem) = MAX(sst_m(:,:),0.0_wp) * rnf(:,:) * r1_rau0
+            rnf_tsc(:,:,jp_tem) = MAX( sst_m(:,:), 0.0_wp ) * rnf(:,:) * r1_rau0
          ENDIF
          !                                                           ! use runoffs salinity data
          IF( ln_rnf_sal )   rnf_tsc(:,:,jp_sal) = ( sf_s_rnf(1)%fnow(:,:,1) ) * rnf(:,:) * r1_rau0
          !                                                           ! else use S=0 for runoffs (done one for all in the init)
-         IF( iom_use('runoffs') )        CALL iom_put( 'runoffs'     , rnf(:,:)                         )   ! output runoff mass flux
+                                         CALL iom_put( 'runoffs'     , rnf(:,:)                         )   ! output runoff mass flux
          IF( iom_use('hflx_rnf_cea') )   CALL iom_put( 'hflx_rnf_cea', rnf_tsc(:,:,jp_tem) * rau0 * rcp )   ! output runoff sensible heat (W/m2)
       ENDIF
@@ -238,6 +252,6 @@
       REAL(wp), DIMENSION(jpi,jpj,2) :: zrnfcl    
       !!
-      NAMELIST/namsbc_rnf/ cn_dir            , ln_rnf_depth, ln_rnf_tem, ln_rnf_sal,   &
-         &                 sn_rnf, sn_cnf    , sn_s_rnf    , sn_t_rnf  , sn_dep_rnf,   &
+      NAMELIST/namsbc_rnf/ cn_dir            , ln_rnf_depth, ln_rnf_tem, ln_rnf_sal, ln_rnf_icb,   &
+         &                 sn_rnf, sn_cnf    , sn_i_rnf, sn_s_rnf    , sn_t_rnf  , sn_dep_rnf,   &
          &                 ln_rnf_mouth      , rn_hrnf     , rn_avt_rnf, rn_rfact,     &
          &                 ln_rnf_depth_ini  , rn_dep_max  , rn_rnf_max, nn_rnf_depth_file
@@ -261,9 +275,9 @@
       !                                   ! ============
       !
-      REWIND( numnam_ref )              ! Namelist namsbc_rnf in reference namelist : Runoffs 
+      REWIND( numnam_ref )
       READ  ( numnam_ref, namsbc_rnf, IOSTAT = ios, ERR = 901)
 901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namsbc_rnf in reference namelist' )
 
-      REWIND( numnam_cfg )              ! Namelist namsbc_rnf in configuration namelist : Runoffs
+      REWIND( numnam_cfg )
       READ  ( numnam_cfg, namsbc_rnf, IOSTAT = ios, ERR = 902 )
 902   IF( ios >  0 )   CALL ctl_nam ( ios , 'namsbc_rnf in configuration namelist' )
@@ -295,4 +309,19 @@
          IF( sn_rnf%ln_tint ) ALLOCATE( sf_rnf(1)%fdta(jpi,jpj,1,2) )
          CALL fld_fill( sf_rnf, (/ sn_rnf /), cn_dir, 'sbc_rnf_init', 'read runoffs data', 'namsbc_rnf', no_print )
+         !
+         IF( ln_rnf_icb ) THEN                      ! Create (if required) sf_i_rnf structure
+            IF(lwp) WRITE(numout,*)
+            IF(lwp) WRITE(numout,*) '          iceberg flux read in a file'
+            ALLOCATE( sf_i_rnf(1), STAT=ierror  )
+            IF( ierror > 0 ) THEN
+               CALL ctl_stop( 'sbc_rnf_init: unable to allocate sf_i_rnf structure' )   ;   RETURN
+            ENDIF
+            ALLOCATE( sf_i_rnf(1)%fnow(jpi,jpj,1)   )
+            IF( sn_i_rnf%ln_tint ) ALLOCATE( sf_i_rnf(1)%fdta(jpi,jpj,1,2) )
+            CALL fld_fill (sf_i_rnf, (/ sn_i_rnf /), cn_dir, 'sbc_rnf_init', 'read iceberg flux data', 'namsbc_rnf' )
+         ELSE
+            fwficb(:,:) = 0._wp
+         ENDIF
+
       ENDIF
       !
@@ -337,5 +366,5 @@
                IF( h_rnf(ji,jj) > 0._wp ) THEN
                   jk = 2
-                  DO WHILE ( jk /= mbkt(ji,jj) .AND. gdept_0(ji,jj,jk) < h_rnf(ji,jj) ) ;  jk = jk + 1
+                  DO WHILE ( jk < mbkt(ji,jj) .AND. gdept_0(ji,jj,jk) < h_rnf(ji,jj) ) ;  jk = jk + 1
                   END DO
                   nk_rnf(ji,jj) = jk
@@ -394,5 +423,5 @@
                IF( zrnfcl(ji,jj,1) > 0._wp ) THEN
                   jk = 2
-                  DO WHILE ( jk /= mbkt(ji,jj) .AND. gdept_0(ji,jj,jk) < h_rnf(ji,jj) ) ;  jk = jk + 1
+                  DO WHILE ( jk < mbkt(ji,jj) .AND. gdept_0(ji,jj,jk) < h_rnf(ji,jj) ) ;  jk = jk + 1
                   END DO
                   nk_rnf(ji,jj) = jk
@@ -435,5 +464,5 @@
          !                                      !    - mixed upstream-centered (ln_traadv_cen2=T)
          !
-         IF ( ln_rnf_depth )   CALL ctl_warn( 'sbc_rnf_init: increased mixing turned on but effects may already',   &
+         IF( ln_rnf_depth )   CALL ctl_warn( 'sbc_rnf_init: increased mixing turned on but effects may already',   &
             &                                              'be spread through depth by ln_rnf_depth'               )
          !
Index: NEMO/branches/UKMO/NEMO_4.0.2_ENHANCE-02_ISF_nemo/src/OCE/SBC/sbcssr.F90
===================================================================
--- NEMO/branches/UKMO/NEMO_4.0.2_ENHANCE-02_ISF_nemo/src/OCE/SBC/sbcssr.F90	(revision 12705)
+++ NEMO/branches/UKMO/NEMO_4.0.2_ENHANCE-02_ISF_nemo/src/OCE/SBC/sbcssr.F90	(revision 12706)
@@ -30,7 +30,9 @@
    PUBLIC   sbc_ssr        ! routine called in sbcmod
    PUBLIC   sbc_ssr_init   ! routine called in sbcmod
+   PUBLIC   sbc_ssr_alloc  ! routine called in sbcmod
 
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   erp   !: evaporation damping   [kg/m2/s]
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   qrp   !: heat flux damping        [w/m2]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   coefice   !: under ice relaxation coefficient
 
    !                                   !!* Namelist namsbc_ssr *
@@ -41,4 +43,5 @@
    LOGICAL         ::   ln_sssr_bnd     ! flag to bound erp term 
    REAL(wp)        ::   rn_sssr_bnd     ! ABS(Max./Min.) value of erp term [mm/day]
+   INTEGER         ::   nn_sssr_ice     ! Control of restoring under ice
 
    REAL(wp) , ALLOCATABLE, DIMENSION(:) ::   buffer   ! Temporary buffer for exchange
@@ -97,5 +100,17 @@
                   END DO
                END DO
-               CALL iom_put( "qrp", qrp )                             ! heat flux damping
+            ENDIF
+            !
+            IF( nn_sssr /= 0 .AND. nn_sssr_ice /= 1 ) THEN
+              ! use fraction of ice ( fr_i ) to adjust relaxation under ice if nn_sssr_ice .ne. 1
+              ! n.b. coefice is initialised and fixed to 1._wp if nn_sssr_ice = 1
+               DO jj = 1, jpj
+                  DO ji = 1, jpi
+                     SELECT CASE ( nn_sssr_ice )
+                       CASE ( 0 )    ;  coefice(ji,jj) = 1._wp - fr_i(ji,jj)              ! no/reduced damping under ice
+                       CASE  DEFAULT ;  coefice(ji,jj) = 1._wp + ( nn_sssr_ice - 1 ) * fr_i(ji,jj) ! reinforced damping (x nn_sssr_ice) under ice )
+                     END SELECT
+                  END DO
+               END DO
             ENDIF
             !
@@ -105,4 +120,5 @@
                   DO ji = 1, jpi
                      zerp = zsrp * ( 1. - 2.*rnfmsk(ji,jj) )   &      ! No damping in vicinity of river mouths
+                        &        *   coefice(ji,jj)            &      ! Optional control of damping under sea-ice
                         &        * ( sss_m(ji,jj) - sf_sss(1)%fnow(ji,jj,1) ) * tmask(ji,jj,1)
                      sfx(ji,jj) = sfx(ji,jj) + zerp                 ! salt flux
@@ -110,5 +126,4 @@
                   END DO
                END DO
-               CALL iom_put( "erp", erp )                             ! freshwater flux damping
                !
             ELSEIF( nn_sssr == 2 ) THEN                               !* Salinity damping term (volume flux (emp) and associated heat flux (qns)
@@ -118,4 +133,5 @@
                   DO ji = 1, jpi                            
                      zerp = zsrp * ( 1. - 2.*rnfmsk(ji,jj) )   &      ! No damping in vicinity of river mouths
+                        &        *   coefice(ji,jj)            &      ! Optional control of damping under sea-ice
                         &        * ( sss_m(ji,jj) - sf_sss(1)%fnow(ji,jj,1) )   &
                         &        / MAX(  sss_m(ji,jj), 1.e-20   ) * tmask(ji,jj,1)
@@ -126,5 +142,4 @@
                   END DO
                END DO
-               CALL iom_put( "erp", erp )                             ! freshwater flux damping
             ENDIF
             !
@@ -154,5 +169,6 @@
       CHARACTER(len=100) ::  cn_dir          ! Root directory for location of ssr files
       TYPE(FLD_N) ::   sn_sst, sn_sss        ! informations about the fields to be read
-      NAMELIST/namsbc_ssr/ cn_dir, nn_sstr, nn_sssr, rn_dqdt, rn_deds, sn_sst, sn_sss, ln_sssr_bnd, rn_sssr_bnd
+      NAMELIST/namsbc_ssr/ cn_dir, nn_sstr, nn_sssr, rn_dqdt, rn_deds, sn_sst, &
+              & sn_sss, ln_sssr_bnd, rn_sssr_bnd, nn_sssr_ice
       INTEGER     ::  ios
       !!----------------------------------------------------------------------
@@ -182,9 +198,9 @@
          WRITE(numout,*) '         flag to bound erp term                 ln_sssr_bnd = ', ln_sssr_bnd
          WRITE(numout,*) '         ABS(Max./Min.) erp threshold           rn_sssr_bnd = ', rn_sssr_bnd, ' mm/day'
-      ENDIF
-      !
-      !                            !* Allocate erp and qrp array
-      ALLOCATE( qrp(jpi,jpj), erp(jpi,jpj), STAT=ierror )
-      IF( ierror > 0 )   CALL ctl_stop( 'STOP', 'sbc_ssr: unable to allocate erp and qrp array' )
+         WRITE(numout,*) '      Cntrl of surface restoration under ice nn_sssr_ice    = ', nn_sssr_ice
+         WRITE(numout,*) '          ( 0 = no restoration under ice)'
+         WRITE(numout,*) '          ( 1 = restoration everywhere  )'
+         WRITE(numout,*) '          (>1 = enhanced restoration under ice  )'
+      ENDIF
       !
       IF( nn_sstr == 1 ) THEN      !* set sf_sst structure & allocate arrays
@@ -216,4 +232,5 @@
       ENDIF
       !
+      coefice(:,:) = 1._wp         !  Initialise coefice to 1._wp ; will not need to be changed if nn_sssr_ice=1
       !                            !* Initialize qrp and erp if no restoring 
       IF( nn_sstr /= 1                   )   qrp(:,:) = 0._wp
@@ -221,4 +238,18 @@
       !
    END SUBROUTINE sbc_ssr_init
+         
+   INTEGER FUNCTION sbc_ssr_alloc()
+      !!----------------------------------------------------------------------
+      !!               ***  FUNCTION sbc_ssr_alloc  ***
+      !!----------------------------------------------------------------------
+      sbc_ssr_alloc = 0       ! set to zero if no array to be allocated
+      IF( .NOT. ALLOCATED( erp ) ) THEN
+         ALLOCATE( qrp(jpi,jpj), erp(jpi,jpj), coefice(jpi,jpj), STAT= sbc_ssr_alloc )
+         !
+         IF( lk_mpp                  )   CALL mpp_sum ( 'sbcssr', sbc_ssr_alloc )
+         IF( sbc_ssr_alloc /= 0 )   CALL ctl_warn('sbc_ssr_alloc: failed to allocate arrays.')
+         !
+      ENDIF
+   END FUNCTION
       
    !!======================================================================
